拉格朗日插值法答案

【问题标题】：Lagrange interpolation method拉格朗日插值法
【发布时间】：2012-06-14 08:45:05
【问题描述】：

我使用卷积和 for 循环（循环太多）来计算插值 Lagrange's method，这里是主要代码：

function[p] = lagrange_interpolation(X,Y)
L = zeros(n);
p = zeros(1,n);


% computing L matrice, so that each row i holds the polynom L_i
% Now we compute li(x) for i=0....n  ,and we build the polynomial 

for k=1:n
    multiplier = 1;
    outputConv = ones(1,1);
    for index = 1:n
        if(index ~= k && X(index) ~= X(k))
            outputConv = conv(outputConv,[1,-X(index)]);
            multiplier = multiplier * ((X(k) - X(index))^-1);
        end
    end
    polynimialSize = length(outputConv);
    for index = 1:polynimialSize
        L(k,n - index + 1) = outputConv(polynimialSize - index + 1);
    end
    L(k,:) = multiplier .* L(k,:);
end

% continues 


end

这些对于计算 l_i(x) 的循环来说太多了（这是在最后一次计算 P_n(x) = Sigma of y_i * l_i(x) 之前完成的）。

对使它更正式的 matlab 有什么建议吗？

谢谢

【问题讨论】：

标签： matlab interpolation polynomial-math

【解决方案1】：

是的，有几个建议（在下面的版本 1 中实现）：if 循环可以与上面的for 组合（只需通过下面的jr(jr~=j) 使index 跳过k）； polynomialSize 总是等于 length(outputConv) 总是等于 n （因为你有 n 数据点，(n-1)th 多项式和 n 系数），所以最后一个 for 循环和下一行也可以替换用简单的L(k,:) = multiplier * outputConv;

所以我在http://en.wikipedia.org/wiki/Lagrange_polynomial 上复制了示例（并采用了他们的j-m 表示法，但对我来说j 变为1:n 和m 是1:n 和m~=j），因此我的初始化看起来像

clear; clc;
X=[-9 -4 -1 7]; %example taken from http://en.wikipedia.org/wiki/Lagrange_polynomial
Y=[ 5  2 -2 9];
n=length(X); %Lagrange basis polinomials are (n-1)th order, have n coefficients
lj = zeros(1,n); %storage for numerator of Lagrange basis polyns - each w/ n coeff
Lj = zeros(n); %matrix of Lagrange basis polyns coeffs (lj(x))
L  = zeros(1,n); %the Lagrange polynomial coefficients (L(x))

那么 v 1.0 看起来像

jr=1:n; %j-range: 1<=j<=n
for j=jr %my j is your k
    multiplier = 1;
    outputConv = 1; %numerator of lj(x)
    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j
    for m = mr %my m is your index
        outputConv = conv(outputConv,[1 -X(m)]);
        multiplier = multiplier * ((X(j) - X(m))^-1);
    end
    Lj(j,:) = multiplier * outputConv; %jth Lagrange basis polinomial lj(x)
end
L = Y*Lj; %coefficients of Lagrange polinomial L(x)

如果您意识到 l_j(x) 的分子只是具有特定根的多项式，则可以进一步简化 - 因为 matlab 中有一个很好的命令 - poly。类似地，分母就是在 X(j) 处计算的那个多边形——因为有polyval。因此，v 1.9：

jr=1:n; %j-range: 1<=j<=n
for j=jr
    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j
    lj=poly(X(mr)); %numerator of lj(x)
    mult=1/polyval(lj,X(j)); %denominator of lj(x)
    Lj(j,:) = mult * lj; %jth Lagrange basis polinomial lj(x)
end
L = Y*Lj; %coefficients of Lagrange polinomial L(x)

为什么是 1.9 版而不是 2.0 版？好吧，可能有一种方法可以摆脱最后一个 for 循环，并将其全部写在 1 行中，但我现在想不起来 - 这是 v 2.0 的待办事项 :)

而且，对于甜点，如果您想获得与维基百科相同的图片：

figure(1);clf
x=-10:.1:10;
hold on
plot(x,polyval(Y(1)*Lj(1,:),x),'r','linewidth',2)
plot(x,polyval(Y(2)*Lj(2,:),x),'b','linewidth',2)
plot(x,polyval(Y(3)*Lj(3,:),x),'g','linewidth',2)
plot(x,polyval(Y(4)*Lj(4,:),x),'y','linewidth',2)
plot(x,polyval(L,x),'k','linewidth',2)
plot(X,Y,'ro','linewidth',2,'markersize',10)
hold off
xlim([-10 10])
ylim([-10 10])
set(gca,'XTick',-10:10)
set(gca,'YTick',-10:10)
grid on

生产

享受并随意重复使用/改进

【讨论】：

【解决方案2】：

试试：

X=0:1/20:1; Y=cos(X) and create L and apply polyval(L,1). 
polyval(L,1)=0.917483227909543
cos(1)=0.540302305868140

为什么会有很大的不同？

【讨论】：