Romberg 法による数値積分

（ここに空行を入れる）

プログラムの動作

integ_romberg.py は、ロンバーグ積分法（Romberg extrapolation method）を用いて、指定された区間における関数の定積分を数値的に計算するPythonプログラムです。

このプログラムの主な機能は以下の通りです。

  • 積分対象関数の定義: デフォルトで \(f(x) = e^x\) を積分対象としています。

  • 解析解との比較: 積分対象関数の解析的な積分（原始関数）も定義されており、数値計算結果との比較に使用されます。デフォルトでは \(F(x) = e^x\) です。

  • ロンバーグ積分の実行: 台形則（Trapezoidal rule）をベースに、リチャードソンの外挿法（Richardson extrapolation）を繰り返し適用することで、数値積分の精度を向上させます。

  • 収束判定: 計算された積分値の差が、事前に設定された許容誤差 eps より小さくなった場合に収束と判断し、計算を終了します。

  • 進捗表示: 各ステップでの台形則の結果 (\(S_{k,0}\)) およびロンバーグ外挿による高精度な近似値 (\(S_{k,d}\)) を標準出力に逐次表示します。

  • パラメータ: 積分区間（デフォルトで \(1.0\) から \(2.0\)）、最大反復回数（nhmax）、収束判定の許容誤差（eps）は、プログラム内で定義されています。

このプログラムは、解析的な方法で積分が困難な関数や、高速かつ高精度な数値積分が必要な場合に、その精度と効率を実演・検証するために利用できます。

原理

integ_romberg.py は、ロンバーグ積分法に基づいて定積分を計算します。ロンバーグ積分法は、台形則による数値積分の結果をリチャードソンの外挿法で補正し、より高精度な近似を得る手法です。

  1. 台形則（Trapezoidal Rule）: 最初に、積分区間 \([a, b]\) を \(n\) 個の小区間に分割し、各小区間で関数を直線で近似する台形則を用いて積分値を計算します。ステップ幅を \(h = (b-a)/n\) とすると、台形則 \(T_n\) は以下の式で表されます。 $\(T_n = \frac{h}{2} [f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)]\)\( プログラムでは、`S[k][0]` が \)2^k\( 分割での台形則の結果に対応します。具体的には、各小区間 \)[x_i, x_{i+1}]\( での台形面積 \)\frac{f(x_i) + f(x_{i+1})}{2} h$ を合計することで計算されます。

  2. ロンバーグ外挿法（Romberg Extrapolation）: 台形則の計算結果 \(T_n\)（プログラムでは \(S_{k,0}\) と表記）の誤差は、刻み幅 \(h\) の偶数べき乗で表される漸近展開を持ちます。この性質を利用し、異なる刻み幅で得られた台形則の結果を線形結合することで、より高次の誤差項を打ち消し、収束を加速させます。 ロンバーグ積分の結果は、以下の漸化式で計算されます。 $\(S_{k,d} = \frac{4^d S_{k,d-1} - S_{k-1,d-1}}{4^d - 1}\)$ ここで、

    • \(k\) は分割回数を表し、ステップが進むごとに積分区間の分割数が \(2^k\) に増加します。

    • \(d\) は外挿の次数を表します。\(d=0\) の場合は単純な台形則の結果（\(S_{k,0}\)）であり、\(d\) が増加するにつれて精度が向上します。

    • \(S_{k,d}\) は \(O(h^{2(d+1)})\) の精度を持つとされます。

プログラム中の S[k][d] は、この \(S_{k,d}\) に対応します。S は2次元配列として実装され、S[k][0] に台形則の結果が格納され、S[k][d] に外挿された結果が格納されます。

  1. 収束条件: プログラムは、現在の最高精度のロンバーグ近似値 \(S_{k,k}\) と前回の最高精度近似値 \(S_{k-1, k-1}\) の差の絶対値が、ユーザー定義の許容誤差 eps より小さくなったときに収束したと判断し、計算を終了します。 $\(|S_{k,k} - S_{k-1,k-1}| < \text{eps}\)$

必要な非標準ライブラリとインストール方法

