応用数学研究~数値積分(6章)
ロンバーグ積分とリチャドソン補外
誤差評価
今回はロンバーグ積分とニュートンコーツ法を利用して、プログラムを組み誤差を評価する
ロンバーグ積分とリチャドソン補外
今まで考えてきた、Newton-Cotes法とは違う考え方を使います。
定積分 $ \int_{a}^{b} f(x) dx $ を $n = 2^i$ 分割し、Newton-Cotes法の台形則で求めた値を$t^{0}_{i}$と表す。
T^{0}_{i} = h_{i} \sum^{2^i-1}_{m=1}f(x_m)+h_{i}\frac{f(a)+f(b)}{2}\\\
h_{i} = \frac{b-a}{2^{i}}
すると、 $ T^0_0,T^0_1, \cdots $ となり、理論的には $ i \rightarrow \infty $ とすれれば厳密な値に収束しますが、実際は累積誤差が生じるためそのような結果にはなりません。
ロンバーグ積分のアルゴリズム
T^{j}_{i}=\frac{4^jT^{j-1}_i-T^{j-1}_{i-1}}{4^j-1}\\
n = 2^{i}
によって外装を行い$T^{n}_{n}$を求める
リチャドソン補外を用いて補外をしている。
補間法と補外法に関して
補間は二点間の間の点を埋めるようなものであるが、補外は二点間の領域外の点を埋めることである。
リチャードソン補外の概論
非常に申し訳ないが、いろいろな本でロンバーグ積分に関して調べたが、ほとんど資料が見つからず、概説しか話すことができない。
具体的には、台形公式に関して考えてみる。刻み幅をhとし、誤差は非常に大雑把に言って$h^2$になるはずである。（教科書の台形則の誤差参照）
つまり、刻み幅が$h_1,h_2$のに種類を取ったとすると、誤差は共通な定数Cを用いて、以下のように書ける。（ここが、厳密な証明でない）
e_{h_1} = Ch^{2}_1\\
e_{h_1} = Ch^{2}_1
と書ける。
言い換えると積分の近似値を$S(h)$で表すと、真値$S$を用いて
S(h_1) = S + Ch^2_1\\
S(h_2) = S + Ch^2_2
と表すことができる。
この式からSを求めようとすると、上式を変形すると
h^2_2S(h_1) = Sh^2_2 + Ch^2_1h^2_2\\
h^2_1S(h_2) = Sh^2_1 + Ch^2_1h^2_2\\
(上の式から下式を引いて、一つにまとめる)\\
h^2_2S(h_1) - h^2_1S(h_2) = Sh^2_2 - Sh^2_1\\
S = \frac{h^2_2 S(h_1) -h^2_1 S(h_2)}{h^2_2 -h^2_1}
とすると、真値が定まる。
この議論は誤差が$h^2$に対し1次式で表せると考え、$h \rightarrow 0$での積分の値を求めるわけである。
まあ、誤差を$Ch^2$と仮定しているのが厳密には正しくないので、真値は求まらない。
ここで、Nevilleのアルゴリズムを使う。
申し訳ないが、Nevilleのアルゴリズムついて調べようとしたら、詳しく乗っている文献が見つからず今回は、補外法に関しての説明はここまでにしていただく。
プログラムと誤差評価
以下の表は真値から近似値を引いた絶対誤差を表している。
$$
\int _{0}^{1} \frac{1}{1+x^2}dx = \frac{\pi}{4}
$$
となっており、今回は$\frac{\pi}{4}$との絶対誤差を表している。
4種類の方法と真値からの絶対誤差に関して、以下のように表にした。
| 分割数 | Retangle | Trapezoidal | Simpuson | Romberg |
|---|---|---|---|---|
| $2^0$ | +2.85e-01 | +3.54e-02 | +2.06e-03 | +7.85e-01 |
| $2^1$ | +1.02e-01 | +1.04e-02 | +6.01e-06 | +3.54e-02 |
| $2^2$ | +3.65e-02 | +2.60e-03 | +3.78e-08 | +2.06e-03 |
| $2^3$ | +1.29e-02 | +6.51e-04 | +5.91e-10 | +1.31e-04 |
| $2^4$ | +4.58e-03 | +1.63e-04 | +9.24e-12 | +1.72e-06 |
| $2^5$ | +1.62e-03 | +4.07e-05 | +1.44e-13 | +2.92e-09 |
| $2^6$ | +5.74e-04 | +1.02e-05 | +2.33e-15 | +1.21e-11 |
| $2^7$ | +2.03e-04 | +2.54e-06 | +4.44e-16 | +1.75e-14 |
| $2^8$ | +7.18e-05 | +6.36e-07 | +2.22e-16 | +1.11e-16 |
| $2^9$ | +2.54e-05 | +1.59e-07 | +3.33e-16 | +7.77e-16 |
| $2^{10}$ | +8.97e-06 | +3.97e-08 | +8.88e-16 | +1.11e-16 |
| $2^{11}$ | +3.17e-06 | +9.93e-09 | +1.78e-15 | +1.44e-15 |
| $2^{12}$ | +1.12e-06 | +2.48e-09 | +2.11e-15 | +2.55e-15 |
| $2^{13}$ | +3.97e-07 | +6.21e-10 | +1.11e-16 | +2.78e-15 |
下のグラフは上記の表をグラフ化したものである。
次の表はどの分割数で誤差が最小になったのかを考慮する。
| Retangle | Trapezoidal | Simpuson | Romberg |
|---|---|---|---|
| $2^{13}$ | $2^{13}$ | $2^{8}$ | $2^{8}$ |
誤差評価でe−15~e−16はIEEE754の64bitの浮動小数点では仮数部が52bitのため桁数に直すとっだいたい、15から16桁くらいである。そのため丸め誤差が生じるため、simppusonとRombergは、たかだか$2^8$の分割数で収束したと考えられる。
最後に、$f(x) = \sqrt{1-x^2}$より$ \Pi $を求める。
実際の誤差のグラフは下記のようになる。