ルジャンドル多項式・陪関数

ルジャンドル多項式・陪関数

ルジャンドル多項式

 ルジャンドルの微分方程式
 
\[\frac{d}{dx}\left[(1-x^2)\frac{df(x)}{dx}\right]+n(n+1)f(x)=0\tag{1}\]
の解 $P_n(x)$ を ルジャンドル多項式 (Legendre polynomial) とよび、
 
\[P_0(x),\ P_1(x),\ P_2(x),\ P_3(x),\ ...\]
は区間 [$-1,\ 1$] で直交関数系をなします:
 
\[\int_{-1}^{1}P_m(x)P_n(x)dx=\frac{2}{2n+1}\delta_{mn}\tag{2}\]
 $P_n(x)$ は漸化式
 
\[\begin{align*}
&P_0(x)=1\\[6pt]
&P_1(x)=x\\[6pt]
&(n+1)P_{n+1}(x)=(2n+1)xP_n(x)-nP_{n-1}(x)\end{align*}\tag{3}\]
から逐次生成することができます。具体的表式を並べると以下のようになります。
 
\[\begin{align*}
&P_0(x)=1\\[6pt]
&P_1(x)=x\\[6pt]
&P_2(x)=\frac{1}{2}(x^2-1)\\[6pt]
&P_3(x)=\frac{1}{2}(x^3-3x)\\[6pt]
&P_4(x)=\frac{1}{8}(x^4-30x^2+3)\end{align*}\]

 ルジャンドル多項式はロドリゲスの公式
 
\[P_n(x)=\frac{1}{2^nn!}\frac{d^n}{dx^n}(x^2-1)^n\tag{4}\]
から求めることもできます。
 

scipy.special.eval_legendre()

  scipy.special.eval_legendre(n, x) は n 次のルジャンドル多項式の点 x における値を評価します。

# scipy.special.eval_legendre-1

x = np.linspace(-1, 1, 5)

# 4次のルジャンドル多項式のxにおける値を計算
p = eval_legendre(4, x)

print(p)
[ 1.        -0.2890625  0.375     -0.2890625  1.       ]

 以下のコードを実行すると、ルジャンドル多項式のグラフを描画します。

# scipy.special.eval_legendre-2
# ルジャンドル多項式のグラフ

# FigureとAxesを設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_title("Legendre polynomials", size = 15)
ax.grid()
ax.set_xlim(-1, 1)
ax.set_ylim(-1.25, 1.25)
ax.set_xlabel("x", size = 14, labelpad = 8)
ax.set_ylabel("Pn(x)", size = 14, labelpad = 8)

# x座標データ
x = np.linspace(-1, 1, 129)

# 色のリスト
c = ["red", "blue", "green", "darkorange", "purple"]

# P0,P1,P2,P3,P4のグラフを描画
for i in range(5):
    ax.plot(x, eval_legendre(i, x),
            label = "n = {}".format(i), color = c[i-1])

# 凡例を表示
ax.legend()

plt.show()

 SciPy ルジャンドル多項式 (Legendre polynomial) の計算
 

scipy.special.legendre()

 scipy.special.legendre(n) はルジャンドル多項式オブジェクトを生成します。

# scipy.special.legendre-1

# 4次のルジャンドル多項式を生成
p4 = legendre(4)

print(p4)
       4             3        2
4.375 x + 4.857e-16 x - 3.75 x + 2.429e-16 x + 0.375

 多項式オブジェクトを積分する integ()メソッドを使って、$P_3(x)P_5(x)$ の [$-1,\ 1$] にわたる積分が $0$ になることを確認してみます。

# scipy.special.legendre-2

# 3次のルジャンドル多項式
p3 = legendre(3)

# 5次のルジャンドル多項式
p5 = legendre(5)

# p3*p5の不定積分
p3p5_i = (p3 * p5).integ()

# p3*p5の[-1,1]にわたる積分
p3p5_id = p3p5_i(1)-p3p5_i(-1)

print(p3p5_id)
1.1102230246251565e-15

 近似計算なので若干の誤差はありますが、ほぼ $0$ になっています。
 

現場で使える Django の教科書《基礎編》

新品価格
¥1,728から
(2019/8/21 23:36時点)

ルジャンドル陪関数

 ルジャンドルの陪微分方程式
 
\[\left[(1-x^2)\frac{d^2}{dx^2}-2x\frac{d}{dx}+n(n+1)-\frac{m^2}{1-x^2}\right]f(x)=0\tag{5}\]
の解を ルジャンドル陪関数 (Associated Legendre polynomials) とよびます:
 
\[P_n^{m}(x)=\frac{1}{2^mm!}(1-x^2)^{m/2}\frac{d^{m+n}}{dx^{m+n}}(x^2-1)^n\tag{6}\]
 右辺の $(x^2-1)^n$ が $2n$ 次多項式なので、$2n+1$ 回の微分で $0$ になります。
 したがって、ルジャンドル陪関数は
 
\[-n\leq m\leq n\]
の範囲で意味をもちます。いくつかの $n$ と $m$ について $P_n^{m}(x)$ を書き並べてみます。
 
\[\begin{align*}
&P_0^{0}(x)=1\\[6pt]
&P_1^{0}(x)=x\\[6pt]
&P_1^{1}(x)=\sqrt{1-x^2}\\[6pt]
&P_2^{0}(x)=\frac{1}{2}(3x^2-1)\\[6pt]
&P_2^{1}(x)=3x\sqrt{1-x^2}\\[6pt]
&P_2^{2}(x)=3(1-x^2)\end{align*}\]

scipy.special.lpmv()

 scipy.special.lpmv(x1, x2, x3) はルジャンドル陪関数の値を返します。
 引数を (6) の表式に合わせると lpmv(m, n, x) です。

# scipy.special.lpmv-1

# P(m=1,n=2,x=0.5)を計算
p = lpmv(1,2,0.5)

print(p)
-1.299038105676658

 すべての引数に配列を渡すことができます。たとえば、$m=0,\ x=0.5$ に固定して、$P_{0}^{0}(0.5),\ P_{1}^{0}(0.5),\ P_{2}^{0}(0.5)$ をまとめて計算させることができます。

# scipy.special.lpmv-2

n = [0, 1, 2]
p = lpmv(0, n, 0.5)

print(p)
[ 1.     0.5   -0.125]

 以下のコードを実行すると、$P_{3}^{0}(x),\ P_{3}^{1}(x),\ P_{3}^{2}(x),\ P_{3}^{3}(x)$ のグラフをプロットします。

# scipy.special.lpmv-3
# ルジャンドル陪関数のグラフ

# FigureとAxesを設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_title("Associated Legendre polynomials", size = 15)
ax.grid()
ax.set_xlim(-1, 1)
ax.set_ylim(-4, 4)
ax.set_xlabel("x", size = 14, labelpad = 8)
ax.set_ylabel("$P_n^m(x)$", size = 14, labelpad = 8)

# x座標データ
x = np.linspace(-1, 1, 129)

# 色のリスト
c = ["red", "blue", "green", "darkorange"]

# ルジャンドル陪関数のグラフをプロット
for i in range(4):
    ax.plot(x, lpmv(i, j, x),
            label = "m = {},  n = 3".format(i), color = c[i])

# 凡例を表示
ax.legend()

plt.show()

 Scipy ルジャンドル陪関数 (Associated Legendre polynomials)