"""
ロンバーグ補外法を用いた数値積分を実行するスクリプト。
このスクリプトは、指定された関数の定積分をロンバーグ補外法を用いて近似します。
台形則による近似を繰り返し計算し、ロンバーグ補外によってその精度を高めます。
収束基準 (eps) が満たされるか、最大繰り返し回数 (nhmax) に達するまで計算を継続します。
:doc:`integ_romberg_usage`
"""
import numpy as np
from math import exp
# define function to be integrated
[ドキュメント]
def func(x):
"""
積分対象の関数を定義します。
この関数は、数値積分で近似される対象の関数 `f(x)` を定義します。
現在の実装では `e^x` を返します。
:param x: float: 関数の入力値。
:returns: float: xにおける関数の値。
"""
return exp(x)
# define analytical integration of f(x)
[ドキュメント]
def F(x):
"""
積分対象の関数の解析的な不定積分を定義します。
この関数は、`func(x)` の解析的な不定積分 `F(x)` を定義します。
現在の実装では `e^x` の不定積分である `e^x` を返します。
数値積分の結果と比較するために使用されます。
:param x: float: 不定積分の入力値。
:returns: float: xにおける不定積分の値。
"""
return exp(x)
#===================
# parameters (Default values)
#===================
nhmax_def = 15
eps_def = 1.0e-6
xmin_def, xmax_def = 1.0, 2.0
[ドキュメント]
def main():
"""
ロンバーグ補外法による数値積分を実行するメイン関数です。
この関数は、`nhmax_def`, `eps_def`, `xmin_def`, `xmax_def` で定義された
デフォルトパラメータを使用して、`func(x)` の `xmin` から `xmax` までの
定積分をロンバーグ補外法で計算します。
計算は以下の手順で行われます。
1. 初期区間での台形近似 `S[0][0]` を計算します。
2. 反復ごとに区間を2倍に分割し、新たな台形近似 `S[k][0]` を計算します。
3. `S[k][0]` と以前の近似値 `S[k-1][d-1]` を用いてロンバーグ補外 `S[k][d]` を計算します。
4. 連続するロンバーグ近似値 `S[k][k]` と `S[k-1][k-1]` の差が `eps` より小さくなったら収束とみなし、
計算を終了します。
5. 最大反復回数 `nhmax` に達しても収束しない場合は、その旨を通知します。
計算の進捗と結果は標準出力に表示されます。
:returns: None
"""
# パラメータのローカル変数化
nhmax = nhmax_def
eps = eps_def
xmin, xmax = xmin_def, xmax_def
print("Numerical integration using by Romberg extrapolation method")
fx = func(xmin)
Fx = F(xmax) - F(xmin)
print("")
print("eps for convergnece: {}".format(eps))
print("Analytical values:")
print(" f({})={}".format(xmin, fx))
print(" integ(f)[{}, {}]={}".format(xmin, xmax, Fx))
print("")
# create 2D array by specifing the size to be nhmax * nhmax
S = []
for i in range(nhmax + 1): # インデックス範囲外エラーを避けるためサイズを調整
S.append([])
for j in range(nhmax + 1):
S[i].append(0.0)
n = 1
h0 = xmax - xmin
y0 = func(xmin)
y1 = func(xmax)
S[0][0] = (y0 + y1) / 2.0 * h0
print("S00 = {}".format(S[0][0]))
kfin = 0
for k in range(1, nhmax + 1):
n *= 2
h = h0 / n
# 台形則の計算
S[k][0] = 0.0
for i in range(n):
xi = xmin + i * h
xi1 = xi + h
yi = func(xi)
yi1 = func(xi1)
S[k][0] += (yi + yi1) / 2.0 * h
print(" S[{}][0] = {}".format(k, S[k][0]))
# 外挿の計算
for d in range(1, k + 1):
S[k][d] = (4.0**d * S[k][d-1] - S[k-1][d-1]) / (4.0**d - 1.0)
print(" S[{}][{}] = {}".format(k, d, S[k][d]))
diff = S[k][k] - S[k-1][k-1]
print(" diff = S[{}][{}] - S[{}][{}] = {} - {} = {}"
.format(k, k, k-1, k-1, S[k][k], S[k-1][k-1], diff))
if abs(diff) < eps:
kfin = k
break
else:
print("Not converged within {} iterations".format(nhmax))
"""
元コードのコメントに含まれる別実装案やPerl版のロジックもそのまま保持します。
"""
if __name__ == "__main__":
main()