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

ルジャンドル多項式
1. scipy.special.eval_legendre()
2. scipy.special.legendre()
ルジャンドル陪関数
1. scipy.special.lpmv()

ルジャンドル多項式

ルジャンドルの微分方程式
$\begin{matrix} (1) & \frac{d}{d x} [(1 - x^{2}) \frac{d f (x)}{d x}] + n (n + 1) f (x) = 0 \end{matrix}$
の解 $P_{n} (x)$ を ルジャンドル多項式 (Legendre polynomial) とよび、
$P_{0} (x), P_{1} (x), P_{2} (x), P_{3} (x), \dots$
は区間 [ $- 1, 1$ ] で直交関数系をなします：
$\begin{matrix} (2) & \int_{- 1}^{1} P_{m} (x) P_{n} (x) d x = \frac{2}{2 n + 1} δ_{m n} \end{matrix}$
$P_{n} (x)$ は漸化式
$\begin{matrix} (3) & \begin{aligned} P_{0} (x) = 1 \\ P_{1} (x) = x \\ (n + 1) P_{n + 1} (x) = (2 n + 1) x P_{n} (x) - n P_{n - 1} (x) \end{aligned} \end{matrix}$
から逐次生成することができます。具体的表式を並べると以下のようになります。
$\begin{aligned} P_{0} (x) = 1 \\ P_{1} (x) = x \\ P_{2} (x) = \frac{1}{2} (x^{2} - 1) \\ P_{3} (x) = \frac{1}{2} (x^{3} - 3 x) \\ P_{4} (x) = \frac{1}{8} (x^{4} - 30 x^{2} + 3) \end{aligned}$
ルジャンドル多項式はロドリゲスの公式
$\begin{matrix} (4) & P_{n} (x) = \frac{1}{2^{n} n!} \frac{d^{n}}{d x^{n}} (x^{2} - 1)^{n} \end{matrix}$
から求めることもできます。

scipy.special.eval_legendre()

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

# PYTHON_LEGENDRE_POLYNOMIAL

# In[1]

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

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

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

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

# In[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.special.legendre()

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

# SCIPY_SPECIAL_LEGENDRE

# In[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$ になることを確認してみます。

# In[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$ になっています。

ルジャンドル陪関数

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

scipy.special.lpmv()

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

# PYTHON_ASSOCIATED_LEGENDRE

# In[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)$ をまとめて計算させることができます。

# In[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)$ のグラフをプロットします。

# In[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()