9. 数値積分 (2) 変数変換による積分公式、誤差の推定、自動積分

種々の基本公式を見てきたが、実際に数値積分を行う際には基本公式をそのまま使うことはなく、さまざまな工夫を加えたものが実用ルーチンとして使用される。以下、いくつかの話題について掘り下げてみる。

(1) 台形則の誤差

台形則では、分点数を倍(刻み幅 h を1/2)にすると誤差はおおよそ1/4になった。これは、被積分関数が解析的な場合、台形則の誤差が次のように表されることによる (Euler-Maclaurinの公式)。

    Im - I = c2h2[f'(b) - f'(a)] + c4h4[f'''(b) - f'''(a)] + ... + c2ph2p[f(2p-1)(b) - f(2p-1)(a)] + O(h2p+2)

ここで、ci は h や f に無関係な定数である。

h4以上の項が十分小さいとすると誤差はほぼh2に比例すると見なすことができ、h が1/2になると誤差は1/4になる。しかし、微分不可能な特異点があったりするとこれが成り立たなくなることは基本公式の例題で見た。

次に、特別な場合の例として、次の関数の[0, π/2]での積分を求めてみる。

    I = ∫1/sqrt(1 - sin2x/2) dx [0, π/2] = 1.85407 46773 01372

ex5f

図からもわかるように、これは周期関数である。種々の公式で計算し、誤差をプロットすると次のようになった。

ex5

意外にも台形則が最も速く収束した。下表のように桁違い(1/4と比較すると)に収束が速くなっている。

分点数 積分値(台形則の計算値) 誤差 収束率
2 1.85495913108563 8.8445378e-4
4 1.85407522776731 5.5046594e-7 1/1.607e+3
8 1.85407467730167 2.9465319e-13 1/1.868e+6
16 1.85407467730137 0.0

これは、周期性のある関数はその両端で微分係数が一致し、誤差を表す式(Euler-Maclaurinの公式)が0に近くなるためである。上の例では積分範囲が半周期であるが対称性があるためこの条件にあてはまる。このような場合には台形則が有効であることが知られている。

もうひとつ、無限積分の例を見てみる。

    ∫exp(-x^2)/2 dx [0, ∞] = 1.25331 41373 15500

関数形は次のとおりである。関数値は、x=8で-14乗のオーダー、x=10で-22乗のオーダーになっている。

ex6f

ここでは積分範囲を[0, ∞]の代わりに[0, 8]として計算してみると、次のようになった。

ex6

ここでも台形則が一番速く収束し、N=12で最高精度に達した。

減衰が速い関数では端ですぐに微分係数が0に近づくため、上の誤差を表す式が速く0に近づくと考えられる。このような無限積分に対しては、台形則が他の公式よりも精度が良いことが知られている。

(2) 変数変換型積分公式

台形則は周期性のある関数の積分や減衰の速い無限積分に対して有効であることがわかったが、実際には台形則がそのまま使われることはない。しかし、台形則のこの特長は以下に示す変数変換による積分公式に応用されている。

一般の積分を台形則が有効な積分になるように変数変換し、それに台形則を適用すれば性能の良い積分公式となることが期待できる。このアイデアに基づき、変数変換型積分公式が開発されている。

積分変数を変換する公式は次のとおりである。

    ∫ f(x)dx [a, b] = ∫ g(t)dt [c, d]
    ただし、x=φ(t), g(t) = f(φ(t))φ'(t)

ここで、φ(t)は区間[c, d]で単調で、かつ a=φ(c), b=φ(d) である。これにより、積分値が不変のまま被積分関数を変えることができる。適切なφ(t)を見つけることができれば効果的な公式を作ることができる。

IMT公式 (文献[1], [2])

φ(t) として次式を使う。

    φ(t) = (1/Q)∫ψ(r)dr [0, t], Q = ∫ψ(r)dr [0, 1]
    ただし、ψ(r) = exp(-1/r - 1/(1-r))

