9. 数値積分 (1) 数値積分の基本公式

有限区間あるいは無限区間における定積分の近似計算を数値積分という。代表的な基本公式について文献[1]の例題を参考に数値実験を行う。

数値積分では、積分区間 [a, b] を n分割する適当な分点列 a = x0 < x1 < ⋯ < xn = b をとり、その重み付きの和を積分近似値とする。

    ∫ f(x)dx [a, b] = Σwif(xi) (i=0~n)

あるいは、分点が積分区間の端点に一致しない場合には、N個の分点 a < x1 < x2 < ⋯ < xN < b により次のように表す。

    ∫ f(x)dx [a, b] = Σwif(xi) (i=1~N)

ここで、wi は xi における重みである。

(1) 補間型積分公式

被積分関数 f(x) をLagrangeの補間多項式 pn(x) で近似し、その積分値が被積分関数の積分値に等しくなるように重み wi を決めた公式を補間型積分公式という。

すなわち、

    pn(x) = Σli(x)f(xi) (i=0~n)
    ただし、li(x) = Π(x-xj)/(xi-xj) (j=0~n, j≠i)

として、これを積分すると、

    ∫ pn(x)dx [a, b]
    = Σ(∫ li(x)dx [a, b])f(xi) (i=0~n)

となるから、

    wi = ∫ li(x)dx [a, b]

として重みを決めることができる。

補間型積分公式では、重みと分点のとり方により多くの公式が導かれる。

公式 nまたはN 重みwi 分点xi 参考
台形則 n=1 w0 = w1 = (b – a)/2 x0=a, x1=b Newton-Cotes公式で n=1 とした場合に一致する
Simpson則 n=2 w0 = w2 = h/3
w1 = 4h/3
x0=a, x1=a+h, x2=a+2h=b (等間隔 (h=(b – a)/2)) Newton-Cotes公式で n=2 とした場合に一致する
Newton-Cotes公式 n=1~7, 9 h×有理数で表される
(n=8 および n≥10 では負の値が現れる)
x0=a, x1=a+h, … , xn=b
(等間隔 (h=xi+1– xi))
nが奇数のときはn次、nが偶数のときはn+1次までの多項式について正確な積分値を与える
n=3 の場合をSimpsonの3/8則とも呼ぶ
重みに負の値が現れると桁落ちを引き起こす原因になるので高次の(n=8 および n≥10の)公式は推奨されない
Chebyshev公式 N=2~7, 9 一定 (になるように分点を選ぶ) a < x1 < ⋯ < xN < b
不等間隔 (値は無理数になる)
N=8 および N≥10 では分点の値が複素数になるため通常は使われない
Gauss型公式 N=2~ 積分区間上の直交多項式のゼロ点を分点に選び、それから重みを決める (重みは一定ではなく、分点は不等間隔: a < x1 < ⋯ < xN < b)
直交多項式の種類により、Gauss-Legendre公式(*)、Gauss-Laguerre公式、Gauss-Hermite公式など多くの種類がある
(*) 単にGauss公式というとGauss-Legendre公式を指す
直交多項式のゼロ点を分点に選ぶことにより、2N-1次までの多項式について正確な積分値を与える

補間型積分公式では、Newton-Cotes公式およびChebyshev公式に比べ、Gauss型公式がおおよそ倍の精度を持つことがわかる。

それぞれの公式の重みおよび分点の具体的な値を下の付表に示す(符号を除いて対称であるから実際には半分だけ記憶しておけばよい)。

(2) 複合公式

積分区間を分割し、その小区間ごとに公式を適用するものを複合公式という。 例えば、積分区間をm分割し台形則を適用すると、

    Im = h (f(x0)/2 + f(x1) + f(x2) + ・・・ + f(xm)/2), h = xi+1- xi (i=0~m)

となり、これを複合台形則というが、単に台形則ということが多い。同様に、複合Simpson則は単にSimpson則と呼ばれることが多い。台形則やSimpson則では、実用上、複合公式を使うのが普通である。

Newton-Cotes公式においても、高次の公式では桁落ちにより精度が悪くなる恐れがあるため積分区間を分割して低次の公式を適用するのが一般的である。

(3) Romberg積分

数列ak が定数aに収束するものとして、収束率λを

    λ = lim (ak+1 - a)/(ak - a)  [k→∞]

と定義する。λの値が事前にわかっているとして、これをaについて解いて求めた値

    bk = (ak+1 - λak)/(1 - λ)

による数列{bk}は、元の数列{ak}よりも速く収束すると予想される。これをRichardsonの加速法という。

