ええ、いくつかの提案 (以下のバージョン 1 で実装):if
ループはその上と組み合わせることができます(以下のような方法でスキップするfor
だけです)。常に等しい(データポイント、係数のある多項式があるため)常に等しいため、最後のループと次の行も単純に置き換えることができますindex
k
jr(jr~=j)
polynomialSize
length(outputConv)
n
n
(n-1)th
n
for
L(k,:) = multiplier * outputConv;
だから私はhttp://en.wikipedia.org/wiki/Lagrange_polynomialの例を複製しました(そしてそれらのj-m
表記法を採用しましたが、私にとってj
は go 1:n
and m
is 1:n
and 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) で評価された polyn です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 の 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
生産する
楽しんで、自由に再利用/改善してください