3

次のように定義されたルジャンドル多項式をプロットしようとしています。

P0(x) = 1
P1(x) = x
Pn+1(x) = ((2n+1)/(n+1)) * x * Pn(x) - (n / (n+1)) * Pn-1(x)

私はそれを簡単でゆっくりとした方法で行い、直接的でもう少し複雑な方法で行いました。どちらも同様の数値になりますが、まったく同じではありません。振幅が異なります。以下は、敬意を表したコードです (定義の添字 n+1 を n に調整していることに注意してください)。

xi = linspace(-1,1,500);
n = 10;

方法 1:

Pn = cell(n+1,1);
Pn{1} = @ (x) 1;
Pn{2} = @ (x) x;
for i=3:(n+1)
    Pn{i} = @ (x) ((2*(i-1)+1)/(i)).*x.*Pn{i-1}(x) - ((i-1)/i) .* Pn{i-2}(x);
end
plot(xi,Pn{1}(xi),'--r',xi,Pn{2}(xi),'--g',xi,Pn{3}(xi),'--b',...
    xi,Pn{4}(xi),'--m',xi,Pn{5}(xi),'--c',xi,Pn{6}(xi),'--k');
legend('P0','P1','P2','P3','P4','P5');

図1：

方法 1

方法 2:

%Notice here that the results of Pj get stored into YI(j+1)
YI = zeros(length(xi),6);
YI(:,1) = ones(size(xi))';
YI(:,2) = xi';
for i=3:6;
    Pn1 = 1;
    Pn2 = xi;
    for j=2:(i-1) 
        Pn3 = ((2*(j-1)+1)/j) .* xi .* Pn2 - ((j-1) / j) .* Pn1;
        Pn1 = Pn2;
        Pn2 = Pn3;
    end
    YI(:,i) = Pn3';
end
figure('Name','direct method');
plot(xi,YI(:,1)','--r',xi, YI(:,2)', '--g', xi, YI(:,3)', '--b', ...
    xi,YI(:,4)','--m', xi,YI(:,5)', '--c', xi,YI(:,6)', '--k');

図 2:

方法 2

控えめに言っても、これは奇妙です。方法 1 では、P11 までのすべての多項式を計算していますが、最初の 6 つだけを使用してプロットします。何が起こっているのか誰か知っていますか？

4

1 に答える 1

2

方法 2 はもっと簡単にできます。

X = linspace(-1,1,500); X = X(:);
N = 10;

Y = zeros(numel(X),N);
Y(:,1) = 1;
Y(:,2) = X;
for n=1:(N-1)
    Y(:,n+2) = ((2*n+1) .* X .* Y(:,n+1) - n .* Y(:,n)) / (n+1);
end

figure, plot(X, Y(:,1:6))
legend(num2str((1:6)'-1,'P_%d(x)'))
xlabel('x'), ylabel('P_n(x)'), title('Legendre Polynomials')

legendre_polynomials

これは、ウィキペディアのページに示されているプロットと同等です。

編集：

配列のインデックス付けに 1 つずつずれているバグがありました。MATLAB は 1 から始まるインデックスを使用しますが、式は 0 から始まる方法で定義されます。現在は修正されています。混乱して申し訳ありません;)

方式

確認するP(n=2,x=0)には-1/2：

>> interp1(X, Y(:,3), 0)
ans =
   -0.5000
于 2014-06-06T09:32:55.920 に答える