ガウス求積で数値積分を行う-Pythonプログラム

ガウス求積(ガウス則)とは、$2n-1$次の多項式のある閉区間における定積分を$2n$個の値を使って比較的少ない演算で数値的に求める手法です。

1. ガウス求積のアルゴリズム
2. ガウス求積のPythonプログラム

ガウス求積のアルゴリズム

$n$ を正の整数とし、$f(x)$を任意の$2n-1$次の多項式関数とします。ガウス求積では、$f(x)$ の $[−1, 1]$ における定積分値$ I$ を、$n$個の積分点(ガウス点、ガウスノード)$x_i$と$n$個の重み$\omega_i$を用いた総和で正確に表現します。

$$ I = \int^1_{-1} f(x) dx = \sum_{i=1}^n \omega_i f(x_i) $$

上式では積分範囲を$[−1, 1]$ に設定していますが、変数変換により有限な閉区間に持っていくことが可能です。

具体的な導出方法については、下記の記事が詳しいです。

参考URL：
ガウス求積法の導出 – シキノート
(PDF版はこちら)

※ガウス求積の導出方法については上記のサイトが一番わかりやすかったです。丁寧にまとめていただき、とても助かりました。

ここでは、ガウス求積で必要になる積分点$x_i$と重み$\omega_i$の具体的な表式だけを示しておきます。

$P_n(x)$を$n$次のルジャンドル多項式として、積分点$x_i $は$P_n(x)$の零点、すなわち

$$ P_n(x_i) = 0 \ \ (i=1,2,\cdots,n)$$

であり、重み$\omega_i$は下記のように表せます。

$$ \omega_i^2 = \frac{2 (1- x_i^2) }{[n P_{n-1}(x_i)]^2} \ \ (i=1,2,\cdots,n)$$

ガウス求積の式のイメージは、区分求積法の場合には積分を均等な刻み幅の矩形領域に区切って足し合わせるのに対して、ガウス求積では不均等な刻み幅を用いてその総和を求めています。

ガウス求積のPythonプログラム

ガウス求積により区間$[−1, 1]$における$f(x) = x^4$の定積分を実行するPythonコードを以下に記載します。

# -*- coding: utf-8 -*-
"""
ガウス求積(ガウス則)で多項式の[-1,1]における積分を計算する
"""

import sympy as sy
import numpy as np

z = sy.Symbol('z')
f = z**4  # 積分したい関数. 2n-1次以下にする
n = 3  # 積分点数 (整数値)


# ルジャンドル多項式のゼロ点x_iと重みw_iを求める
def leg_weights_roots(n):
    x = sy.Symbol('x')
    roots = sy.Poly(sy.legendre(n, x)).all_roots()  # n次ルジャンドル多項式の根を求める
    x_i = [rt.evalf(20) for rt in roots] 
    w_i = [(2*(1-rt**2)/(n*sy.legendre(n-1, rt))**2).evalf(20) for rt in roots]  # 重み
    return x_i, w_i


# ガウス求積を実行する
def gauss_quadrature(f,n):
    xx, ww = leg_weights_roots(n)  # xx:零点のリスト, ww:重みのリスト
    integ = 0
    for i in range(n):
        integ += ww[i]*f.subs(z,xx[i]).evalf()
    return integ


print(gauss_quadrature(f, n))

実行結果は以下のようになりました。