integ_romberg.py は、以下の非標準ライブラリを使用します。

  • numpy: 数値計算を効率的に行うためのライブラリです。本プログラムでは直接的な配列操作には標準のリストのリストを使用していますが、数値計算プロジェクトで一般的に推奨されるライブラリです。

インストール方法: コマンドラインまたはターミナルで以下のコマンドを実行してインストールしてください。

pip install numpy

なお、math モジュール（exp 関数）はPythonの標準ライブラリのため、別途インストールは不要です。

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。 積分する関数 func(x)、その解析解 F(x)、積分区間 xmin と xmax、最大反復回数 nhmax、および収束判定の許容誤差 eps は、すべてプログラムのソースコード内に直接記述（ハードコード）されています。

これらの設定値を変更する場合は、integ_romberg.py のソースコードをテキストエディタで開き、func 関数、F 関数、およびデフォルトのパラメータ群（nhmax_def, eps_def, xmin_def, xmax_def）の値を直接編集してください。

生成される出力ファイル

integ_romberg.py は、実行時に明示的に出力ファイルを生成しません。 すべての計算結果およびプログラムの進行状況は、標準出力（STDOUT）にテキスト形式で表示されます。

出力内容をファイルに保存したい場合は、シェルが提供するリダイレクト機能を使用してください。

例: プログラムの実行結果を romberg_output.txt というファイルに保存する場合。

python integ_romberg.py > romberg_output.txt

コマンドラインでの使用例 (Usage)

このプログラムは、コマンドライン引数を受け付けません。 実行するには、Pythonインタプリタにスクリプトファイルを指定するだけです。

python integ_romberg.py

コマンドラインでの具体的な使用例

以下のコマンドは、デフォルトのパラメータ（\(f(x) = e^x\) を区間 \([1.0, 2.0]\) で積分、許容誤差 \(1.0 \times 10^{-6}\)、最大反復回数 \(15\)）でプログラムを実行します。

python integ_romberg.py

実行結果の例:

Numerical integration using by Romberg extrapolation method

eps for convergnece: 1e-06
Analytical values:
  f(1.0)=2.718281828459045
  integ(f)[1.0, 2.0]=4.670774270471604

S00 = 5.43656365691809
  S[1][0] = 4.881775796013629
     S[1][1] = 4.6821798363487055
  diff = S[1][1] - S[0][0] = 4.6821798363487055 - 5.43656365691809 = -0.7543838205693845
  S[2][0] = 4.729119934185792
     S[2][1] = 4.670908006240245
     S[2][2] = 4.670775080271701
  diff = S[2][2] - S[1][1] = 4.670775080271701 - 4.6821798363487055 = -0.011404756077004489
  S[3][0] = 4.680608034515568
     S[3][1] = 4.670781067711831
     S[3][2] = 4.670774270597332
     S[3][3] = 4.670774270471607
  diff = S[3][3] - S[2][2] = 4.670774270471607 - 4.670775080271701 = -8.098000000000005e-07

出力の説明:

  • eps for convergnece: 1e-06: 収束判定に使用される許容誤差が表示されます。

  • Analytical values:: 積分対象の関数 \(f(x)\) の \(x=\text{xmin}\) における値と、区間 \([\text{xmin}, \text{xmax}]\) における解析解による定積分値が表示されます。

  • S00 = ...: \(S_{0,0}\) は、最も粗い分割（区間を1つと見なす）での台形則による初期近似値です。

  • S[k][0] = ...: 各ステップ k における台形則の結果です。k が増えるごとに分割数が \(2^k\) に倍増し、より精度の高い台形則の結果が得られます。

  • S[k][d] = ...: 各ステップ k におけるロンバーグ外挿の結果です。d が増えるほど、より高次の外挿が適用され、精度が向上します。S[k][k] がそのステップで得られる最も精度の高いロンバーグ近似値です。

  • diff = S[k][k] - S[k-1][k-1] = ...: 現在のステップで得られた最高精度近似値と、前回のステップで得られた最高精度近似値との差を示します。この diff の絶対値が eps より小さくなると、プログラムは収束したと判断し、計算を終了します。上記の例では、\(k=3\) のステップで diff の絶対値が eps (1e-06) より小さくなったため、ここで計算が終了します。最終的な S[3][3] の値が解析値に非常に近いことが確認できます。