この変換を区間[0, 1]の積分に適用すると、変換後の関数g(t)は g(j)(0) = g(j)(1) = 0 (j=0, 1, …) を満たす。従って、台形則の誤差(Euler-Maclaurinの公式)は0に近くなり、効率よく積分を計算することができる。

分点数をNとして、計算式は次のようになる。

    SN = (1/N)Σwi(N)f(xi(N)) (i=1~N-1)
    重み wi(N) = φ'(i/N) = (1/Q)ψ(i/N)
    分点 xi(N) = φ(i/N) = (1/Q)∫ψ(r)dr [0, i/N]

任意の区間[a, b]の積分については、x = (b-a)t + a により区間[0, 1]の積分に変換してからこの公式を適用すればよい。

二重指数型積分公式 (DE公式) (文献[2], [3])

DE公式では、対象とする積分を減衰の速い無限積分に変換してから台形則を適用することにより効率よく積分を計算する。

(i) 有限区間[-1, 1]の積分

φ(t) として次を使う。

    φ(t) = tanh((π/2)sinh(t))
    φ'(t) = (π/2)cosh(t) / cosh2((π/2)sinh(t))

きざみをh, 分点数をN (= N++N+1)とすると、計算式は次のようになる。

    Sh(N) = hΣwi(N)f(xi(N)) (i=-N-~N+)
    重み wi(N) = hφ'(ih)
    分点 xi(N) = φ(ih)

無限積分の計算は、下側を-N、上側をN+で打ち切ることとした。これらの値は許容誤差限界をもとに適切な値を設定する。

任意の区間[a, b]の積分については、x = ((b-a)/2)t + (b+a)/2 により区間[-1, 1]の積分に変換してからこの公式を適用する。

(ii) 半無限区間[0,+∞]の場合

半無限区間の積分において被積分関数がゆるやかに減衰するような場合についても、有限区間の場合と同様の考え方でDE公式を導くことができる。文献[3]に示されているプログラムDEHINTでは、被積分関数により3種類の公式を使い分けるようになっている。

(a) 通常の関数の場合
    φ(t) = exp(2sinh(t))
    φ'(t) = 2cosh(t)exp(2sinh(t))

(b) 被積分関数が指数関数的に減衰する場合 (f(x) = f1(x)exp(-x) の形をしている場合)
    φ(t) = exp(t - exp(-t))
    φ'(t) = (1 + exp(-t))exp(t - exp(-t))

(c) 被積分関数の減衰がさらに急激な場合 (f(x) = f2(x)exp(-x^2) の形をしている場合)
    φ(t) = exp(t/2 - exp(-t))
    φ'(t) = (1/2 + exp(-t))exp(t/2 - exp(-t))

(iii) 全無限区間[-∞,+∞]の場合

被積分関数がゆるやかに減衰するような全無限区間の積分についても同様にDE公式を導くことができる。

    φ(t) = sinh((π/2)sinh(t))
    φ'(t) = (π/2)cosh(t)cosh((π/2)sinh(t))

数値実験

まずは、前に実験した例を用いて、変数変換型公式の性能を見てみる。
次の積分を計算してみる。

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

DE公式ではプログラムの作りによって関数評価回数が変わってくるので、DE公式だけ横軸をXLPackのDefinルーチンを使用したときの実際の関数評価回数とした。

ex10

IMT公式、DE公式ともに特別に高性能ではなく、Gauss公式には及ばないようである。なお、変数変換型公式としてはIMT公式が最初に提案され、DE公式は後で開発された。DE公式の性能が良いので、現在では変数変換型公式としては主にDE公式が用いられている。

次に、同様に前に実験した例を用いて、端点に特異点がある場合について見てみる。

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

x=1が特異点である。

ex11

前にも見たように補間型公式ではどれもまともな結果が得られないが、IMT公式とDE公式は特異点の影響を受けていないのがわかる。このように、IMT公式とDE公式の特徴は、端点に特異点がある(端点で微係数が定義できない)場合にも適用できることである。なお、積分区間内に特異点がある場合には良い結果が得られない。

(3) 誤差の推定

