数値計算法概論： No.6 数値積分(Richardson 加速)

1 数値積分の収束加速法

数値積分で台形公式を使って、刻み幅を小さくして精度をあげようとすると計算量が増大する。又、無闇に刻み幅を小さくすると今度は丸め誤差の影響がでてくる。そこで刻み幅の荒いデータから収束性を早くする方法を考えよう。

1.1 Euler-Maclaurin 総和公式

f(x) が滑らかな関数の時、次のオイラー・マクローリン ( Euler-Maclaurin ) 総和公式が成立する [1, 2]。

∫_{a}^{b} f(x)dx = h [\frac{1}{2}f(a) + ∑_{i=1}^{i=N−1} f(a + hi) + \frac{1}{2}f(b)]
− \frac{B_{2}h^{2}}{2!} (f′(b) − f′(a)) − ⋯ − \frac{B_{2k}h^{2k}}{(2k)!} (f^{(2k−1)}(b) − f^{(2k−1)}(a)) − ⋯,
h ≡ \frac{(b − a)}{N}  (1)

ここで B_{2k} は Bernoulli 数と呼ばれ、次の関係で定義される。

\frac{t}{e^{t} − 1} = ∑_{n=0}^{∞} B_{n} \frac{t^{n}}{n!}  (2)

具体的には

B_{0} = 1, B_{2} = \frac{1}{6}, B_{4} = −\frac{1}{30}, B_{6} = \frac{1}{42}, ⋯  (3)

(奇数値は B_{1} = −1/2 を除くとすべて消える)。

1.2 収束性

分割数 N の時の台形公式による数値積分の結果を I_{N} とする。

I_{N} = h_{N} [\frac{1}{2}f(a) + ∑_{i=1}^{i=N−1} f(a + h_{N}i) + \frac{1}{2}f(b)]; h_{N} = \frac{b − a}{N}  (4)

台形公式の収束性はオイラー・マクローリンの公式 (1) で与えられ、 h の偶関数である。分割数 N に対して誤差は N^{−2}, N^{−4}, … と振舞う。

I_{N} = I_{∞} + c_{1}N^{−2} + c_{2}N^{−4} + c_{3}N^{−6} ⋯  (5)

ところで I_{N} と I_{∞} の直接比較して収束性評価できるのは結果が解析的に分かっている場合のみである。一般には、 I_{2N} − I_{N} を調べることで収束性を調べられる。例えば誤差が、

I_{N} = I_{∞} + cN^{−α}  (6)

と冪乗的に振舞うなら、 I_{2N} − I_{N} = c′N^{−α}, c′ = c(2^{−α} − 1) となる。さらに

(I_{4N} − I_{2N})/(I_{2N} − I_{N}) = 2^{−α}  (7)

として誤差の収束の早さも評価できる。

1.3 台形公式の精度向上法:Richardson 補外

さらにオイラー・マクローリンの公式 (1) を使うと収束性を改善できる。つまり、

I_{N} = I_{∞} + c_{1}N^{−2} + c_{2}N^{−4} + c_{3}N^{−6} ⋯

と

I_{2N} = I_{∞} + c_{1}(2N)^{−2} + c_{2}(2N)^{−4} + c_{3}(2N)^{−6} ⋯

を比較して、

I_{N}^{(1)} = \frac{I_{2N} − I_{N}/4}{1 − 1/4}  (8)

とすると N^{−2} の誤差の項が打ち消されるので、

I_{N}^{(1)} − I_{∞} = −\frac{c_{2}}{4} N^{−4} − \frac{5c_{3}}{16} N^{−6} ⋯  (9)

となり元の I_{N} より収束が早くなる。この操作を Richardson 加速 [3, 4] と呼ぶ。台形公式による数値積分を Richardson 加速すると、誤差は一般に N^{−4} のように振舞う。

上のようにして得られた I_{N}^{(1)} の数列を、さらに

I_{N}^{(2)} = \frac{I_{2N}^{(1)} − I_{N}^{(1)}/16}{1 − 1/16}  (10)

とすると、誤差が N^{−6} のように振舞う。特に分割数 N を 2, 4, 8, 16, 32, ⋯ と 2 の冪乗 N = 2^{m} にとると、

I_{2^{m}} = I_{∞} + ∑_{j=1}^{∞} c_{j}2^{−2jm}  (11)

であるので、上の操作を繰り返すことができて、

I_{N}^{(k)} = \frac{I_{2N}^{(k−1)} − I_{N}^{(k−1)}/4^{k}}{1 − 1/4^{k}}  (12)

とするとさらに収束を早くできる。

2 プログラム

積分の上限、下限と分割数を与えた時、台形公式で数値積分の結果を出すサブルーチンを利用すると、分割数を 2 の冪乗とした時の結果を出すプログラムは以下のようになる。

c     Main routine
c
      implicit real*8 (a-h,o-z)
      dimension dInt(10,3)
c
      do m=1,6
         n=2**m
         dInt(m,1)=dIntegral(n,a,b)
      end do
c
      do m=1,6-1
         write(*,*) 2**m, dInt(m,1),dInt(m+1,1)-dInt(m,1)
      end do
      end
c
c
c     Integral routine(trapezoidal integral)
c
      real*8 function dIntegral(n,a,b)
      implicit real*8 (a-h,o-z)
      integer n,i
c
      ...........
      dIntegral=......
      return
      end

また、これを利用すると、 Richardson 補外のプログラムができる。複数回 Richardson 補外を繰り返すには、 subroutine 文を使うと有効である。

3 課題

(1-1) No.5 の積分に付いて、分割数 N を変えて I_{N} を計算し O(h^{2})(= O(1/N^{2})) の誤差の収束性を確かめよ。この中で、収束性が予期されるものよりずっと遅いものがある。その原因を考察して見よ。

(1-2) さらに Richardson 補外をして I_{N}^{(1)}, I_{N}^{(2)} … の収束性を見よ。以上のものの内、収束性が予想より早いものと、ずっと遅いものとがある。その原因を考察して見よ。

参考文献

[1] 「Numerical Recipes in C」、技術評論社、第 4 章 2 節

[2] 「数値計算法の数理」、杉原正顕、室田一雄著、岩波書店、第 12 章 3 節 [3] 「数値計算の常識」、伊理正夫、藤野和建著、共立出版、第 7 章

[4] 「数値計算法の数理」、杉原正顕、室田一雄著、岩波書店、第 12 章 4 節