9. 数値積分 (2) 変数変換型積分公式, 自動積分

台形則の誤差

台形則では分点数を倍(刻み幅 h を1/2)にすると誤差はおおよそ1/4になった. これは, 被積分関数が解析的な場合台形則の誤差が次のように表されることによる(オイラー・マクローリン展開).

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

 
ここで c_i は h や f に無関係な定数である.

h⁴ 以上の項が十分小さいと仮定すると誤差はほぼ h² に比例すると見なすことができ, h が1/2になると誤差は1/4になる. しかし, 微分不可能な特異点があったりするとこれが成り立たなくなることは「数値積分(1)の数値実験」で見た.

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

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

図からもわかるようにこれは周期関数である. 種々の積分公式で分点数N(横軸)を変えて計算し, その時の相対誤差(縦軸)を対数グラフにすると次のようになった.

意外にも台形則が最も速く収束している. 台形則の結果を表にしてみると, 収束率 = 1/4 と比較すると桁違いに小さく(収束が速く)なっていることがわかる. こうなると図のようにロンバーグ積分はかえって遅くなる.

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

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

 
関数形は次のとおりである. 関数値は x = 8 で 10^-14 のオーダー, x = 10 で 10^-22 のオーダーと急速に0に近づく関数になっている.

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

ここでも台形則が一番速く収束し, N = 15 で最高精度に達した.

減衰が速い関数では端ですぐに微分係数が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公式

φ(t) として次式を使うのがIMT公式である.

  
  φ(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, …) となる. すなわち, 積分区間の両端で微分係数が0になっていて, 台形則の誤差が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]

 
IMT公式の詳細は文献[1], [2]に詳しい.

二重指数型積分公式(DE公式)

対象とする積分を減衰の速い無限積分に変換して台形則の誤差が0に近くなる条件を満たすことにより効率よく積分を計算するのが二重指数型積分公式(DE公式)である. 変換関数 φ(t) は積分区間などによっていくつかの種類が提案されている.

DE公式の詳細は文献[2], [3], [4]に詳しい. 変数変換型積分公式としてはIMT公式がはじめに提案され, その後DE公式が提案された. 現在では主に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₊ で打ち切る. これらの値は許容誤差限界をもとに適切な値を設定する.

端点に特異点がある関数の例で取り上げた I = ∫ √(1-x²)dx [0, 1] = π/4 の関数のグラフを再掲する. この関数は下の数値実験(2)でも使用している.

これを上の変換関数 φ(t) を使って無限積分に変換すると次のようになる.

上側も下側も減衰が速いことがわかる. きざみ幅を0.2としてやると, 上側・下側それぞれ16点で 10^-16 のオーダーまで減衰した. これは台形則で精度よく積分することができる.

(ii) 半無限区間 [0, +∞] の積分

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

(a) 通常の関数の場合

  
  φ(t) = exp(2sinh(t))
  φ'(t) = 2cosh(t)exp(2sinh(t))

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

  
  φ(t) = exp(t - exp(-t))
  φ'(t) = (1 + exp(-t))exp(t - exp(-t))

 
(c) 被積分関数の減衰がさらに急激な場合 (f(x) = f₂(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))

(iv) フーリエ型積分

次の半無限区間の積分(フーリエ型積分)を求める.

  
  Ic = ∫f(x)cos(ωx)dx [a, +∞] または Is = ∫f(x)sin(ωx)dx [a, +∞]

 
次の変換関数を使用して無限区間 [-∞, +∞] の積分に変換する.

  
  Icの場合: x = Mφ(t - π/2M)/ω
  Isの場合: x = Mφ(t)/ω

 
φ(t) は t → -∞ のときに二重指数関数的に φ'(t) → 0, t → +∞ のときに二重指数関数的に φ(t) → t となるような関数である. きざみ幅 h と定数 M は, Mh = π が成り立つように選ばれる. この変換により, xが大きくなるときに公式の分点が cos(ωx) あるいは sin(ωx) のゼロ点に二重指数関数的に近づくようになり, 大きなxについては関数値を計算しなくてもよくなる.

φ(t)はいくつか提案されているが, 例えば次式がある.

  
φ(t) = t/(1 - exp(-ksinh(t))), k = 6 (または 2π) (文献[2], [3], [7])
  または
φ(t) = t/(1 - exp(-2t -α(1 - exp(-t)) - β(exp(t) - 1))), β = 1/4, α = β/sqrt(1 + Mlog(1 + M)/4π) (文献[8])

 
下の数値実験(3)で使用する I = ∫log(x)sin(x)dx [0, +∞] = -γ の関数のグラフを示す.

文献[7]の公式(K = 2π)を使い, Tol = 10^-5 としたときの分点を赤丸で示した. xが大きくなると分点が関数のゼロ点に一致していくのがわかる.