「数値積分(1) 基本公式」で見たように、関数評価回数が同程度であれば種々の公式の中でGauss公式が最も精度が良かった。単に計算するのであればこれで良いが、実務的には積分値と同時に推定誤差も知りたい。また、後で示す自動積分のベースとして使う場合には誤差推定ができなければならない。Gauss公式の弱点は誤差の推定がしにくいことである。

台形則の誤差推定

計算結果の誤差を推定するためには、例えば、N分点で計算したものと2N分点で計算したものを比較するなどすればよい。次の例(滑らかな関数)で、台形則の計算例を見てみる。

    I = ∫01 ex cos x dx = 1.37802 46135 47364

分点数 計算値(台形則) 真値 誤差 誤差推定値
2 1.34061800327106 1.37802461354736 0.037406610 0.035423678
4 1.36858238253106 1.37802461354736 0.009442231 0.009321460
8 1.37565843490021 1.37802461354736 0.002366179 0.002358684
16 1.37743271822098 1.37802461354736 0.000591895 0.000591428
32 1.37787661780930 1.37802461354736 0.000147996 0.000147967

誤差推定値はRomberg積分の時と同じように次の式から計算した。実際の誤差と推定値はよく一致している。

    推定誤差 = (4I - Iprev)/3

ここでは、分点数を倍々にするような順序(分割幅は等間隔で半分)で計算を行った。Iは積分の計算値、Iprevは分点数が半分の時の(つまり1つ前の)計算値である。この場合、分点数を倍にした時にもう一回最初から全部の分点で計算する必要はなく、半分はすでに計算済の値がそのまま使える。ところが、Gauss公式では分点が等間隔でないためそのようなことができないのが弱点となっている。

誤差推定機能付きのGauss公式

これを解決するために、N点Gauss公式の計算値を生かしたまま2N+1点を使った積分公式がGauss-Kronrod公式である。

N点Gauss公式による積分をGNで表す(積分区間を[-1, 1]とした場合)。

    GN = Σ wi f(xi) (i=1~N)

この時、Kronrod公式は同じN個の点xiと新たに計算が必要なN+1個の点yiを使って、

    K2N+1 = Σ ai f(xi) + Σ bj f(yj)  (i=1~N, j=1~N+1)

と表される。なお、同じN個の点における係数aiはGauss公式の係数wiとは異なる値を持つ。

N点Gauss公式GNは2N-1次多項式まで正しい解を与える(2N-1次である)が、Kronrod公式K2N+1は3N+1次となる。同じ関数評価回数で計算できるG2N+1は4N+1次となるので、Kronrod公式は精度的にはGauss公式に劣るが、誤差の推定ができる利点の方が大きい場合に有効である。

Gauss-Kronrod公式では、2N+1回の関数評価で2N+1次のGNと3N+1次のK2N+1の両方を求めることができるので、この2つを比べることにより誤差を推定することができる。推定誤差は経験的に次式で表される。

    推定誤差 = (200*|GN - K2N+1|)1.5

この式による推定誤差は安全側に求められるとされている。

Gauss-Kronrod公式の重みおよび分点の具体的な値を下の付表に示す(Gauss公式への追加分のみ)。

数値実験

Gauss-Kronrod公式を用いて次の積分を求めてみる。

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

これは前にも例題に使用した性質のよい関数の例であり、計算結果は次のようになった。

ex20

横軸に関数評価回数、縦軸に誤差をプロットした。緑色のデータは前項の式で求めた推定誤差である。各データ点に表示した数字はNの値である。

こうしてみるとKronrod公式でも思ったよりも精度が落ちていないように見える。なお、推定誤差はかなり安全側に算出されているのがわかる。

この例題では収束が速すぎてうまく比較できないので、もう少し難しい次の積分を計算してみる。

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

これも前に例題に使用した関数で、積分区間外であるが近くに特異点を持つ例である。

ex21

計算回数がかなり増えたが傾向は同じであった。Kronrod公式を使うと、それほど大きなペナルティを払わずに誤差を推定できることが確かめられた。

(4) 自動積分

