ロンバーグ積分のソースを表示

'''ロンバーグ積分'''<ref>{{Cite web|和書|url=http://www.ipc.tohoku-gakuin.ac.jp/nken/java/romberg/romberg.html |title=ロンバーグ積分 |accessdate=2021-01-18}}</ref><ref>{{Cite web|和書|url=http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/xml/193.xml |title=ＲＯＭＢＧＳ／Ｄ　（ロンバーグ積分） |accessdate=2021-01-18}}</ref>（ロンバーグせきぶん、Romberg integration<ref name="#1">Stoer & Bulirsch, p. 161.</ref><ref name="#2">Press et al., p. 140.</ref>）または'''ロンベルク積分'''は、[[関数 (数学)|関数]]の[[数値積分]]アルゴリズムのひとつである。この方法では[[台形公式]]と[[リチャードソンの補外]]を組み合わせ離散化幅 <math>h</math> をゼロとする極限として数値積分を評価する。他の数値積分法に比べ、少ない回数の被積分関数の評価によって高精度の結果が得られる。

1955年に{{仮リンク|ヴェルナー・ロンベルク|en|Werner Romberg}}によって考案された<ref name="Romberg1955"/>。

== 概要 ==
ロンバーグ積分は区間 <math>[ a, b ]</math> で定義された関数 <math>f</math> の定積分
{{Indent|<math>\int_a^b f ( x ) dx</math>}}
を数値的に求めるアルゴリズムである。積分区間 <math>[ a, b ]</math> を幅 <math>h = ( b - a ) / N</math> の <math>N</math> 個の小区間 <math>[ x_i, x_{i+1} ]</math> (<math>x_i = a + i h</math>, <math>i = 0, 1, \dots, N-1</math>) に分割するとき,
数値積分アルゴリズムのひとつである[[台形公式]]は求める定積分を
{{Indent|<math>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]</math>}}
により近似する<ref>Stoer & Bulirsch, p. 148.</ref><ref>Press et al., p. 133.</ref>。台形公式による値 <math>T ( h )</math> と真の定積分の値の差は <math>\mathcal{O} ( h^2 )</math> であるが、オイラー・マクローリンの和公式（[[オイラーの和公式]]）によりそれをより詳しく評価することができる。すなわち、<math>f</math> が <math>C^{2m+2}</math>-級関数であるならば
{{Indent|<math>T ( h ) = \tau_0 + \tau_1 h^2 + \tau_2 h^4 + \dots + \tau_m h^{2m} + \alpha_{m+1} ( h ) h^{2m+2}</math>}}
{{Indent2|<math>\tau_0 = \int_a^b f ( x ) dx</math>}}
{{Indent2|<math>\tau_k = \frac{ B_{2k} }{ ( 2 k )! } \left[ f^{(2k-1)} ( a ) - f^{(2k-1)} ( b ) \right] \ \ (k = 1, 2, \dots, m)</math>}}
{{Indent2|<math>\alpha_{m+1} ( h ) = \frac{ B_{2m+2} }{ ( 2 m + 2 )! } ( b - a ) f^{(2 m - 2)} ( \xi ( h ) )</math>}}
が成立する<ref>Stoer & Bulirsch, p. 160.</ref><ref>Press et al., p. 141-142.</ref>。ここに <math>B_{2k}</math> は[[ベルヌーイ数]] であり、<math>\xi ( h )</math> は <math>a < \xi ( h ) < b</math> を満足するある数である。それ故に台形公式を <math>h</math> に関する漸近級数と見るとき、誤差項としては <math>h</math> の偶数次の項だけが現れる<ref>Stoer & Bulirsch, pp. 160-161.</ref>。そこで複数の <math>h</math> について台形公式を適用し、それらの結果について誤差項を打ち消すように[[リチャードソンの補外]]を行うことで、より精度の高い定積分の評価が得られる。これがロンバーグ積分である<ref name="#1"/><ref name="#2"/>。

== アルゴリズム ==
ロンバーグ積分のアルゴリズムは以下のように記述される<ref name="#1"/>。
# いくつかの幅 <math>h_0</math>, <math>h_1 = h_0 / N_1</math>, ..., <math>h_m = h_0 / N_m</math> に対して台形公式を適用することにより定積分の近似値 <math>T_{0, 0} = T ( h_0 )</math>, ..., <math>T_{m, 0} = T ( h_m )</math> を求める.
# 次の漸化式（[[ネヴィルのアルゴリズム]]）<br/><math>\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 }</math><br/>により帰納的に <math>T_{i, k}</math> を計算し、得られた <math>T_{m, m}</math> が求める補外値を与える。

幅 <math>h_0, \cdots, h_m</math> としては、ロンバーグ<ref name="Romberg1955"/>は
{{Indent|<math>h_0 = b - a, h_1 = \frac{ h_0 }{ 2 } , \dots, h_i = \frac{ h_{i-1} }{ 2 } , \dots</math>}}
を用いていたが、[[:en:Roland Bulirsch|Roland Bulirsch]]は次数が急激に増加し計算コストが増大することを抑えるために
{{Indent|<math>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</math>}}
を提案している<ref name="Bulirsch1964"/><ref name="#3">Stoer & Bulirsch, p. 163.</ref>。

補外ステップにおいてリチャードソン補外テーブルの各項の差分に基づいて積分の収束性を判定することで、計算をどこまで進めればよいか判定することも可能である<ref>Stoer & Bulirsch, p. 164.</ref><ref>{{Cite book |first1=Germund |last1=Dahlquist |first2=Åke |last2=Björck |title=Numerical Methods in Scientific Computing, Volume I |publisher=Society for Industrial and Applied Mathematics |date=2008 |isbn=978-0-89871-644-3 |doi=10.1137/1.9780898717785 |page=551}}</ref>。

== 特徴 ==
ロンバーグ積分 <math>T_{i, k}</math> の誤差は、評価に用いたステップ <math>h_{i-k}, \dots, h_i</math> から <math>\mathcal{O} ( h_{i-k}^2 h_{i-k+1}^2 \cdots h_i^2 )</math> と見積もられる<ref name="#4">Stoer & Bulirsch, p. 165.</ref>。特に <math>h_i</math> をロンバーグの提案に従って選ぶときには誤差の厳密な評価
{{Indent|<math>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_{2 k + 2} }{ ( 2 k + 2 )! } f^{(2 k + 2)} ( \xi )</math>}}
が得られている<ref name="#4"/>。

