Legendre多項式的概念以及正交特性在此不多作描述，可以參考數學物理方程相關教材，本文主要討論在數值計算中對於Legendre多項式以及其導數的計算方法。
Legendre多項式的計算
遞推公式
\[\begin{align} (n+1)P_{n+1}(x)=(2n+1)x \cdot P_{n}(x)-nP_{n-1}(x) \qquad (n\ge2) \end{align} \]
通式可以用冪級數表示為以下形式：
\[\begin{align} P_{n}(x)=\sum\limits_{k=0}^{[\frac{n}{2}]}{\frac{(-1)^k(2n-2k)!}{2^nk!(n-k)!(n-2k)!}}x^{n-2k} \qquad n=0,1,2,\cdots \end{align} \]
Legendre多項式前幾項
| \(n\) | \(0\) | \(1\) | \(2\) | \(3\) | \(4\) |
|---|---|---|---|---|---|
| \(P_{n}(x)\) | \(1\) | \(x\) | \(\frac{1}{2}(3x^2-1)\) | \(\frac{1}{2}(5x^3-3x)\) | \(\frac{1}{8}(35x^4-30x^2+3)\) |
在實際數值計算中，按照通項公式計算\(P_n(x)\)會涉及到大量的乘除法運算，同時由於數據字節長度的限制，對於基數較大階乘的運算會導致數據的溢出，因此，在實際的計算中，使用遞推公式計算\(P_n(x)\)更為合適。
%% Legendre多項式Pi(x)
function L=myLegendre(N,x)
% 返回Pi(x),i=0~N ，x為標量
if N<2
error("N is less then 2 !")
end
L=zeros(1,N);%生成Legendre多項式矩陣
L(1,1)=1;
L(1,2)=x;
for i=2:N-1
L(1,i+1)=((2*i-1)*x*L(1,i)-(i-1)*L(1,i-1))/i;
end
end
Legendre多項式的導數
對於\(x\)的一階導數\(P_{n}'(x)\)存在類似的遞推公式
\[\begin{align} (n+1)P'_{n+1}(x)=(2n+1)\{x \cdot P'_{n-1}(x)+P_{n-1}(x)\}-nP'_{n-1}(x) \qquad (n\ge2) \end{align} \]
化簡得到一階導數項的遞推公式與之前形式類似，但包含額外一項\((2n+1)P_{n-1}(x)\)
\[\begin{align} (n+1)P'_{n+1}(x)=(2n+1)x \cdot P'_{n-1}(x)-nP'_{n-1}(x) +(2n+1)P_{n-1}(x) \qquad (n\ge2) \end{align} \]
Legendre多項式前幾項
| \(n\) | \(0\) | \(1\) | \(2\) | \(3\) | \(4\) |
|---|---|---|---|---|---|
| \(P'_{n}(x)\) | \(0\) | \(1\) | \(3x\) | \(\frac{1}{2}(15x^2-3)\) | \(\frac{1}{2}(35x^3-15x)\) |
%% Legendre多項式Pi(x)對於x的一階導數
function L_x=myLegendre_x(N,x)
L=myLegendre(N,x);
L_x=zeros(1,N);%生成Legendre_x多項式矩陣
L_x(1,1)=0;
L_x(1,2)=1;
for i=2:N-1
L_x(1,i+1)=((2*i-1)*x*L_x(1,i)-(i-1)*L_x(1,i-1)+(2*i-1)*L(1,i))/i;
end
end
計算結果
%%demo
x=0.5;N=100;
c=myLegendre(N,x);
c_x=myLegendre_x(N,x);
figure(2)
subplot(2,1,1)
plot(c,'-*');xlabel('N');title("P(x)")
subplot(2,1,2)
plot(c_x,'-o');xlabel('N');title("P'(x)")
