この記事をざっくりと
ガウス求積で必要になるルジャンドル多項式の零点と重みを求めるプログラムです。
ルジャンドル多項式は漸化式を用いて再帰的に、零点はニュートン法で求めます。
数学的な部分の解説はありません。(式の導出とかの詳細は私も知りたい…)
環境
GoogleColaboratory
ライブラリはnumpyとnumba jitを使いました。
import numpy as np
from numba import jit
ルジャンドル多項式
ルジャンドル多項式は
$$
P_n(x) = \frac{(-1)^n}{2^nn!}\frac{d^n}{dx^n}(1-x^2)^n
$$で表される関数のようです。
これを解くのは骨が折れそうなので、今回は漸化式を利用します。
ルジャンドル多項式は以下の漸化式を満たします。
\begin{align}
P_0(x) &= 1\\
P_1(x) &= x\\
P_n(x) &= \frac{(2n-1)xP_{n-1}(x)-(n-1)xP_{n-2}(x)}{n}
\end{align}
これをPythonで実装してみました。
@jit
def legendre(n,x):
if n == 0:
return 1
elif n == 1:
return x
else:
Pox = (x*(2*n-1)*legendre(n-1,x)-(n-1)*legendre(n-2,x))/(n)
return Pox
こんな感じ。
※関数を再帰的に呼び出しているので、$n$が大きくなるとめちゃめちゃ時間かかります。jitを使うとだいぶ早くなります。
零点を求める
ニュートン法で零点を求めます。
ニュートン法は、わかりやすい解説があちらこちらにあるので詳しいことは割愛します。
以下の式
x_{i+1} = x_i - \frac{P_n(x_i)}{P'_n(x_i)}
を、$ x_{i+1} - x_i $が十分小さくなるまで繰り返し解くことで零点が求まります。
なお、$x_{i+1}$を求めるのに$P_n'(x_i)$が必要です。ルジャンドル多項式の場合、どうやら
$$
P_n'(x)=\frac{n(P_{n-1}(x)-xP_n(x))}{1-x^2}
$$で求められるようです。
初期値は、
$$
x_i = \cos\biggl(\frac{i+0.75}{n+0.5}\pi\biggr), (i=0,1,2,...,n-1)
$$を用いると良いようです。
プログラムは以下のようにしました。
#ルジャンドル多項式の微分
@jit
def diff_legendre(n,x):
if (x == 1) or (x == -1):
x += 0.001
return n*(legendre(n-1,x) - x*legendre(n,x))/(1-x**2)
#ニュートン法で零点を一点求める
def ZeroPont(x0, n, N=50): #Nは零点を探すための最大試行回数
for i in range(N):
x1 = x0 - legendre(n,x0)/diff_legendre(n,x0)
if abs(x1 - x0) < 1e-6: #差が十分に小さいことを評価
break
x0 = x1
if i == 50:
print("zero point not found")
return x1
#Pn(x)の全ての零点を取得する
def search_ZeroPoint_list(n):
zero_point_list = []
for i in range(n):
x = np.cos(((i+0.75)/(n+0.5))*np.pi)#初期値
ZP = ZeroPont(x, n)
if abs(ZP) < 1e-10:
ZP = 0.0
return zero_point_list
重みを求める
分点を零点とすると、ルジャンドル多項式の重みは
$$
w_i = \frac{2(1-x_i^2)}{(nP_n(x_i))^2}
$$で求められるようです。
#重みを計算する
def weight(x,n):
return 2*(1-(x**2)) / ((n*legendre(n-1, x))**2)
作成したsearch_ZeroPoint_list関数を少し変えて、零点と同時に重みを計算できるようにしました。
#零点と重みを計算する
def search_ZeroPoint_and_Weight(n):
ZP_list = []
weight_list = []
for i in range(n):
x = np.cos(((i+0.75)/(n+0.5))*np.pi)#初期値
ZP = ZeroPont(x, n)
if abs(ZP) < 1e-10:
ZP = 0.0
W = weight(ZP,n)
ZP_list.append(x)
weight_list.append(w)
return ZP_list, weight_list
値を出してみる
for i in range(4):
zps, ws = search_ZeroPoint_and_Weight(i)
print("n = " , str(i) , "\nzero_points = " , zps, "\nweights = " , ws, "\n" )
n = 0
zero_points = []
weights = []
n = 1
zero_points = [0.0]
weights = [2.0]
n = 2
zero_points = [0.5773502691896257, -0.5773502691896257]
weights = [1.0000000000000004, 1.0000000000000004]
n = 3
zero_points = [0.7745966692414841, 0.0, -0.7745966692414841]
weights = [0.5555555555555541, 0.8888888888888888, 0.5555555555555541]
解析解に近い値になっていますね。
まとめ
漸化式を再帰関数で再現することで、ルジャンドル多項式を実装しました。
零点はニュートン法により求めました。特に工夫せずとも、解析解に近い値(Wikipedia参照)を得ることができました。