Quadratic Lagrange Interpolation in MATLAB not working

조회 수: 6 (최근 30일)
Tobias Frederiksen
Tobias Frederiksen 2021년 10월 4일
답변: Animesh 2024년 2월 28일
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
  댓글 수: 1
Jan
Jan 2021년 10월 4일
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);

댓글을 달려면 로그인하십시오.

답변 (1개)

Animesh
Animesh 2024년 2월 28일
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!

카테고리

Help CenterFile Exchange에서 Polynomials에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by