博士以前

人間です

Gaussian quadrature で遊んでみた

はじめに

最近そこそこ大変な多重積分が必要になり Gaussian quadrature という数値積分法の威力を知りました。日本語で書かれた記事は専門的なものしかなさそうなので、概要をメモしておこうと思います。 なお数学的な厳密さは考慮していないですが、許してください。

基本的な定理

Gaussian quadrature というのは、積分  \int_a^b dx\, f(x) を、被積分関数  f(x) の重み付き和によって近似する手法です。 これの基盤になっているのは多項式積分に関する次の定理です。

 h(x)積分したい 2n-1 次の多項式とする。積分区間 ( [-1,1], [0,\infty]など) 上の、関数  \omega(x) についての直交多項式 p_n(x) とする。つまり  p_n(x)

 \displaystyle{\int_a^b dx\, \omega(x) p_n(x)p_m(x) \propto \delta_{n,m}}

を満たす多項式とする。ただし  n, m多項式の次数を表す。この時次の式が厳密に成立する:

 \displaystyle{\int_a^b dx\, \omega (x) h(x) = \sum_{i=1}^n w_i h(x_i)}

ここで  x_i p_n(x_i)=0 を満たす n 個の点 (ノードと呼ばれる) で、  w_i は各ノードに対応したウェイトと呼ばれる量である。

証明とウェイトの具体的な形は、英語版のウィキペディアに書いてあります。

Gaussian quadrature - Wikipedia

この定理のおかげで、多項式で良く近似できる関数  f(x) については積分が重み付きの和

 \displaystyle{\int_a^b dx\,\omega(x) f(x) \simeq \sum_{i=1}^n w_i f(x_i)}

によって近似できることがわかります。なお考える直交多項式によって、呼び方が微妙に変わります。 例えば  [-1,1] で直交多項式としてルジャンドル多項式を使うときは Gauss–Legendre quadrature と呼ばれ、 [0, \infty] 上でラゲール多項式を使うときは Gauss–Laguerre quadrature と呼ばれたりします。

具体例

具体的に次の3重積分を考えてみましょう*1:

\begin{aligned}
 I &= \int_0^\infty dx \int_{-\infty}^\infty dx_1 \int_{-\infty}^\infty dx_2\, x^3 f(x_1)f(x_2)f(x-x_1-x_2)\,,\\
 &\text{where}\, f(x) = \frac{1}{e^x+1}
\end{aligned}

これを Gauss–Laguerre quadrature で計算するために次のように書き換えます:

\begin{aligned}
 I &= \int_0^\infty dx \int_{0}^\infty dx_1 \int_{0}^\infty dx_2\, \sum_{j_1=\pm 1}\sum_{j_2=\pm 1} x^3 f(j_1x_1)f(j_2x_2)f(x-j_1x_1-j_2x_2)\,,\\
   &= \int_0^\infty dx \omega(x)\int_{0}^\infty dx_1\omega(x_1) \int_{0}^\infty dx_2\omega(x_2)\, 
e^{x+x_1+x_2}
\sum_{j_1=\pm 1}\sum_{j_2=\pm 1} x^3 f(j_1x_1)f(j_2x_2)f(x-j_1x_1-j_2x_2)\,,
\end{aligned}

ここで \omega(x) = e^{-x} としました。この \omega(x) [0,\infty] に対してはラゲール多項式が直交多項式になります。

あとは適当な次数のラゲール多項式を持ってきて、そのノードとウェイトを計算し、公式に当てはめれば積分計算ができます。 ノードとウェイトの計算は次数が上がるほど大変になるのですが、世の中には色々なパッケージが存在するのでありがたく使わせてもらいましょう。

解析解

なお、実はこの積分は解析的に実行することができます;

\begin{aligned}
 I = \frac{457 \pi^6}{5040} = 87.1735836238\dots
\end{aligned}

数値積分

実際に Gauss–Laguerre quadrature を使って上の解析表式を数値的に再現してみます。

コード

以下では Julia の FastGaussQuadrature モジュールを使って計算します。 試しに1次から20次までの quadrature を実行して、結果を比べてみます。

GitHub - ajt60gaibb/FastGaussQuadrature.jl: Gauss quadrature nodes and weights in Julia.

using FastGaussQuadrature

f(x) = 1.0/(exp(x)+1.0)

function integrand(x1, x2, x)
    integrand_tmp = 0.0
    for j1=[-1,1], j2=[-1,1]
        integrand_tmp += x^3 * f(j1*x1) * f(j2*x2) * f(x - j1*x1 - j2*x2)
    end
    return integrand_tmp
end

function main()
    for n in 1:20
        integ_tmp = 0.0
        nodes, weights = gausslaguerre(n)
        
        for i1=1:n, i2=1:n, i=1:n
            x1=nodes[i1]
            x2=nodes[i2]
            x=nodes[i]
            integ_tmp += weights[i1]*weights[i2]*weights[i] * exp(x1+x2+x) * integrand(x1, x2, x)
        end
        
        println("n=", n, ", I=", integ_tmp)
    end
end

main()

実行結果

n=1, I=3.695294445091779
n=2, I=61.88325249650755
n=3, I=87.32157370172439
n=4, I=87.35456973636666
n=5, I=87.22502771278617
n=6, I=87.19838277730469
n=7, I=87.1816826909689
n=8, I=87.1739737777013
n=9, I=87.17318617155999
n=10, I=87.17367651240122
n=11, I=87.17389771350166
n=12, I=87.17387417770796
n=13, I=87.17377359537485
n=14, I=87.17368240366083
n=15, I=87.17362525065808
n=16, I=87.17359768692653
n=17, I=87.17358772170228
n=18, I=87.17358557612305
n=19, I=87.1735857713404
n=20, I=87.17358607174498

20次のラゲール多項式を使うことで7桁の精度で積分値が合うことが確かめられました。 3桁程度の精度でよければ、たった n=5 の計算でも十分なことがわかります。 さらに多重の積分になった時には、台形公式を使うよりずっと効率的であろうと予測できます。*2

まとめ

まあ FORTRAN の QUADPACK やそのラッパーを使えばよいという話だけど、中身をなんとなく理解しておくのは重要だと思います。

*1:見る人が見たらわかるかもしれませんが、これは一応物理的なモチベーションがある式です

*2:本当は台形公式でもやって計算時間を比較したかった