<math>h_2 = h_1 / 2</math> に選ぶとき、リチャードソン補外において現れる <math>T_{2, 1}</math> や <math>T_{3, 1}</math> は[[シンプソンの公式]]によって問題の定積分を評価した値に一致する<ref name="#3"/>。その意味で、より多くの積分幅 <math>h</math> での評価を行うことは[[ニュートン・コーツの公式]]において積分公式を拡張することに、漸化式を解き補外を行うことは積分の次数を上げることに対応していると解釈できる<ref name="#3"/>。

Press らは著書「Numerical Recipes in C」においてロンバーグ積分が異なる次数の数値積分法を組み合わせた一般的な方法であることから数値積分においてこの方法を用いることを強く推奨している<ref>Press et al., p. 130.</ref>。特に、被積分関数が積分区間において十分に滑らかであれば、ロンバーグ積分はニュートン・コーツ型の積分公式に比べて極めて少ない回数の被積分関数の評価だけで十分な精度の結果を与える<ref>Press et al., pp. 140-141.</ref>。

なお、ネヴィルのアルゴリズムは[[ラグランジュ多項式]]による補外であるが、[[有理関数]]補外を行うこともでき、この場合ロンバーグ積分における漸化式は
{{Indent|<math>\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 }</math>}}
へと変更される<ref>Stoer & Bulirsch, pp. 164-165.</ref>。

