Pythonでルジャンドル多項式を使ってガウス積分をする方法

数値解析の分野ではガウス積分という求積手法がよく用いられます。ガウス積分はルジャンドル多項式を使って積分点と重みを算出しますが、初学者はまずこれらの概念を理解するのが難しいです。この記事では簡単な関数を題材にガウス積分を計算する方法をPythonコードと共に紹介します。

こんにちは。wat（@watlablog）です。今回は数値解析の主要技術の1つ、ガウス積分をPythonで書いてみます！

目次（項目クリックでジャンプできます）

動作環境

この記事ではPythonコードを挟みながら説明しますが、以下のPython環境で動作確認をしました。

Python	Python 3.9.6
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.23.0
matplotlib	3.7.1

ガウス積分の概要

ガウス積分の定義式

積分点の求め方：ルジャンドル多項式を利用する

式（1）のガウス積分では積分点の位置を知る必要があります。積分点はルジャンドル多項式（Legendre polynomials）のゼロ点を用います（ルジャンドル多項式については以下の外部サイトがわかりやすいと思います）。

・高校数学の美しい物語：ルジャンドル多項式の性質と計算

百聞は一見にしかずということで、ルジャンドル多項式を可視化してみましょう。以下のPythonコードを実行することで各次数の多項式プロットとゼロ点を可視化できます。

import numpy as np
import matplotlib.pyplot as plt

# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの入れ物を用意する。
fig = plt.figure()
ax1 = fig.add_subplot(111)

# 軸のラベルを設定する。
ax1.set_xlabel('x')
ax1.set_ylabel('y')

# グラフの上下左右に目盛線を付ける。
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')

x = np.linspace(-1, 1, 400)
for n in range(5):  # 最初の4つのルジャンドル多項式をプロット
    # ルジャンドル多項式を求める
    y = np.polynomial.legendre.legval(x, [0]*n + [1])
    ax1.plot(x, y, label=f'P$_{n}$(x)')

    # 交点を求める
    roots = np.polynomial.legendre.legroots([0] * n + [1])
    ax1.scatter(roots, np.zeros_like(roots), color='yellow', edgecolors='black', s=50)
ax1.legend()

# グラフを表示する。
fig.tight_layout()
plt.show()

import numpy as np

import matplotlib.pyplot as plt

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの入れ物を用意する。

fig = plt.figure()

ax1 = fig.add_subplot(111)

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('y')

# グラフの上下左右に目盛線を付ける。

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

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

for n in range(5): # 最初の4つのルジャンドル多項式をプロット

# ルジャンドル多項式を求める

y = np.polynomial.legendre.legval(x, [0]*n + [1])