この節の参考文献:
[7] Ooura and Mori, “The double exponential formula for oscillatory functions over the half infinite interval”, J. Comput. Appl. Math., 38 (1991), 353-360.
[8] Ooura and Mori, “A robust double exponential formula for Fourier-type integrals”, J. Comput. Appl. Math., 112 (1999), 229-241.

数値実験 (1)

  
  I = ∫ ex cos x dx [0, 1] = 1.37802461354736

数値実験 (2)

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

数値実験 (3)

  
  I = ∫log(x)sin(x)dx [0, +∞] = -γ

誤差の推定

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

台形則の誤差推定

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

  
  I = ∫ ex cos x dx [0, 1] = 1.37802461354736

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

 
ここでは分点数を倍々にするような順序(分割幅は等間隔で半分)で計算を行った. Iは積分の計算値, I_prevは分点数が半分の時の(つまり1つ前の)計算値である. この場合, 分点数を倍にした時にもう一回最初から全部の分点で計算する必要はなく, 半分はすでに計算済の値がそのまま使える.

誤差推定機能付きのガウス公式

台形則とは異なり, ガウス公式では分点が等間隔でないため, 分点数を倍にするともう一回最初から全部の分点で計算する必要があるのが弱点である.

これを解決するために, N点ガウス公式の計算値を生かしたまま 2N+1点を使った積分公式がガウス・クロンロッド公式である.

N点ガウス公式による積分を G_N で表す(積分区間を [-1, 1] とする).

  
  GN = Σwif(xi) (i = 1～N)

 
ガウス・クロンロッド公式では同じN個の分点 x_i と新たに計算が必要な N+1個の分点 y_i を使って求める積分 K_2N+1 を次式で表す.

  
  K2N+1 = Σaif(xi) + Σbjf(yj)  (i = 1～N, j = 1～N+1)

 
N個の分点 x_i における係数 a_i はガウス公式の係数 w_i とは異なる値になる.

N点ガウス公式 G_N は 2N-1次多項式まで正しい解を与える(すなわち 2N-1次である)が, ガウス・クロンロッド公式 K_2N+1 は 3N+1次(Nが偶数のときは 3N+2次)となるように分点と重みを選ぶ.

x_i は N次のルジャンドル多項式 P_N(x) のゼロ点であるが, 残りの分点 y_i は n+1次の多項式 E_N+1(x) のゼロ点になるとすると, E_N+1(x) は次の直交関係を満たす.

  
  ∫ PN(x)EN+1(x)xkdx = 0 (k = 0, ..., N)

 
上の条件を満たす多項式 E_N+1(x) のゼロ点は連立方程式を解いて求めることができ, 区間(-1, 1)に存在し, P_N(x) のゼロ点と交互に現れることがわかっている. クロンロッドはこの連立方程式を解いて求めたが, この方程式は非常に悪条件で解きにくいため, 後にいくつかの計算法が提案されている.

N = 2～10 のときの分点および重みの値を「積分公式の係数表」に掲載した. ここでは, 下記文献[9]に従い, ヤコビ・クロンロッド行列とよばれる対称三重対角行列を作りその固有値を求める方法により計算した.

ガウス・クロンロッド公式では, 2N+1回の関数評価で 2N+1次の G_N と 3N+1次の K_2N+1 の両方を求めることができ, この2つを比べることにより誤差を推定することができる. 推定誤差は経験的に次式で表される.

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

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

ガウス・クロンロッド公式 K_2N+1 は 3N+1次であるが, 同じ関数評価回数で計算できる G_2N+1 は 4N+1次となる. ガウス・クロンロッド公式は精度的にはガウス公式に劣るが, 誤差の推定ができる利点の方が大きい場合に有効である.

この節の参考文献:
[9] D. P. Laurie, “Calculation of Gauss-Kronrod quadrature rules”, Mathematics of Computation, Vol.66 No.219 (1997).

数値実験 (4)

  
  I = ∫ ex cos x dx [0, 1] = 1.37802461354736

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

自動積分

自動積分ルーチンとは, 積分区間と要求精度を与えると誤差推定値が要求精度を満たすまできざみ幅を順次小さくするなどして積分を繰り返すプログラムである.

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

適応自動積分としては, 誤差推定を行うためと計算済の関数値を再利用するために, 分点が等間隔なシンプソン則やニュートン・コーツ公式が使われてきた. しかし, 現在では精度を追求して, 誤差推定機能を持つガウス・クロンロッド公式のような不等間隔分点ではあるが高精度の公式を使ったものが使われている.

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

適応自動積分ルーチンの例として QUADPACK (文献[6]) の QAG について簡単に説明する. QAG は要求精度を満たすまで自動的に積分区間を部分区間に分割しながらガウス・クロンロッド公式の固定分点のルーチン(QK15, …, QK61)を呼び出すことにより自動積分を行う. XLPackにはQagとして収録されており, 通常の数値積分プログラムとしては最も頻繁に使用される.

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

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

