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}\)の場合)。
\( N \le 2^{20} \) くらいまでは順調そうだがそれより大きくなると丸め誤差が無視できなくなるようだ。\(I_{2^{20}} = 3.141592653589158068\)であった。
\( N \le 2^9 \) くらいでもう怪しくなっている。\( I_{2^9}^{1} = 3.141592653589792228 \)であった。
ダメダメ。
結果、\(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 |