ガウス求積法(Gauss–Legendre 公式)

定積分を近似計算する話題です。問題設定から解説して,後半ではガウス求積法(Gaussian quadrature)を紹介します。直交多項式が出てきておもしろいです。

問題設定

  • abf(x)dxiwif(xi)\displaystyle\int_{a}^bf(x)dx\fallingdotseq\sum_{i}w_if(x_i)
    となるような wi,xiw_i,x_i を求める問題を考えます。

  • つまり,定積分を被積分関数の値 f(xi)f(x_i) の重みつき和で近似する問題を考えます。

  • このような形の近似公式はたくさんあります。最終目標はガウス求積法ですが,まずは簡単な例として台形公式を紹介します。

台形公式

abf(x)dx(ba)2f(a)+(ba)2f(b)\displaystyle\int_a^bf(x)dx\fallingdotseq\dfrac{(b-a)}{2}f(a)+\dfrac{(b-a)}{2}f(b)

特に,f(x)f(x) が1次以下の多項式なら両辺は等しい。

ガウス求積と台形公式

  • 定積分を台形の面積:(上底+下底)×高さ÷2で近似しています。f(x)f(x) が1次関数のときに両者が一致することは図形的に明らかです。

最適な重みを決める

abf(x)dxiwif(xi)\displaystyle\int_{a}^bf(x)dx\fallingdotseq\sum_{i}w_if(x_i)
となる wi,xiw_i,x_i を探しましょう。wiw_ixix_i を同時に探すのは大変なので,まずは xix_i が与えられたときに最適な wiw_i を計算する方法を考えます。

定理1

nn 個の点 x1,...,xn[a,b]x_1,...,x_n\in[a,b] が与えられたとき,

wi=abLi(x)Li(xi)dxw_i=\displaystyle\int_a^b\dfrac{L_i(x)}{L_i(x_i)}dx

とすると,任意の (n1)(n-1) 次以下の多項式 f(x)f(x) に対して
abf(x)dx=i=1nwif(xi)\displaystyle\int_{a}^bf(x)dx=\sum_{i=1}^{n}w_if(x_i)

ただし,Li(x)=ji(xxj)L_i(x)=\displaystyle\prod_{j\neq i}(x-x_j) です。つまり,(xxj)(x-x_j)n1n-1 個の積です。後の証明でも登場しますがラグランジュの補間公式を知っていると Li(x)L_i(x) の意味がわかりやすいです。一般形は複雑なので例を見てみましょう。

例1

n=2,x1=a,x2=bn=2,x_1=a,x_2=b としてみる。

w1=abxabadx=ba2w_1=\displaystyle\int_a^b\dfrac{x-a}{b-a}dx=\dfrac{b-a}{2}
w2=abxbabdx=ba2w_2=\displaystyle\int_a^b\dfrac{x-b}{a-b}dx=\dfrac{b-a}{2}

つまり,定積分の近似式は ba2f(a)+ba2f(b)\dfrac{b-a}{2}f(a)+\dfrac{b-a}{2}f(b) となり台形公式を得る。

より詳しくは,台形公式を用いた積分の近似とその誤差

例2

n=3,a=x1=1,x2=0,b=x3=1n=3,a=x_1=-1,x_2=0,b=x_3=1 としてみる。

w1=11(x0)(x1)(10)(11)dx=13w_1=\displaystyle\int_{-1}^1\dfrac{(x-0)(x-1)}{(-1-0)(-1-1)}dx=\dfrac{1}{3}

同様に w2=43,w3=13w_2=\dfrac{4}{3},w_3=\dfrac{1}{3} となりシンプソンの公式を得る:

11f(x)dx13{f(1)+4f(0)+f(1)}\displaystyle\int_{-1}^1f(x)dx\fallingdotseq\dfrac{1}{3}\{f(-1)+4f(0)+f(1)\}

f(x)f(x) が2次以下なら定理1よりこの近似は正確な等式になります。実はより強く,f(x)f(x) が3次関数でも等号が成立します! より詳しくは,シンプソンの公式の証明と例題

定理1の証明

f(x)f(x)n1n-1 次以下の多項式なら,

f(x)=i=1nLi(x)f(xi)Li(xi)f(x)=\displaystyle\sum_{i=1}^{n}\dfrac{L_i(x)f(x_i)}{L_i(x_i)}

が成立する(※)。この両辺を [a,b][a,b] で積分すると定理1を得る。

※の理由:ラグランジュの補間公式そのものである。n1n-1 次以下の多項式は nn 点での値が分かれば1つに決まる。これを明示的に表すのがラグランジュの補間公式である。