自動積分ルーチンは、積分区間と要求精度を与えると、誤差推定値が要求精度を満たすまで刻み幅を順次小さくして積分を繰り返すものである。

適応自動積分ルーチンは、これを一歩進めて、積分区間を部分区間に分割し、誤差推定値が要求された精度を満たすように部分区間の幅を自動的に調節する。被積分関数がなめらかにゆっくり変化するような区間では大きな幅を用い、誤差推定値が大きい部分では細かい区間に分割して刻み幅を小さくして計算することにより、最小の計算量で要求精度の結果を得ようとするものである。

適応自動積分では、誤差推定を行うためと計算済の関数値を再利用するために、分点が等間隔なSimpson公式やNewton-Cotes公式が使われてきた (例、QUANC8 (文献[4]))。しかし、よく使われているQAG (文献[5])のように誤差推定機能を持つGauss-Kronrod公式のような不等間隔分点の公式を使ったものもある。

適応自動積分プログラム QAG

適応自動積分ルーチンの例としてQAGについて簡単に説明する。QAGは、要求精度を満たすまで自動的に積分区間を部分区間に分割しながらGauss-Kronrod公式の固定分点のルーチンを呼び出すことにより自動積分を行う。XLPackには倍精度版(Dqag)が収録されている。

このプログラムの動作の概要は次のとおりである。

与えられた収束条件を「推定誤差 <= Tol」とする。

  1. まず、全区間をGauss-Kronrodルーチンにより積分を行い、積分値Sと推定誤差Eを求める。
  2. E > Tol である限り、以下を繰り返す
    1. 推定誤差が最も大きい部分区間を半分に分割してそれぞれをGauss-Kronrodルーチンにより積分する
    2. 積分値Sと推定誤差Eを、新たに構成されたものを含めたすべての部分区間における積分値と推定誤差それぞれの総和に更新する
  3. E <= Tol になったら、その時の積分値の総和Sを最終結果として返す

QAGではベースのGauss-KronrodルーチンとしてQK15(N=7)~QK61(N=30)のうちユーザーが指定したものを使用し、これを自動的に変更することはしない。また、区間を分割する際に計算済の関数値を再利用することもしない。シンプルなアルゴリズムではあるがよく機能する。

例として、区間[0,1]において次の関数を積分した時のQAGの挙動を下図に示す (文献[4]の例)。

    f(x) = 1/((x - 0.3)2 + 0.01) + 1/((x - 0.9)2 + 0.04) - 6

ex22

ここでは、XLPackに収録されているDqagを使用し、目的誤差を10-5以下としてベースとしてQK15を指定して計算した。区間は5つに分割された (0~0.25, 0.25~0.375, 0.375~0.5, 0.5~0.75, 0.75~1)。トータルの関数評価回数は135回であった。最も急峻な変化のところで分割されているのがわかる。

試しに同じ条件でQK61を指定してみると次のようになった。

ex23

2分割で計算しており (0~0.5, 0.5~1)、トータルの関数評価回数は183回になった。変化が激しい場合にはQK15を選択し分割数を多くし、そうでない場合にはQK61のように高精度のものを選び分割数を少なくするとよさそうである。

付表. 積分公式の係数

Gauss-Kronrod公式

    K2N+1 = Σ ai f(xi) + Σ bj f(yj)  (i=1~N, j=1~N+1)

において、係数 ai, bjおよび yj を下表に示す(xiはGauss公式のuiと同じ)。

