Richardson補外とか丸め誤差について

数値計算の常識』(伊理正夫・藤野和建著)を読んでいろいろ試したくなったので、取り敢えず簡単にできるものをやってみる。

積分\( I = \int_{a}^{b} f(x) dx\)をN等分した台形\(I_N\)で近似するとき、
\( I_{1} = \frac{1}{2}(b - a)(f(a) + f(b)) \)
\( I_{2^{N+1}} = \frac{1}{2} (I_{2^N} + h \sum_{i = 1}^{2^N} f(a + (i - \frac{1}{2})h) ) \), 但し\( h = \frac{b - a}{2^N}\)
で近似値を(逐次的に)出すことができる。

また、f(x)が\(C^{2p}\)級のとき、hに関係しない係数\( d_{2j}\)を用いて
(*)\( I_N - I = \sum_{i=1}^{p} d_{2i} h^{2i} + O(h^{2p+2})\)
と表すことができるらしい(間違っていたらすいません)。

(*)より、上手く数値積分ができているならば、k=1,...,pとして
\( I_{N}^{k} = \frac{ I_{N}^{k - 1} - 4^{-k} I_{N / 2}^{k - 1}}{ 1 - 4^{-k}} \)
但し、\( I_{N}^{0} = I_{N}\)と定めたとき[Richardson補外]、少し計算すると\( I_{N}^{k} - N\)はNが2倍になると\(1 / 4^{k+1}\)倍になること、即ち\( I_{N}^{k} - I_{N/2}^{k}\)はNが2倍になると\(1 / 4^{k+1}\)倍になることが分かる。これより、もしも\( I_{N}^{k}\)の階差を対数目盛でみたときに直線にのっていなければ丸め誤差等の影響で健全な計算が行えていないことが分かる。ちなみに、k=1の場合に得られる近似式はSimpsonの式に一致する。


例として次の(よく例題に出てくる)積分数値計算で求めることを考える。
\( I = \int_{-1}^{1} \frac{2}{1 + x^2} dx = \pi \)

2^34(=17179869184)等分まで倍精度で求めた。更に、ここから\( I_{N}^{1} \)と \(I_{N}^{2}\)を計算した。

\(I_{N}^{k}\)階差をプロットしたものが下図(順に\(I_{N}, I_{N}^{1}, I_{N}^{2}\)の場合)。

f:id:lan496:20150906011606p:plain
\( N \le 2^{20} \) くらいまでは順調そうだがそれより大きくなると丸め誤差が無視できなくなるようだ。\(I_{2^{20}} = 3.141592653589158068\)であった。

f:id:lan496:20150906011615p:plain
\( N \le 2^9 \) くらいでもう怪しくなっている。\( I_{2^9}^{1} = 3.141592653589792228 \)であった。

f:id:lan496:20150906011620p:plain
ダメダメ。


結果、\(I_{N}^{1}\)を活用することで小数点桁数13桁くらい(3.1415926535897...)まで求めることができた。

実際には直線に乗らない範囲で円周率により近い値を得ているが、本来真値を知っている訳ではないのでそれをより真値に近いものと判断する根拠は乏しい。

[参考]

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <cstdio>
using namespace std;

double a = -1.0;
double b = 1.0;

double func(double x){
  return 2.0 / (1.0 + x * x);
}

int main() {
  int M = 34;
  vector<double> integral(M);
  integral[0] = (b - a) * (func(a) + func(b)) / 2.0;
  
  for(int i=1;i<M;++i){
    double middle = 0.0;
    long long N = pow(2, i - 1);
    double h = (b - a) / N;
    
    for(long long i=1;i<=N;i++){
      middle += h * func(a + (i - 0.5) * h);
    }
    
    integral[i] = (integral[i-1] + middle) / 2.0;
  }
  for(int i=0;i<M;i++){
    printf("%lld %.18lf\n", (long long)pow(2, i), integral[i]);
  }
  return 0;
}
N \(I_{N}\) \(I_{N}^{1}\) \( I_{N}^{2}\)
1 2.000000000000000000 0.000000000000000000 0.000000000000000000
2 3.000000000000000000 3.333333333333333037 0.000000000000000000
4 3.100000000000000089 3.133333333333333304 3.555555555555555358
8 3.131176470588235450 3.141568627450980422 3.119999999999999662
16 3.138988494491089298 3.141592502458707248 3.142117647058823682
32 3.140941612041388886 3.141592651224821786 3.141594094125888859
64 3.141429893174974453 3.141592653552836456 3.141592661142562637
128 3.141551963485655019 3.141592653589214912 3.141592653708037641
256 3.141582481063752041 3.141592653589784234 3.141592653591640083
512 3.141590110458282403 3.141592653589792228 3.141592653589822426
1024 3.141592017806918768 3.141592653589797557 3.141592653589792672
2048 3.141592494644073863 3.141592653589791784 3.141592653589798001
4096 3.141592613853361193 3.141592653589790007 3.141592653589791340
8192 3.141592643655684025 3.141592653589791784 3.141592653589790007
16384 3.141592651106267731 3.141592653589795336 3.141592653589791784
32768 3.141592652968914656 3.141592653589796669 3.141592653589795781
65536 3.141592653434578608 3.141592653589800221 3.141592653589796669
131072 3.141592653550970837 3.141592653589768247 3.141592653589800221
262144 3.141592653580073335 3.141592653589774020 3.141592653589766471
524288 3.141592653587365724 3.141592653589796669 3.141592653589774464
1048576 3.141592653589158068 3.141592653589754924 3.141592653589798001
2097152 3.141592653589570183 3.141592653589707851 3.141592653589752260
4194304 3.141592653589582618 3.141592653589587059 3.141592653589704742
8388608 3.141592653589698081 3.141592653589736273 3.141592653589578621
16777216 3.141592653589746931 3.141592653589763362 3.141592653589746487
33554432 3.141592653589535988 3.141592653589465378 3.141592653589765582
67108864 3.141592653589374340 3.141592653589320605 3.141592653589445838
134217728 3.141592653588935136 3.141592653588788586 3.141592653589311279
268435456 3.141592653589452944 3.141592653589625250 3.141592653588753059
536870912 3.141592653589918793 3.141592653590073780 3.141592653589681206
1073741824 3.141592653590653761 3.141592653590898454 3.141592653590103534
2147483648 3.141592653591862572 3.141592653592265805 3.141592653590953521
4294967296 3.141592653595981055 3.141592653597353735 3.141592653592356843
8589934592 3.141592653606938512 3.141592653610590702 3.141592653597693019