(複合)台形則において分点数を倍々にした系列 I1, I2, ⋯ を作れば、hは半々になってゆき、誤差は1段進むごとに1/4になる(収束率が1/4)。 これにRichardsonの加速法を適用した補外を行い、

    Ik(1) = (Ik+1 - λIk)/(1 - λ) = (4Ik+1 - Ik)/(4 - 1)

とした系列 {I1(1), I2(1), ⋯ } を得ることができる。これはSimpson則に一致し、収束率λ(1)は1/42になる。更に収束を加速するために、新しい系列に第2回補外を行うと

    Ik(2) = (Ik+1(1) - λ(1)Ik(1))/(1 - λ(1)) = (42Ik+1(1) - Ik(1))/(42 - 1)

とした系列 {I1(2), I2(2), ⋯ } を得る。このように、補外を繰り返して収束を加速する方法をRomberg積分法という。

(4) 数値実験

性質の良い関数の積分

次の積分を考える。

    I = ∫ ex cos x dx [0, 1] = 1.37802 46135 47364

これは性質の良い(解析的な)関数で、グラフに表すと次のような形をしている。

ex1f

これを種々の積分公式で分点数N(横軸)を変えて計算し、その時の誤差(縦軸)を対数グラフにすると次のようになった。

ex1

台形則は小区間数をN-1とした複合公式の結果を示した。また、Simpson則は小区間数を(N-1)/2とした複合公式の結果を示した。Newton-Cotes公式では、N=n+1 としてプロットした。

Chebyshev公式(N≤7およびN=9のみ計算)はNewton-Cotes公式とほぼ同等の結果であるが、不等間隔分点が必要なこと、N≤7と N=9 以外では係数が複素数になることなど取扱いがやや不便なため使われることは少ない。

Newton-Cotes公式(n≤10のみ計算)はGauss公式の半分程度の精度であるが、等間隔分点というメリットがある。被積分関数値が任意の点では計算できず、表形式で等間隔に与えられているような場合にも使える。また、適応自動積分に適用しやすく、しばしば使われている。

Gauss公式は他の公式に比べ精度がよく N=7 で最高精度(15桁)に達している。係数の記憶量が他の公式より多くなるが、できればこれを使用するとよい。

台形則とSimpson則は収束が遅くこのままで実用性に乏しい。

分点を倍々に増やしていったときの台形則の計算結果を見てみると次のようになった。

分点数 積分値(計算値) 積分値(真値) 誤差 収束率
2 1.34061800327106 1.37802461354736 0.037406610
4 1.36858238253106 1.37802461354736 0.009442231 1/3.961
8 1.37565843490021 1.37802461354736 0.002366179 1/3.990
16 1.37743271822098 1.37802461354736 0.000591895 1/3.997
32 1.37787661780930 1.37802461354736 0.000147996 1/3.999

収束率はほぼ1/4になっており、Romberg積分の項で示した仮定が成りたっていることがわかる。図でもRomberg積分は台形則に比べて明らかに収束が加速しており、期待通りの動きをしている。ただし、実用性は他の公式に比べもの足りない。

特異点の影響

積分区間外に特異点を持つ次の関数の積分を考える。

    I = ∫ 1/(1+x)dx [0, 4] = ln(1+4) = ln(5)

次のように積分区間 [0, 4] では滑らかな関数に見えるが、近くの x=-1 に特異点を持つ。

ex2f

これを上と同じように種々の積分公式で計算し誤差を対数グラフにする。

ex2

図のように、どの公式でも精度が落ちた。Gauss公式でも N=18 以上でないと最高精度は得られない。

次に以下の積分を考える。

    I = ∫ 1/(1+x2)dx [0, 4] = arctan(4)

これは、実数軸上には特異点がないが、複素平面では近く(x=i)に特異点がある。

ex3f

ex3

これも一つ前の例と同じような傾向を示した。Romberg積分は台形則に対して改善になっていない。

次に端点に特異点があるとどうなるか見てみる。

    I = ∫ √(1-x2)dx [0, 1] = π/4

x=1が特異点である。

ex4f

ex4

このように端点とはいえ区間内に特異点がある場合にはどの公式でもまともな積分値は得られなかった。

この例における台形則の計算結果を見てみると次のようになった。

分点数 積分値(台形則の計算値) 誤差 収束率
2 0.683012701892219 1.0238546e-001
4 0.748927267025610 3.6470896e-002 1/2.807
8 0.772454786089293 1.2943377e-002 1/2.818
16 0.780813259456935 4.5849039e-003 1/2.823
32 0.783775605719283 1.6225577e-003 1/2.826

収束率は1/4ではなく1/2.8にしかなっておらず、正常な動きをしていないことを示唆している。Romberg積分の効果も小さい。

付表. 積分公式の係数

Newton-Cotes公式

    ∫ f(x)dx [a, b] = hΣaif(xi) (i=0~n)
    h = (b - a)/n