ax1.plot(x, y, label=f'P

_{n}

(x)')

# 交点を求める

roots = np.polynomial.legendre.legroots([0] * n + [1])

ax1.scatter(roots, np.zeros_like(roots), color='yellow', edgecolors='black', s=50)

ax1.legend()

# グラフを表示する。

fig.tight_layout()

plt.show()

こちらがグラフです。ルジャンドル多項式は[-1, 1]の範囲で定義される正規直交多項式であり、次数に応じて異なる数のゼロ点を持つ特徴があります。

変数変換

実際の被積分関数は必ずしも[-1, 1]の範囲にあるわけではないので、式（2）の関係式を使って変数変換する必要があります。

\begin{matrix} (2) & x = \frac{b - a}{2} ξ + \frac{b + a}{2} \end{matrix}

式（2）は $ξ =$ -1の時 $x = a$ 、 $ξ =$ 1の時 $x = b$ となります。

実際には先にルジャンドル多項式から積分点を計算しておき、式（1）でガウス積分を計算する段階で変数変換を反映させます。[-1, 1]の範囲におけるガウス積分式は式（3）です。

\begin{matrix} (3) & \int_{- 1}^{1} f (ξ) d x \approx \sum_{i = 1}^{n} w_{i} f (ξ_{i}) d x \end{matrix}

式（2）の変換式を考慮すると、積分の演算子は式（4）の関係を持ちます。

\begin{matrix} (4) & d x = \frac{b - a}{2} d ξ \end{matrix}

式（4）を考慮すると変数変換を考慮したガウス積分は式（5）になります。

\begin{array}{rcl} \int_{a}^{b} f (x) d x & = & \int_{- 1}^{1} f (\frac{b - a}{2} ξ + \frac{b + a}{2}) \frac{b - a}{2} d ξ \\ (5) & \approx & \frac{b - a}{2} \sum_{i = 1}^{n} w_{i} f (\frac{b - a}{2} ξ_{i} + \frac{b + a}{2}) d ξ \end{array}

それではこの内容をPythonでコーディングしてみましょう！

Pythonでガウス積分するコード(イメージ図付き)

理論や証明は難しいですが、使うだけなら簡単です。以下に任意の関数をガウス積分するコードを示します。f関数のreturnに関数形状を書きます。ルジャンドル多項式はnp.polynomial.legendre.leggaussを使って求めます。積分区間a, bと次数nを設定して実行してみましょう。

import numpy as np
from matplotlib import pyplot as plt

def f(x):
    """ 被積分関数 """
    return x ** 2


def gaussian_quadrature(f, start, end, n):
    """ ガウス積分をする関数 """

    # 次数nにおけるガウス点と重みを取得する
    xi, weights = np.polynomial.legendre.leggauss(n)
    print(xi)

    # 積分区間の変換をする
    x = 0.5 * (end - start) * xi + 0.5 * (end + start)
    print(x)

    # 積分する
    integral_result = 0.5 * (b - a) * sum(weights * f(x))
    print('Integral value=', integral_result)

    # プロットで可視化する
    plot(x, weights, n)

    return integral_result


def plot(x, weights, n):
    """ グラフ描画する関数 """

    x_val = np.linspace(-0.2, 2.2, 100)
    y_val = f(x_val)

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # 目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # グラフの上下左右に目盛線を付ける。
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    # 軸のラベルを設定する。
    ax1.set_xlabel('x')
    ax1.set_ylabel('f(x)')

    # データをプロットする。
    ax1.plot(x_val, y_val, label='f(x)')
    for i in range(n):
        # 積分点
        ax1.scatter(x[i], f(x[i]), color='yellow', edgecolors='black', s=100, label='Gauss points' if i == 0 else "")

        # 重みが反映された矩形
        width = weights[i] * (b - a) / 2
        ax1.fill_between([x[i] - width / 2, x[i] + width / 2], 0, f(x[i]), color='skyblue', alpha=0.5,
                         label='Weighted area' if i == 0 else "")
    ax1.legend()

    # グラフを表示する。
    fig.tight_layout()
    plt.show()
    plt.close()

    return


if __name__ == '__main__':
    """ Main """

    # 積分区間[a, b]と積分点の数を指定してガウス積分をする
    a, b = 0, 2
    n = 2
    result = gaussian_quadrature(f, a, b, n)

import numpy as np

from matplotlib import pyplot as plt

def f(x):

""" 被積分関数 """

return x ** 2

def gaussian_quadrature(f, start, end, n):

""" ガウス積分をする関数 """

# 次数nにおけるガウス点と重みを取得する

xi, weights = np.polynomial.legendre.leggauss(n)

print(xi)

# 積分区間の変換をする

x = 0.5 * (end - start) * xi + 0.5 * (end + start)

print(x)

# 積分する

integral_result = 0.5 * (b - a) * sum(weights * f(x))

print('Integral value=', integral_result)

# プロットで可視化する

plot(x, weights, n)

return integral_result

def plot(x, weights, n):

""" グラフ描画する関数 """

x_val = np.linspace(-0.2, 2.2, 100)

y_val = f(x_val)

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('f(x)')

# データをプロットする。

ax1.plot(x_val, y_val, label='f(x)')

for i in range(n):

# 積分点

ax1.scatter(x[i], f(x[i]), color='yellow', edgecolors='black', s=100, label='Gauss points' if i == 0 else "")

# 重みが反映された矩形

width = weights[i] * (b - a) / 2

ax1.fill_between([x[i] - width / 2, x[i] + width / 2], 0, f(x[i]), color='skyblue', alpha=0.5,

label='Weighted area' if i == 0 else "")

ax1.legend()

# グラフを表示する。

fig.tight_layout()

plt.show()

plt.close()

return

if __name__ == '__main__':

""" Main """

# 積分区間[a, b]と積分点の数を指定してガウス積分をする

a, b = 0, 2

n = 2

result = gaussian_quadrature(f, a, b, n)

まずはn=2の結果です。黄色点で積分点（ガウス点）、矩形エリアに重みを考慮した面積をプロットしています。

出力ウィンドウには「Integral value= 2.666666666666667」と出ていますが、式（6）から理論値と一致していることが確認できます。

\begin{array}{rcl} \int_{0}^{2} x^{2} d x & = & \frac{x^{3}}{3} |_{0}^{2} \\ (6) & = & \frac{8}{3} = 2.66666 . . . \end{array}

以下は次数違いの積分結果です。ガウス積分の特性から、積分値はほぼ変わりませんが矩形形状の幅寸法が重みによって変わる様子を確認できると思います。

まとめ

本日はここまで。
この記事ではガウス積分を使って積分値を近似計算するための関連式を紹介し、Pythonで計算する方法を紹介しました。

ガウス積分をするためにはルジャンドル多項式を使い積分点と重みを算出する必要があること、そして変数変換が必要であることを学びました。

最後はmatplotlibを使い可視化を行うことで、ガウス積分のイメージを掴むことができました。
今後は多重積分における使用等に進みたいと思います。

また、以下の記事で有限要素法を学び始めたので、有限要素法の計算に取り入れるにはどうしたら良いかも考えていきたいと思います。

・Pythonで基礎から学ぶ1次元梁モデルの有限要素法

難しいと思っていたガウス積分も使うだけなら簡単でした！
Xでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

Pythonでルジャンドル多項式を使ってガウス積分をする方法

動作環境

ガウス積分の概要

ガウス積分の定義式

積分点の求め方：ルジャンドル多項式を利用する

変数変換

Pythonでガウス積分するコード(イメージ図付き)

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

動作環境

ガウス積分の概要

ガウス積分の定義式

積分点の求め方：ルジャンドル多項式を利用する

変数変換

Pythonでガウス積分するコード(イメージ図付き)

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル