この記事をざっくりと

ガウス求積で必要になるルジャンドル多項式の零点と重みを求めるプログラムです。
ルジャンドル多項式は漸化式を用いて再帰的に、零点はニュートン法で求めます。
数学的な部分の解説はありません。(式の導出とかの詳細は私も知りたい…)

環境

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参照)を得ることができました。

参考

ガウス求積-Wikipedia