== 具体例 ==
次の定積分をロンバーグ積分により求めることを考える<ref name="#5">Stoer & Bulirsch, p. 162.</ref>。
{{Indent|<math>\int_0^1 x^5 dx = \frac{ 1 }{ 6 }</math>}}
積分ステップとして <math>h_0 = 1</math>, <math>h_1 = 1/2</math>, <math>h_2 = 1/4</math> を取り補外を行う。その結果、次の表に示すようにロンバーグ積分によって正確な積分値を得ることができる。

{| class="wikitable"
|+ ロンバーグ積分におけるリチャードソン補外<ref name="#5"/>
! <math>h_0 = 1</math>
| <math>T_{0, 0} = 0.500 000 = \frac{1}{2}</math> || || 
|-
! <math>h_1 = 1/2</math>
| <math>T_{1, 0} = 0.265 625 = \frac{17}{64}</math> || <math>T_{1, 1} = 0.187 500 = \frac{3}{16}</math> || 
|-
! <math>h_2 = 1/4</math>
| <math>T_{2, 0} = 0.192 383 = \frac{197}{1024}</math> || <math>T_{2, 1} = 0.167 969 = \frac{43}{256}</math> || <math>T_{2, 2} = 0.166 667 = \frac{1}{6}</math>
|}

== ライブラリ ==
[[C言語]]により実装されている[[オープンソース]]の科学技術計算ライブラリ[[GNU Scientific Library]] (GSL) にはロンバーグ積分を行う関数が含まれている<ref>{{Cite web |url=https://www.gnu.org/software/gsl/doc/html/integration.html |title=Numerical Integration - GSL 2.6 documentation |accessdate=2021-01-18}}</ref>。また、[[Python]]の数値解析ライブラリ[[SciPy]]にはロンバーグ積分により数値積分を行う <code>scipy.integrate.romberg</code> および <code>scipy.integrate.romb</code> というふたつの関数がある<ref>{{Cite web |url=https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html |title=scipy.integrate.romberg - SciPy v1.6.0 Reference Guide |accessdate=2021-01-18}}</ref><ref>{{Cite web |url=https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romb.html |title=scipy.integrate.romb - SciPy v1.6.0 Reference Guide |accessdate=2021-01-18}}</ref>.

== 脚注 ==
{{Reflist |2 |refs=
<ref name="Romberg1955">{{Cite journal|last1=Romberg|first1=W.|title=Vereinfachte numerische Integration|journal=Det Kongelige Norske Videnskabers Selskab Forhandlinger|volume=28|year=1955|location=Trondheim|pages=30–36|issue=7}}</ref>
<ref name="Bulirsch1964">{{cite journal|last1=Bulirsch|first1=Roland|title=Bemerkungen zur Romberg-Integration|journal=Numerische Mathematik|volume=6|issue=1|year=1964|pages=6–16|issn=0029-599X|doi=10.1007/BF01386048}}</ref>
}}

== 参考文献 ==
* {{Cite book |last1=Press |first1=William H. |first2=Saul A. |last2=Teukolsky |first3=William T. |last3=Vetterling |first4=Brian P. |last4=Flannery |date=1992 |title=Numerical Recipes in C: The Art of Scientific Computing |edition=2nd |publisher=Cambridge University Press |doi=10.2277/0521431085 |isbn=978-0-521-43108-8}}
* {{Cite book |last1=Stoer |first1=Josef |last2=Bulirsch |first2=R. |title=Introduction to Numerical Analysis |publisher=Springer |date=2002 |isbn=978-0-387-21738-3 |doi=10.1007/978-0-387-21738-3}}

== 関連項目 ==
* [[台形公式]]
* [[シンプソンの公式]]
* [[ニュートン・コーツの公式]]
* {{仮リンク|ブリルシュ・ストア法|en|Bulirsch–Stoer algorithm}} - [[常微分方程式の数値解法]]のひとつで、[[修正中点法]]の結果に対してロンバーグ積分と同様にリチャードソンの補外を行う。

{{DEFAULTSORT:ろんはあくせきふん}}
[[Category:数値積分]]
[[Category:数学のエポニム]]
[[Category:数学に関する記事]]