ロンバーグ積分 [ 1] [ 2] [ 3] [ 4] ロンベルク積分 は、関数 の数値積分 アルゴリズムのひとつである。この方法では台形公式 とリチャードソンの補外 を組み合わせ離散化幅
h
{\displaystyle h}
1955年にヴェルナー・ロンベルク （英語版 ） [ 5]
ロンバーグ積分は区間
[
a
,
b
]
{\displaystyle [a,b]}
f
{\displaystyle f}
∫
a
b
f
(
x
)
d
x
{\displaystyle \int _{a}^{b}f(x)dx}
を数値的に求めるアルゴリズムである。積分区間
[
a
,
b
]
{\displaystyle [a,b]}
h
=
(
b
−
a
)
/
N
{\displaystyle h=(b-a)/N}
N
{\displaystyle N}
[
x
i
,
x
i
+
1
]
{\displaystyle [x_{i},x_{i+1}]}
x
i
=
a
+
i
h
{\displaystyle x_{i}=a+ih}
i
=
0
,
1
,
…
,
N
−
1
{\displaystyle i=0,1,\dots ,N-1}
台形公式 は求める定積分を
T
(
h
)
=
h
[
1
2
f
(
x
0
)
+
f
(
x
1
)
+
⋯
+
f
(
x
N
−
1
)
+
1
2
f
(
x
N
)
]
{\displaystyle T(h)=h\left[{\frac {1}{2}}f(x_{0})+f(x_{1})+\dots +f(x_{N-1})+{\frac {1}{2}}f(x_{N})\right]}
により近似する[ 6] [ 7]
T
(
h
)
{\displaystyle T(h)}
O
(
h
2
)
{\displaystyle {\mathcal {O}}(h^{2})}
オイラーの和公式 ）によりそれをより詳しく評価することができる。すなわち、
f
{\displaystyle f}
C
2
m
+
2
{\displaystyle C^{2m+2}}
T
(
h
)
=
τ
0
+
τ
1
h
2
+
τ
2
h
4
+
⋯
+
τ
m
h
2
m
+
α
m
+
1
(
h
)
h
2
m
+
2
{\displaystyle T(h)=\tau _{0}+\tau _{1}h^{2}+\tau _{2}h^{4}+\dots +\tau _{m}h^{2m}+\alpha _{m+1}(h)h^{2m+2}}
τ
0
=
∫
a
b
f
(
x
)
d
x
{\displaystyle \tau _{0}=\int _{a}^{b}f(x)dx}
τ
k
=
B
2
k
(
2
k
)
!
[
f
(
2
k
−
1
)
(
a
)
−
f
(
2
k
−
1
)
(
b
)
]
(
k
=
1
,
2
,
…
,
m
)
{\displaystyle \tau _{k}={\frac {B_{2k}}{(2k)!}}\left[f^{(2k-1)}(a)-f^{(2k-1)}(b)\right]\ \ (k=1,2,\dots ,m)}
α
m
+
1
(
h
)
=
B
2
m
+
2
(
2
m
+
2
)
!
(
b
−
a
)
f
(
2
m
−
2
)
(
ξ
(
h
)
)
{\displaystyle \alpha _{m+1}(h)={\frac {B_{2m+2}}{(2m+2)!}}(b-a)f^{(2m-2)}(\xi (h))}
が成立する[ 8] [ 9]
B
2
k
{\displaystyle B_{2k}}
ベルヌーイ数 であり、
ξ
(
h
)
{\displaystyle \xi (h)}
a
<
ξ
(
h
)
<
b
{\displaystyle a<\xi (h)<b}
h
{\displaystyle h}
h
{\displaystyle h}
[ 10]
h
{\displaystyle h}
リチャードソンの補外 を行うことで、より精度の高い定積分の評価が得られる。これがロンバーグ積分である[ 3] [ 4]
ロンバーグ積分のアルゴリズムは以下のように記述される[ 3]
いくつかの幅
h
0
{\displaystyle h_{0}}
h
1
=
h
0
/
N
1
{\displaystyle h_{1}=h_{0}/N_{1}}
h
m
=
h
0
/
N
m
{\displaystyle h_{m}=h_{0}/N_{m}}
T
0
,
0
=
T
(
h
0
)
{\displaystyle T_{0,0}=T(h_{0})}
T
m
,
0
=
T
(
h
m
)
{\displaystyle T_{m,0}=T(h_{m})}
次の漸化式（ネヴィルのアルゴリズム ）
T
i
,
k
=
T
i
,
k
−
1
+
T
i
,
k
−
1
−
T
i
−
1
,
k
−
1
[
h
i
−
k
h
i
]
2
−
1
{\displaystyle \quad T_{i,k}=T_{i,k-1}+{\frac {T_{i,k-1}-T_{i-1,k-1}}{\left[{\frac {h_{i-k}}{h_{i}}}\right]^{2}-1}}}
T
i
,
k
{\displaystyle T_{i,k}}
T
m
,
m
{\displaystyle T_{m,m}}
幅
h
0
,
⋯
,
h
m
{\displaystyle h_{0},\cdots ,h_{m}}
[ 5]
h
0
=
b
−
a
,
h
1
=
h
0
2
,
…
,
h
i
=
h
i
−
1
2
,
…
{\displaystyle h_{0}=b-a,h_{1}={\frac {h_{0}}{2}},\dots ,h_{i}={\frac {h_{i-1}}{2}},\dots }
を用いていたが、Roland Bulirsch は次数が急激に増加し計算コストが増大することを抑えるために
h
0
=
b
−
a
,
h
1
=
h
0
2
,
h
2
=
h
0
3
,
…
,
h
i
=
h
i
−
2
2
,
…
{\displaystyle h_{0}=b-a,h_{1}={\frac {h_{0}}{2}},h_{2}={\frac {h_{0}}{3}},\dots ,h_{i}={\frac {h_{i-2}}{2}},\dots }
を提案している[ 11] [ 12]
補外ステップにおいてリチャードソン補外テーブルの各項の差分に基づいて積分の収束性を判定することで、計算をどこまで進めればよいか判定することも可能である[ 13] [ 14]
ロンバーグ積分
T
i
,
k
{\displaystyle T_{i,k}}
h
i
−
k
,
…
,
h
i
{\displaystyle h_{i-k},\dots ,h_{i}}
O
(
h
i
−
k
2
h
i
−
k
+
1
2
⋯
h
i
2
)
{\displaystyle {\mathcal {O}}(h_{i-k}^{2}h_{i-k+1}^{2}\cdots h_{i}^{2})}
[ 15]
h
i
{\displaystyle h_{i}}
T
i
,
k
−
∫
a
b
f
(
x
)
d
x
=
(
b
−
a
)
h
i
−
k
2
h
i
−
k
+
1
2
⋯
h
i
2
(
−
1
)
k
B
2
k
+
2
(
2
k
+
2
)
!
f
(
2
k
+
2
)
(
ξ
)
{\displaystyle T_{i,k}-\int _{a}^{b}f(x)dx=(b-a)h_{i-k}^{2}h_{i-k+1}^{2}\cdots h_{i}^{2}{\frac {(-1)^{k}B_{2k+2}}{(2k+2)!}}f^{(2k+2)}(\xi )}
が得られている[ 15]
h
2
=
h
1
/
2
{\displaystyle h_{2}=h_{1}/2}
T
2
,
1
{\displaystyle T_{2,1}}
T
3
,
1
{\displaystyle T_{3,1}}
シンプソンの公式 によって問題の定積分を評価した値に一致する[ 12]
h
{\displaystyle h}
ニュートン・コーツの公式 において積分公式を拡張することに、漸化式を解き補外を行うことは積分の次数を上げることに対応していると解釈できる[ 12]
Press らは著書「Numerical Recipes in C」においてロンバーグ積分が異なる次数の数値積分法を組み合わせた一般的な方法であることから数値積分においてこの方法を用いることを強く推奨している[ 16] [ 17]
なお、ネヴィルのアルゴリズムはラグランジュ多項式 による補外であるが、有理関数 補外を行うこともでき、この場合ロンバーグ積分における漸化式は
T
i
,
k
=
T
i
,
k
−
1
+
T
i
,
k
−
1
−
T
i
−
1
,
k
−
1
[
h
i
−
k
h
i
]
2
[
1
−
T
i
,
k
−
1
−
T
i
−
1
,
k
−
1
T
i
,
k
−
1
−
T
i
−
1
,
k
−
2
]
−
1
{\displaystyle \quad T_{i,k}=T_{i,k-1}+{\frac {T_{i,k-1}-T_{i-1,k-1}}{\left[{\frac {h_{i-k}}{h_{i}}}\right]^{2}\left[1-{\frac {T_{i,k-1}-T_{i-1,k-1}}{T_{i,k-1}-T_{i-1,k-2}}}\right]-1}}}
へと変更される[ 18]
次の定積分をロンバーグ積分により求めることを考える[ 19]
∫
0
1
x
5
d
x
=
1
6
{\displaystyle \int _{0}^{1}x^{5}dx={\frac {1}{6}}}
積分ステップとして
h
0
=
1
{\displaystyle h_{0}=1}
h
1
=
1
/
2
{\displaystyle h_{1}=1/2}
h
2
=
1
/
4
{\displaystyle h_{2}=1/4}
ロンバーグ積分におけるリチャードソン補外[ 19]
h
0
=
1
{\displaystyle h_{0}=1}
T
0
,
0
=
0.500000
=
1
2
{\displaystyle T_{0,0}=0.500000={\frac {1}{2}}}
h
1
=
1
/
2
{\displaystyle h_{1}=1/2}
T
1
,
0
=
0.265625
=
17
64
{\displaystyle T_{1,0}=0.265625={\frac {17}{64}}}
T
1
,
1
=
0.187500
=
3
16
{\displaystyle T_{1,1}=0.187500={\frac {3}{16}}}
h
2
=
1
/
4
{\displaystyle h_{2}=1/4}
T
2
,
0
=
0.192383
=
197
1024
{\displaystyle T_{2,0}=0.192383={\frac {197}{1024}}}
T
2
,
1
=
0.167969
=
43
256
{\displaystyle T_{2,1}=0.167969={\frac {43}{256}}}
T
2
,
2
=
0.166667
=
1
6
{\displaystyle T_{2,2}=0.166667={\frac {1}{6}}}
C言語 により実装されているオープンソース の科学技術計算ライブラリGNU Scientific Library (GSL) にはロンバーグ積分を行う関数が含まれている[ 20] Python の数値解析ライブラリSciPy にはロンバーグ積分により数値積分を行う scipy.integrate.romberg および scipy.integrate.romb というふたつの関数がある[ 21] [ 22]
^ “ロンバーグ積分 ”. 2021年1月18日 閲覧。 ^ “ＲＯＭＢＧＳ／Ｄ （ロンバーグ積分） ”. 2021年1月18日 閲覧。 ^ a b c Stoer & Bulirsch, p. 161.
^ a b Press et al., p. 140.
^ a b Romberg, W. (1955). “Vereinfachte numerische Integration”. Det Kongelige Norske Videnskabers Selskab Forhandlinger (Trondheim) 28 (7): 30–36.
^ Stoer & Bulirsch, p. 148.
^ Press et al., p. 133.
^ Stoer & Bulirsch, p. 160.
^ Press et al., p. 141-142.
^ Stoer & Bulirsch, pp. 160-161.
^ Bulirsch, Roland (1964). “Bemerkungen zur Romberg-Integration”. Numerische Mathematik 6 (1): 6–16. doi :10.1007/BF01386048 . ISSN 0029-599X . ^ a b c Stoer & Bulirsch, p. 163.
^ Stoer & Bulirsch, p. 164.
^ Dahlquist, Germund; Björck, Åke (2008). Numerical Methods in Scientific Computing, Volume I . Society for Industrial and Applied Mathematics. p. 551. doi :10.1137/1.9780898717785 . ISBN 978-0-89871-644-3 ^ a b Stoer & Bulirsch, p. 165.
^ Press et al., p. 130.
^ Press et al., pp. 140-141.
^ Stoer & Bulirsch, pp. 164-165.
^ a b Stoer & Bulirsch, p. 162.
^ “Numerical Integration - GSL 2.6 documentation ”. 2021年1月18日 閲覧。 ^ “scipy.integrate.romberg - SciPy v1.6.0 Reference Guide ”. 2021年1月18日 閲覧。 ^ “scipy.integrate.romb - SciPy v1.6.0 Reference Guide ”. 2021年1月18日 閲覧。
![{\displaystyle [a,b]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9c4b788fc5c637e26ee98b45f89a5c08c85f7935)
![{\displaystyle [x_{i},x_{i+1}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7138606cdb8eb7dddaba59b5aefe1ade6bc05a1b)
![{\displaystyle T(h)=h\left[{\frac {1}{2}}f(x_{0})+f(x_{1})+\dots +f(x_{N-1})+{\frac {1}{2}}f(x_{N})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1b5ac3c0380af4819b4f5a336ebedafcf04abdd2)
![{\displaystyle \tau _{k}={\frac {B_{2k}}{(2k)!}}\left[f^{(2k-1)}(a)-f^{(2k-1)}(b)\right]\ \ (k=1,2,\dots ,m)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/39502243f73ae5ec267d1424eb56272393814ac8)
![{\displaystyle \quad T_{i,k}=T_{i,k-1}+{\frac {T_{i,k-1}-T_{i-1,k-1}}{\left[{\frac {h_{i-k}}{h_{i}}}\right]^{2}-1}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/843267bf869aaab8970dadf838fa04238986d913)
![{\displaystyle \quad T_{i,k}=T_{i,k-1}+{\frac {T_{i,k-1}-T_{i-1,k-1}}{\left[{\frac {h_{i-k}}{h_{i}}}\right]^{2}\left[1-{\frac {T_{i,k-1}-T_{i-1,k-1}}{T_{i,k-1}-T_{i-1,k-2}}}\right]-1}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cf8277dbe890116b317ea2309891d15d2187377d)
