Quadratic Lagrange Interpolation in MATLAB not working

Hope you guys can help me here.

I cannot get this code working. I am able to do a Linear Langrange Interpolation, but the script seems to fail when i try to make the Quadratic Lagrange Interpolation.

% Numerical Analysis in Civil Engineering
% Exercise 2b: Quadratic Lagrange polynomial
% -------------------------------------------------------------------------
clear ; close all ; clc ;

%% Input
n = 8 ; % number of intervals in interpolation
ab = [ 0 8 ] ; % limits for the x interval
fx = 'sin(0.5*pi*x)' ; % function to be evaluated

%% Plot the function for visual inspection
x = linspace(ab(1),ab(2),1001) ; % define "continuous" x-axis
f = F(fx,x) ; % "continuous" user-defined function
figure(1) ; plot(x,f,'k-','linewidth',4) ; % plot f(x)
xlabel('x') ; ylabel('f(x)') ; % put labels om the two axes
grid on ; hold on ; % show grid and allow more plotting in the same axes

%% Discrete points
xj = linspace(ab(1),ab(2),n+1) ; % define discrete x-axis
fj = F(fx,xj) ; % discrete function values
plot(xj,fj,'k.','markersize',24)

%% Lagrange polynomial
[p1,x1] = LagrangePoly1(xj,fj,20) ; % linear Lagrange polynomial
[p2,x2] = LagrangePoly2(xj,fj,20) ; % quadratic Lagrange polynomial
plot(x1,p1,'b-','linewidth',2)
plot(x2,p2,'r-','linewidth',2)
legend('f(x)','Discrete points','Lagrange1','Lagrange2','location','best')
hold off

%% Functions

function f = F(fx,x)
f = eval(fx) ;
end

function [P,X] = LagrangePoly1(xj,fj,ns)
% Input:
%   xj : discrete x values
%   fj : discrete function values, fj = f(xj)
%   ns : number of subintervals in each interval
% Output:
%   P  : Lagrange polynomial values, P(X)
%   X  : x values where the polynomial is evaluated

% Number of intervals
n = length(xj)-1 ;

% Loop over intervals
X = [] ; P = [] ; % reset X and P
for j = 1:n
    x = xj(j) + (xj(j+1)-xj(j+0))*linspace(0,1,ns+1) ;
    L0 = (x-xj(j+1))/(xj(j+0)-xj(j+1)) ;
    L1 = (x-xj(j+0))/(xj(j+1)-xj(j+0)) ;
    p = fj(j+0)*L0 + fj(j+1)*L1 ;
    X = [X x] ; P = [P p] ;
end
end

function [P,X] = LagrangePoly2(xj,fj,ns)
% Input:
%   xj : discrete x values
%   fj : discrete function values, fj = f(xj)
%   ns : number of subintervals in each interval
% Output:
%   P  : Lagrange polynomial values, P(X)
%   X  : x values where the polynomial is evaluated

% Number of intervals
n = length(xj)-1 ;

% Loop over intervals
X = [] ; P = [] ; % reset X and P
for j = 1:n
    x = xj(j) + (xj(j+1)-xj(j+0))*linspace(0,1,ns+1) ;
    L0 = ((x-xj(j+1)*(x-xj(j+2))))/(((xj(j+0)-xj(j+1))*(xj(j+0)-xj(j+2)))) ;
    L1 = ((x-xj(j+0)*(x-xj(j+2))))/(((xj(j+1)-xj(j+0))*(xj(j+1)-xj(j+2)))) ;
    L2 = ((x-xj(j+0)*(x-xj(j+1))))/(((xj(j+2)-xj(j+0))*(xj(j+2)-xj(j+1)))) ;
    p = fj(j+0)*L0 + fj(j+1)*L1 + fj(j+2)*L2 ;
    X = [X x] ; P = [P p] ;
end
end

[Jan]

A hint, not solving your problem:

Use an anonymous function instead of a char vector:

% fx = 'sin(0.5*pi*x)'
fx = @(x) sin(0.5 * pi * x);

Answer 1

[Animesh]

It seems you are facing an issue to interpolate using quadratic lagrange interpolation.