において、係数 ai を下表に示す。

 n  a0 a1  a2 a3  a4 a5  a6
 1  1/2 1/2
 2  1/3 4/3 1/3
 3  3/8 9/8 9/8 3/8
 4  14/45 64/45 24/45 64/45 14/45
 5  95/288 375/288 250/288 250/288 375/288 95/288
 6  41/140 216/140 27/140 272/140 27/140 216/140 41/140
 7  5257/17280 25039/17280 9261/17280 20923/17280 20923/17280 9261/17280 25039/17280
 8  3956/14175 23552/14175 -3712/14175 41984/14175 -18160/14175 41984/14175 -3712/14175
 9  25713/89600 141669/89600 9720/89600 174096/89600 52002/89600 52002/89600 174096/89600
n a7 a8 a9
7 5257/17280
8 23552/14175 3956/14175
9 9720/89600 141669/89600 25713/89600

Chebyshev公式

    ∫ f(x)dx [a, b] = (b - a)/N Σf((a + b)/2 + (b - a)ui/2) (i=1~N)

において、係数 ui を下表に示す。

N u1 u2 u3 u4 u5 u6
2 -0.577350269189626 0.577350269189626
3 -0.707106781186547 0 0.707106781186547
4 -0.794654472291766 -0.187592474085080 0.187592474085080 0.794654472291766
5 -0.832497487000982 -0.374541409553581 0 0.374541409553581 0.832497487000982
6 -0.866246818107821 -0.422518653761112 -0.266635401516705 0.266635401516705 0.422518653761112 0.866246818107821
7 -0.883861700758049 -0.529656775285157 -0.323911810519908 0 0.323911810519908 0.529656775285157
9 -0.911589307728434 -0.601018655380238 -0.528761783057880 -0.167906184214804 0 0.167906184214804
N u7 u8 u9
7 0.883861700758049
9 0.528761783057880 0.601018655380238 0.911589307728434

Gauss公式

    ∫ f(x)dx [a, b] = (b - a)/2 Σwif((a + b)/2 + (b - a)ui/2) (i=1~N)

において、係数 wiおよび ui を下表に示す。

N w1 w2 w3 w4 w5 w6
2 1 1
3 0.555555555555556 0.888888888888889 0.555555555555556
4 0.347854845137454 0.652145154862546 0.652145154862546 0.347854845137454
5 0.236926885056189 0.478628670499367 0.568888888888889 0.478628670499367 0.236926885056189
6 0.171324492379170 0.360761573048139 0.467913934572691 0.467913934572691 0.360761573048139 0.171324492379170
7 0.129484966168870 0.279705391489277 0.381830050505119 0.417959183673469 0.381830050505119 0.279705391489277
8 0.101228536290376 0.222381034453374 0.313706645877887 0.362683783378362 0.362683783378362 0.313706645877887
9 0.0812743883615744 0.180648160694857 0.260610696402935 0.312347077040003 0.330239355001260 0.312347077040003
10 0.0666713443086881 0.149451349150580 0.219086362515982 0.269266719309996 0.295524224714753 0.295524224714753
N w7 w8 w9 w10
7 0.129484966168870
8 0.222381034453374 0.101228536290376
9 0.260610696402935 0.180648160694857 0.0812743883615744
10 0.269266719309996 0.219086362515982 0.149451349150580 0.0666713443086881
N u1 u2 u3 u4 u5 u6
2 -0.577350269189626 0.577350269189626
3 -0.774596669241483 0 0.774596669241483
4 -0.861136311594053 -0.339981043584856 0.339981043584856 0.861136311594053
5 -0.906179845938664 -0.538469310105683 0 0.538469310105683 0.906179845938664
6 -0.932469514203152 -0.661209386466264 -0.238619186083197 0.238619186083197 0.661209386466264 0.932469514203152
7 -0.949107912342758 -0.741531185599394 -0.405845151377397 0 0.405845151377397 0.741531185599394
8 -0.960289856497536 -0.796666477413627 -0.525532409916329 -0.183434642495650 0.183434642495650 0.525532409916329
9 -0.968160239507626 -0.836031107326636 -0.613371432700590 -0.324253423403809 0 0.324253423403809
10 -0.973906528517172 -0.865063366688985 -0.679409568299024 -0.433395394129247 -0.148874338981631 0.148874338981631
N u7 u8 u9 u10
7 0.949107912342758
8 0.796666477413627 0.960289856497536
9 0.613371432700590 0.836031107326636 0.968160239507626
10 0.433395394129247 0.679409568299024 0.865063366688985 0.973906528517172

参考文献

[1] 森口繁一「数値計算工学」(1989) 岩波書店

Top