1. まず, 全区間をガウス・クロンロッド公式の固定分点ルーチンにより積分を行い積分値 S と推定誤差 E を求める.
2. E > Tol である限り以下を繰り返す.
   1. 推定誤差が最も大きい部分区間を半分に分割してそれぞれを固定分点ルーチンにより積分し, 積分値 S_i と推定誤差 E_i を求める(_i は部分区間ごとの値を示す).
   2. 新たに構成されたものを含めたすべての部分区間における積分値と推定誤差のそれぞれの総和 S と E を更新する.
3. E <= Tol になったらその時の積分値の総和 S を最終結果として返す.

QAG ではベースの固定分点ルーチンとして QK15 ～ QK61 のうちユーザーが指定したものを使用し, これを自動的に変更することはしない. また, 区間を分割する際に計算済の関数値を再利用することもしない.

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

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

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

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

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

積分公式の係数表

ガウス・クロンロッド公式

  
  K2N+1 = Σ ai f(xi) + Σ bj f(yj)  (第1項のΣは i = 1～N, 第2項のΣは j = 1～N+1)

 
ガウス・クロンロッド公式における分点 x_i, y_j および重み a_i, b_j を以下に示す(x_i はガウス公式の x_i と同じ). 計算は拡張倍精度(80ビット浮動小数)により行い, 結果の16桁を表示した.

N は対応するガウス公式の N の値を示す. ガウス・クロンロッド公式としての分点数は 2N + 1 となる (例. N = 7 は 15分点).

  
N = 2
  xi: -5.773502691896258E-0001  5.773502691896258E-0001
  ai:  4.909090909090909E-0001  4.909090909090909E-0001
  yi: -9.258200997725515E-0001  0.000000000000000E+0000  9.258200997725515E-0001
  bi:  1.979797979797980E-0001  6.222222222222222E-0001  1.979797979797980E-0001
N = 3
  xi: -7.745966692414834E-0001  0.000000000000000E+0000  7.745966692414834E-0001
  ai:  2.684880898683334E-0001  4.509165386584741E-0001  2.684880898683334E-0001
  yi: -9.604912687080203E-0001 -4.342437493468026E-0001  4.342437493468026E-0001
       9.604912687080203E-0001
  bi:  1.046562260264673E-0001  4.013974147759622E-0001  4.013974147759622E-0001
       1.046562260264673E-0001
N = 4
  xi: -8.611363115940526E-0001 -3.399810435848563E-0001  3.399810435848563E-0001
       8.611363115940526E-0001
  ai:  1.700536053357227E-0001  3.269491896014516E-0001  3.269491896014516E-0001
       1.700536053357227E-0001
  yi: -9.765602507375731E-0001 -6.402862174963100E-0001  0.000000000000000E+0000
       6.402862174963100E-0001  9.765602507375731E-0001
  bi:  6.297737366547301E-0002  2.667983404522844E-0001  3.464429818901364E-0001
       2.667983404522844E-0001  6.297737366547301E-0002
N = 5
  xi: -9.061798459386640E-0001 -5.384693101056831E-0001  0.000000000000000E+0000
       5.384693101056831E-0001  9.061798459386640E-0001
  ai:  1.152333166224734E-0001  2.410403392286476E-0001  2.829874178574912E-0001
       2.410403392286476E-0001  1.152333166224734E-0001
  yi: -9.840853600948425E-0001 -7.541667265708492E-0001 -2.796304131617832E-0001
       2.796304131617832E-0001  7.541667265708492E-0001  9.840853600948425E-0001
  bi:  4.258203675108183E-0002  1.868007965564927E-0001  2.728498019125589E-0001
       2.728498019125589E-0001  1.868007965564927E-0001  4.258203675108183E-0002
N = 6
  xi: -9.324695142031520E-0001 -6.612093864662645E-0001 -2.386191860831969E-0001
       2.386191860831969E-0001  6.612093864662645E-0001  9.324695142031520E-0001
  ai:  8.369444044690663E-0002  1.810719943231376E-0001  2.337708641169944E-0001
       2.337708641169944E-0001  1.810719943231376E-0001  8.369444044690663E-0002
  yi: -9.887032026126789E-0001 -8.213733408650279E-0001 -4.631182124753046E-0001
       0.000000000000000E+0000  4.631182124753046E-0001  8.213733408650279E-0001
       9.887032026126789E-0001
  bi:  3.039615411981977E-0002  1.373206046344469E-0001  2.132096522719623E-0001
       2.410725801734648E-0001  2.132096522719623E-0001  1.373206046344469E-0001
       3.039615411981977E-0002