Firstly, there seems to be some issue with the creation of lagrange basis polynomials (“L0”, “L1” and “L2”). Here are the modified polynomials that should work.

L0 = ((x-xj(j+1)).*(x-xj(j+2)))./((xj(j)-xj(j+1)).*(xj(j)-xj(j+2))) ; 
L1 = ((x-xj(j)).*(x-xj(j+2)))./((xj(j+1)-xj(j)).*(xj(j+1)-xj(j+2))) ; 
L2 = ((x-xj(j)).*(x-xj(j+1)))./((xj(j+2)-xj(j)).*(xj(j+2)-xj(j+1))) ; 

Also, we need to ensure that the loop in “LagrangePoly2” function does not try to access the element beyond the length of “xj”. Given “n+1” points, we can only form “n-1” quadratic polynomials because each polynomial is based on 3 points. Therefore, the loop should run upto “length(xj)-2” to ensure that we do not exceed the array bounds.

Here is the modified code that should work:

    clear ; close all ; clc ;
    n = 8 ; % number of intervals in interpolation
    ab = [ 0 8 ] ; % limits for the x interval
    fx = 'sin(0.5*pi*x)' ; % function to be evaluated
    %% Plot the function for visual inspection
    x = linspace(ab(1),ab(2),1001) ; % define "continuous" x-axis
    f = F(fx,x) ; % "continuous" user-defined function
    figure(1) ; plot(x,f,'k-','linewidth',4) ; % plot f(x)
    xlabel('x') ; ylabel('f(x)') ; % put labels om the two axes
    grid on ; hold on ; % show grid and allow more plotting in the same axes
    %% Discrete points
    xj = linspace(ab(1),ab(2),n+1) ; % define discrete x-axis
    fj = F(fx,xj) ; % discrete function values
    plot(xj,fj,'k.','markersize',24)
    %% Lagrange polynomial
    [p1,x1] = LagrangePoly1(xj,fj,20) ; % linear Lagrange polynomial
    [p2,x2] = LagrangePoly2(xj,fj,20) ; % quadratic Lagrange polynomial
    plot(x1,p1,'b-','linewidth',2)
    plot(x2,p2,'r-','linewidth',2)
    legend('f(x)','Discrete points','Lagrange1','Lagrange2','location','best')
    hold off

    %% Functions
    function f = F(fx,x)
    f = eval(fx) ;
    end
    function [P,X] = LagrangePoly1(xj,fj,ns)
    n = length(xj)-1 ;
    % Loop over intervals
    X = [] ; P = [] ; % reset X and P
    for j = 1:n
        x = xj(j) + (xj(j+1)-xj(j+0))*linspace(0,1,ns+1) ;
        L0 = (x-xj(j+1))/(xj(j+0)-xj(j+1)) ;
        L1 = (x-xj(j+0))/(xj(j+1)-xj(j+0)) ;
        p = fj(j+0)*L0 + fj(j+1)*L1 ;
        X = [X x] ; P = [P p] ;
    end
    end
    function [P,X] = LagrangePoly2(xj,fj,ns)
    n = length(xj)-2; % Ensure we have two additional points for the quadratic polynomial
    % Loop over intervals using n (length(xj)-2) to avoid going out of bounds
    X = [] ; P = [] ; % reset X and P
    for j = 1:n
        x = xj(j) + (xj(j+2)-xj(j))*linspace(0,1,ns+1) ; % Use j+2 to include the third point
        % Lagrange basis polynomials
        L0 = ((x-xj(j+1)).*(x-xj(j+2)))./((xj(j)-xj(j+1)).*(xj(j)-xj(j+2))) ;
        L1 = ((x-xj(j)).*(x-xj(j+2)))./((xj(j+1)-xj(j)).*(xj(j+1)-xj(j+2))) ;
        L2 = ((x-xj(j)).*(x-xj(j+1)))./((xj(j+2)-xj(j)).*(xj(j+2)-xj(j+1))) ;
        p = fj(j)*L0 + fj(j+1)*L1 + fj(j+2)*L2 ;
        X = [X x] ; P = [P p] ;
    end
    end   
   

Hope this helps!