是的,一些建议(在下面的版本1中实现)：if循环可以与上面的组合(只需通过下面的jr(jr~ = j)使索引跳过k); polynomialSize总是等长(outputConv),它总是等于n(因为你有n个数据点,第n-1个多项式有n个系数),所以最后一个for循环和下一个行也可以用简单的L(k, 🙂 =乘数* outputConv;
所以我在/wiki/Lagrange_polynomial上复制了这个例子(并采用了他们的j-m表示法,但对我来说j是1：n而m是1：n和m~ = j),因此我的初始化看起来像
clear; clc;
X=[-9 -4 -1 7]; %example taken from /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)处评估的多边形 – 因为存在多边形.因此,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？好吧,有可能有办法摆脱这最后的循环,并将其全部写入1行,但我现在想不到它 – 它是v 2.0的todo 🙂
而且,对于甜点,如果你想获得与维基百科相同的图片：
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
产生
享受并随意重用/改进