ガウス求積法

次に,最適な xix_i について考えましょう。

以下,a=1,b=1a=-1,b=1 の場合を考えます。積分区間が一般に [a,b][a,b] の場合は y=2xabbay=\dfrac{2x-a-b}{b-a} と変換すれば [1,1][-1,1] の場合に帰着できます。

定理1では「n1n-1 次以下なら近似は正確」でしたが,xix_i をうまく選ぶことで「2n12n-1 次以下なら近似は正確」にパワーアップできます。

定理2:ガウス求積法
  • nn 次のルジャンドル多項式 Pn(x)P_{n}(x) の零点を x1,...,xnx_1,...,x_n とする。
  • 定理1のように重みを決める。つまり wi=abLi(x)Li(xi)dxw_i=\displaystyle\int_a^b\dfrac{L_i(x)}{L_i(x_i)}dx とする。

すると,任意の 2n12n-1 次以下の多項式 f(x)f(x) に対して abf(x)dx=i=1nwif(xi)\displaystyle\int_{a}^bf(x)dx=\sum_{i=1}^{n}w_if(x_i)

ルジャンドル多項式については ルジャンドル多項式の性質と計算 で紹介しています。Pn(x)P_{n}(x)nn 次の多項式で,Pn(x)=0P_{n}(x)=0 は相異なる nn 個の実数解を [1,1][-1,1] 内に持つことが知られています。

n=2n=2 のとき,

  • P2(x)=32x212P_2(x)=\dfrac{3}{2}x^2-\dfrac{1}{2} である。P2(x)=0P_2(x)=0 の解は x1=13,x2=13x_1=-\dfrac{1}{\sqrt{3}},x_2=\dfrac{1}{\sqrt{3}}

  • w1=11x1323dx=1w_1=\displaystyle\int_{-1}^1\dfrac{x-\frac{1}{\sqrt{3}}}{-\frac{2}{\sqrt{3}}}dx=1,同様に w2=1w_2=1

よって,ガウス求積法は

11f(x)dxf(13)+f(13)\displaystyle\int_{-1}^1f(x)dx\fallingdotseq f\left(-\dfrac{1}{\sqrt{3}}\right)+f\left(\dfrac{1}{\sqrt{3}}\right)

となる。ff が3次以下ならこの近似式は正確。

定理2の証明

f(x)f(x)2n12n-1 次以下の場合を考える。f(x)f(x)Pn(x)P_n(x) で割った商を Q(x)Q(x),余りを R(x)R(x) とおく。

定積分の値は,
11f(x)dx=11Pn(x)Q(x)dx+11R(x)dx=11R(x)dx\displaystyle\int_{-1}^1f(x)dx\\ =\displaystyle\int_{-1}^1P_n(x)Q(x)dx+\displaystyle\int_{-1}^1R(x)dx\\ =\displaystyle\int_{-1}^1R(x)dx

ただし,最後の等号はルジャンドル多項式の直交性を用いた(nn 次のルジャンドル多項式は n1n-1 次以下の多項式と直交する)。

一方,定理1の証明と同じようにラグランジュの補間公式:
R(x)=i=1nLi(x)R(xi)Li(xi)R(x)=\displaystyle\sum_{i=1}^{n}\dfrac{L_i(x)R(x_i)}{L_i(x_i)} の両辺を積分すると,

11R(x)dx=i=1nwiR(xi)\displaystyle\int_{-1}^1R(x)dx=\sum_{i=1}^n w_iR(x_i)

さらに,Pn(xi)=0P_n(x_i)=0 を使うと右辺は

i=1nwiR(xi)=i=1nwi{R(xi)+Pn(xi)Q(xi)}=i=1nwif(xi)\displaystyle\sum_{i=1}^n w_iR(x_i)\\ =\displaystyle\sum_{i=1}^n w_i\{R(x_i)+P_n(x_i)Q(x_i)\}\\ =\displaystyle\sum_{i=1}^n w_if(x_i)

となり定理2は示された。

補足

  • ルジャンドル多項式以外の直交多項式についても,似たような公式が知られています(重み関数がついたり,積分区間が [0,][0,\infty] になるなど少し形は変わります)。

  • ガウス求積法のことをガウスの数値積分公式,ガウス・ルジャンドル (Gauss–Legendre) 公式などとも呼びます。

  • 定理2において,wiw_i を変形していくと,wi=2(1xi2)(Pn(xi))2w_i=\dfrac{2}{(1-x_i^2)(P_n'(x_i))^2} となることが知られています。ルジャンドル多項式の微分が出てきます。

知人のH氏に Inspire されて書いた記事です。