N a1 a2 a3 a4 a5 a6
2 0.490909090909091 0.490909090909091
3 0.268488089868333 0.450916538658474 0.268488089868333
4 0.170053605335723 0.326949189601452 0.326949189601452 0.170053605335723
5 0.115233316622474 0.241040339228648 0.282987417857491 0.241040339228648 0.115233316622474
6 0.0836944404469063 0.181071994323137 0.233770864116994 0.233770864116994 0.181071994323137 0.0836944404469063
7 0.0630920926299790 0.140653259715526 0.190350578064785 0.209482141084728 0.190350578064785 0.140653259715526
8 0.0494393950021390 0.111646370826839 0.156652606168188 0.181400025068035 0.181400025068035 0.156652606168188
9 0.0396318951602617 0.0907906816887259 0.130001406855341 0.156413527788484 0.164896012828349 0.156413527788484
10 0.0325581623079661 0.0750396748109198 0.109387158802298 0.134709217311473 0.147739104901338 0.147739104901338
N a7 a8 a9 a10
7 0.0630920926299790
8 0.111646370826839 0.0494393950021390
9 0.130001406855341 0.0907906816887259 0.0396318951602617
10 0.134709217311473 0.109387158802298 0.0750396748109198 0.0325581623079661
N b1 b2 b3 b4 b5 b6
2 0.197979797979798 0.622222222222222 0.197979797979798
3 0.104656226026467 0.401397414775962 0.401397414775962 0.104656226026467
4 0.0629773736654729 0.266798340452285 0.346442981890136 0.266798340452285 0.0629773736654729
5 0.0425820367510818 0.186800796556493 0.272849801912559 0.272849801912559 0.186800796556493 0.0425820367510818
6 0.0303961541198197 0.137320604634447 0.213209652271962 0.241072580173465 0.213209652271962 0.137320604634447
7 0.0229353220105292 0.104790010322250 0.169004726639268 0.204432940075299 0.204432940075299 0.169004726639268
8 0.0178223833207103 0.0824822989313584 0.136263109255172 0.172070608555211 0.184446405744692 0.172070608555211
9 0.0143047756438389 0.0665181559402741 0.111789134684418 0.145239588384366 0.162862827440115 0.162862827440115
10 0.0116946388673719 0.0547558965743520 0.0931254545836976 0.123491976262066 0.14277593857706 0.149445554002917
N b7 b8 b9 b10 b11
6 0.0303961541198197
7 0.104790010322250 0.0229353220105292
8 0.136263109255172 0.0824822989313584 0.0178223833207103
9 0.145239588384366 0.111789134684418 0.0665181559402741 0.0143047756438389
10 0.14277593857706 0.123491976262066 0.0931254545836976 0.0547558965743520 0.0116946388673719
N y1 y2 y3 y4 y5 y6
2 -0.925820099772551 0 0.925820099772551
3 -0.96049126870802 -0.434243749346802 0.434243749346802 0.96049126870802
4 -0.976560250737573 -0.64028621749631 0 0.64028621749631 0.976560250737573
5 -0.984085360094842 -0.754166726570849 -0.279630413161783 0.279630413161783 0.754166726570849 0.984085360094842
6 -0.988703202612679 -0.821373340865028 -0.463118212475305 0 0.463118212475305 0.821373340865028
7 -0.991455371120813 -0.864864423359769 -0.586087235467691 -0.207784955007898 0.207784955007898 0.586087235467691
8 -0.993379875881716 -0.894120906847456 -0.672354070945159 -0.360701097928132 0 0.360701097928132
9 -0.994678160677340 -0.914963507249678 -0.734486765183934 -0.47546247911246 -0.164223563614987 0.164223563614987
10 -0.995657163025808 -0.930157491355708 -0.780817726586417 -0.562757134668605 -0.29439286270146 0
N y7 y8 y9 y10 y11
6 0.988703202612679
7 0.864864423359769 0.991455371120813
8 0.672354070945159 0.894120906847456 0.993379875881716
9 0.47546247911246 0.734486765183934 0.914963507249678 0.994678160677340
10 0.29439286270146 0.562757134668605 0.780817726586417 0.930157491355708 0.995657163025808

参考文献

[1] 森口繁一「数値計算工学」岩波書店 (1989)
[2] 杉原正顕、室田一雄「数値計算法の数理」岩波書店 (1994)
[3] 森正武「FORTRAN77数値計算プログラミング(増補版)」岩波書店 (1987)
[4] 森正武(訳)「計算機のための数値計算法」科学技術出版社 (1978)
[5] R. Piessens, E. de Doncker, C. Uberhuber, D. Kahaner, “QUADPACK: A subroutine package for automatic integration”, Springer-Verlag (1983)

Top