ロンバーグ積分とリチャドソン補外

# 応用数学研究~数値積分(6章)
### ロンバーグ積分とリチャドソン補外
---
### 誤差評価
#### 今回はロンバーグ積分とニュートンコーツ法を利用して、プログラムを組み誤差を評価する 
---
## ロンバーグ積分とリチャドソン補外

今まで考えてきた、Newton-Cotes法とは違う考え方を使います。
定積分 $ \int_{a}^{b} f\(x\) dx $ を $n = 2^i$ 分割し、Newton-Cotes法の台形則で求めた値を$t^{0}_{i}$と表す。

```math
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 $ とすれれば厳密な値に収束しますが、実際は累積誤差が生じるためそのような結果にはなりません。

---

### ロンバーグ積分のアルゴリズム 

```math
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を用いて、以下のように書ける。（ここが、厳密な証明でない）

```math
e_{h_1} = Ch^{2}_1\\
e_{h_1} = Ch^{2}_1
```
と書ける。

---

言い換えると積分の近似値を$S(h)$で表すと、真値$S$を用いて

```math
S(h_1) = S + Ch^2_1\\
S(h_2) = S + Ch^2_2
```
と表すことができる。
この式からSを求めようとすると、上式を変形すると

```math
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のアルゴリズムを使う。
![](http://mathworld.wolfram.com/images/equations/NevillesAlgorithm/NumberedEquation1.gif)
![](http://mathworld.wolfram.com/images/equations/NevillesAlgorithm/NumberedEquation2.gif)

申し訳ないが、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	|


---

下のグラフは上記の表をグラフ化したものである。
![](https://docs.google.com/spreadsheets/d/e/2PACX-1vSa1_vMM_ZSnWYuThWEwumnRy7Uyra8HrQ5-3zHYaiEO_dQIJOB6fdhxf42JwdSu6dK5MHrqH7BU-0e/pubchart?oid=1903399811&format=image)
![reference](https://docs.google.com/spreadsheets/d/e/2PACX-1vSa1_vMM_ZSnWYuThWEwumnRy7Uyra8HrQ5-3zHYaiEO_dQIJOB6fdhxf42JwdSu6dK5MHrqH7BU-0e/pubchart?oid=1649831740&format=image)

---

#### 次の表はどの分割数で誤差が最小になったのかを考慮する。

|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$の分割数で収束したと考えられる。


<!-- 上述の表より中点則と台形則は明らかに収束していないことがわかるが、シンプソン則とロンバード積分法に関しては、2^8の分割数で一旦収束し、丸め誤差などにより、また発散しようとしている。ことがわかる。 -->

---

最後に、$f(x) = \sqrt{1-x^2}$より$　\Pi $を求める。
実際の誤差のグラフは下記のようになる。
![](https://docs.google.com/spreadsheets/d/e/2PACX-1vSa1_vMM_ZSnWYuThWEwumnRy7Uyra8HrQ5-3zHYaiEO_dQIJOB6fdhxf42JwdSu6dK5MHrqH7BU-0e/pubchart?oid=2068700971&format=image)

![](https://docs.google.com/spreadsheets/d/e/2PACX-1vSa1_vMM_ZSnWYuThWEwumnRy7Uyra8HrQ5-3zHYaiEO_dQIJOB6fdhxf42JwdSu6dK5MHrqH7BU-0e/pubchart?oid=1425048789&format=image)

---
#### 参考文献
<h5>1.[数理解析研究所講究録第889巻1994年26-36　加速法とその漸化式表現　東京大学工学部物理工学科 降旗大介](http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0889-03.pdf)
<h5>2.[わかばテクノロジー-Romberg積分](http://wakaba-technica.sakura.ne.jp/library/integral_romberg.html)
<h5>3.[Richardsonの補外-D. Laurie, “The SIAM 100-Digit Callenge”, SIAM](https://na-inet.jp/nasoft/chap17.pdf)
<h5>4.[プログラマーの為の数学勉強会-中村晃一](http://nineties.github.io/math-seminar/)
<h5>5.[外挿法](http://jun.artcompsci.org/kougi/system_suuri4_1999/note6/node1.html)
<h5>6.[WolframMathWorld-Neville's Algorithm](http://mathworld.wolfram.com/NevillesAlgorithm.html)