ロンバーグ積分

ロンバーグ積分[1][2]（ロンバーグせきぶん、Romberg integration[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,\cdots ,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})+\cdots +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}+\cdots +\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,\cdots ,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]。

  1. いくつかの幅 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})} を求める.
  2. 次の漸化式（ネヴィルのアルゴリズム）
    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}},\ \ \cdots ,h_{i}={\frac {h_{i-1}}{2}},\cdots }

を用いていたが、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}},\ \ \cdots ,\ \ h_{i}={\frac {h_{i-2}}{2}},\cdots }

を提案している[11][12]。

補外ステップにおいてリチャードソン補外テーブルの各項の差分に基づいて積分の収束性を判定することで、計算をどこまで進めればよいか判定することも可能である[13][14]。

特徴

ロンバーグ積分 T i , k {\displaystyle T_{i,k}} の誤差は、評価に用いたステップ h i − k , ⋯ , h i {\displaystyle h_{i-k},\cdots ,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].

脚注

  1. ^ “ロンバーグ積分”. 2021年1月18日閲覧。
  2. ^ “ＲＯＭＢＧＳ／Ｄ　（ロンバーグ積分）”. 2021年1月18日閲覧。
  3. ^ a b c Stoer & Bulirsch, p. 161.
  4. ^ a b Press et al., p. 140.
  5. ^ a b Romberg, W. (1955). “Vereinfachte numerische Integration”. Det Kongelige Norske Videnskabers Selskab Forhandlinger (Trondheim) 28 (7): 30–36. 
  6. ^ Stoer & Bulirsch, p. 148.
  7. ^ Press et al., p. 133.
  8. ^ Stoer & Bulirsch, p. 160.
  9. ^ Press et al., p. 141-142.
  10. ^ Stoer & Bulirsch, pp. 160-161.
  11. ^ Bulirsch, Roland (1964). “Bemerkungen zur Romberg-Integration”. Numerische Mathematik 6 (1): 6–16. doi:10.1007/BF01386048. ISSN 0029-599X. 
  12. ^ a b c Stoer & Bulirsch, p. 163.
  13. ^ Stoer & Bulirsch, p. 164.
  14. ^ 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 
  15. ^ a b Stoer & Bulirsch, p. 165.
  16. ^ Press et al., p. 130.
  17. ^ Press et al., pp. 140-141.
  18. ^ Stoer & Bulirsch, pp. 164-165.
  19. ^ a b Stoer & Bulirsch, p. 162.
  20. ^ “Numerical Integration - GSL 2.6 documentation”. 2021年1月18日閲覧。
  21. ^ “scipy.integrate.romberg - SciPy v1.6.0 Reference Guide”. 2021年1月18日閲覧。
  22. ^ “scipy.integrate.romb - SciPy v1.6.0 Reference Guide”. 2021年1月18日閲覧。

参考文献

  • Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (1992). Numerical Recipes in C: The Art of Scientific Computing (2nd ed.). Cambridge University Press. doi:10.2277/0521431085. ISBN 978-0-521-43108-8 
  • Stoer, Josef; Bulirsch, R. (2002). Introduction to Numerical Analysis. Springer. doi:10.1007/978-0-387-21738-3. ISBN 978-0-387-21738-3 

関連項目