N = 7
  xi: -9.491079123427585E-0001 -7.415311855993944E-0001 -4.058451513773972E-0001
       0.000000000000000E+0000  4.058451513773972E-0001  7.415311855993944E-0001
       9.491079123427585E-0001
  ai:  6.309209262997855E-0002  1.406532597155259E-0001  1.903505780647854E-0001
       2.094821410847278E-0001  1.903505780647854E-0001  1.406532597155259E-0001
       6.309209262997855E-0002
  yi: -9.914553711208126E-0001 -8.648644233597691E-0001 -5.860872354676911E-0001
      -2.077849550078985E-0001  2.077849550078985E-0001  5.860872354676911E-0001
       8.648644233597691E-0001  9.914553711208126E-0001
  bi:  2.293532201052922E-0002  1.047900103222502E-0001  1.690047266392679E-0001
       2.044329400752989E-0001  2.044329400752989E-0001  1.690047266392679E-0001
       1.047900103222502E-0001  2.293532201052923E-0002
N = 8
  xi: -9.602898564975362E-0001 -7.966664774136267E-0001 -5.255324099163290E-0001
      -1.834346424956498E-0001  1.834346424956498E-0001  5.255324099163290E-0001
       7.966664774136267E-0001  9.602898564975362E-0001
  ai:  4.943939500213931E-0002  1.116463708268396E-0001  1.566526061681884E-0001
       1.814000250680346E-0001  1.814000250680346E-0001  1.566526061681884E-0001
       1.116463708268396E-0001  4.943939500213931E-0002
  yi: -9.933798758817162E-0001 -8.941209068474564E-0001 -6.723540709451587E-0001
      -3.607010979281320E-0001  0.000000000000000E+0000  3.607010979281320E-0001
       6.723540709451587E-0001  8.941209068474564E-0001  9.933798758817162E-0001
  bi:  1.782238332071036E-0002  8.248229893135833E-0002  1.362631092551722E-0001
       1.720706085552113E-0001  1.844464057446916E-0001  1.720706085552113E-0001
       1.362631092551722E-0001  8.248229893135833E-0002  1.782238332071036E-0002
N = 9
  xi: -9.681602395076261E-0001 -8.360311073266358E-0001 -6.133714327005904E-0001
      -3.242534234038089E-0001  0.000000000000000E+0000  3.242534234038089E-0001
       6.133714327005904E-0001  8.360311073266358E-0001  9.681602395076261E-0001
  ai:  3.963189516026125E-0002  9.079068168872639E-0002  1.300014068553412E-0001
       1.564135277884839E-0001  1.648960128283494E-0001  1.564135277884839E-0001
       1.300014068553412E-0001  9.079068168872639E-0002  3.963189516026126E-0002
  yi: -9.946781606773402E-0001 -9.149635072496779E-0001 -7.344867651839338E-0001
      -4.754624791124599E-0001 -1.642235636149868E-0001  1.642235636149868E-0001
       4.754624791124599E-0001  7.344867651839338E-0001  9.149635072496779E-0001
       9.946781606773402E-0001
  bi:  1.430477564383894E-0002  6.651815594027414E-0002  1.117891346844183E-0001
       1.452395883843662E-0001  1.628628274401151E-0001  1.628628274401151E-0001
       1.452395883843662E-0001  1.117891346844183E-0001  6.651815594027414E-0002
       1.430477564383894E-0002
N = 10
  xi: -9.739065285171717E-0001 -8.650633666889845E-0001 -6.794095682990244E-0001
      -4.333953941292472E-0001 -1.488743389816312E-0001  1.488743389816312E-0001
       4.333953941292472E-0001  6.794095682990244E-0001  8.650633666889845E-0001
       9.739065285171717E-0001
  ai:  3.255816230796473E-0002  7.503967481091995E-0002  1.093871588022976E-0001
       1.347092173114733E-0001  1.477391049013385E-0001  1.477391049013385E-0001
       1.347092173114733E-0001  1.093871588022976E-0001  7.503967481091995E-0002
       3.255816230796473E-0002
  yi: -9.956571630258081E-0001 -9.301574913557082E-0001 -7.808177265864169E-0001
      -5.627571346686047E-0001 -2.943928627014602E-0001  0.000000000000000E+0000
       2.943928627014602E-0001  5.627571346686047E-0001  7.808177265864169E-0001
       9.301574913557082E-0001  9.956571630258081E-0001
  bi:  1.169463886737187E-0002  5.475589657435200E-0002  9.312545458369761E-0002
       1.234919762620659E-0001  1.427759385770601E-0001  1.494455540029169E-0001
       1.427759385770601E-0001  1.234919762620659E-0001  9.312545458369761E-0002
       5.475589657435200E-0002  1.169463886737187E-0002