Quadratic Lagrange Interpolation in MATLAB not working
조회 수: 6 (최근 30일)
이전 댓글 표시
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
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
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!
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Polynomials에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!