プログラム integ_order_h.py の技術ドキュメント
プログラムの動作
integ_order_h.py は、指定された数学関数に対し、異なる数値積分手法（リーマン和、台形公式、シンプソンの公式、ボーデの公式）を用いて積分を実行し、その精度を評価するPythonプログラムです。積分範囲の分割数 h を変化させながら各手法による積分値と解析解との絶対誤差を計算し、その結果を標準出力に表示します。また、matplotlib ライブラリを使用して、積分値と誤差のグラフを生成し、数値積分の精度が h にどのように依存するかを視覚的に示します。
このプログラムは、数値積分アルゴリズムの挙動と収束性、特にh（積分刻み幅）に対する誤差の次数を理解することを目的としています。
原理
本プログラムでは、以下の数学関数とそれらの解析解（不定積分）が使用されます。
Lorentz関数:
関数: \(f(x) = \frac{1}{1 + x^2}\)
解析積分: \(\int f(x) dx = \arctan(x)\)
指数関数 (Exponential function):
関数: \(f(x) = e^x\)
解析積分: \(\int f(x) dx = e^x\)
ガウス関数 (Gaussian function):
関数: \(f(x) = e^{-x^2}\)
解析積分: \(\int f(x) dx = \frac{\sqrt{\pi}}{2} \text{erf}(x)\)
ここで \(\text{erf}(x)\) は誤差関数 (Error Function) を表します。
これらの関数に対して、以下の数値積分手法を適用します。積分区間 \([a, b]\) を \(N\) 個の小区間 \(h = (b-a)/N\) に分割し、\(x_i = a + i h\) とします。
リーマン和 (左端点): 各小区間で関数値を左端点で評価し、矩形の面積の和で近似します。 $\( \int_a^b f(x) dx \approx \sum_{i=0}^{N-1} f(x_i) h \)$
台形公式: 各小区間で関数値を両端点で評価し、台形の面積の和で近似します。 $\( \int_a^b f(x) dx \approx \sum_{i=0}^{N-1} \frac{f(x_i) + f(x_{i+1})}{2} h \)$
シンプソンの公式 (1/3 公式): 連続する2つの小区間（計3点）に対して、2次関数で近似して積分します。\(N\) は偶数である必要があります。 $\( \int_{x_i}^{x_{i+2}} f(x) dx \approx \frac{h}{3} (f(x_i) + 4f(x_{i+1}) + f(x_{i+2})) \)\( 全体の積分は、この式を \)i=0, 2, \ldots, N-2$ について合計します。
ボーデの公式: 連続する4つの小区間（計5点）に対して、4次関数で近似して積分します。\(N\) は4の倍数である必要があります。 $\( \int_{x_i}^{x_{i+4}} f(x) dx \approx \frac{2h}{45} (7f(x_i) + 32f(x_{i+1}) + 12f(x_{i+2}) + 32f(x_{i+3}) + 7f(x_{i+4})) \)\( 全体の積分は、この式を \)i=0, 4, \ldots, N-4$ について合計します。
各数値積分手法によって得られた積分値 \(S\) と解析解から得られる厳密な積分値 \(S_{ex}\) との絶対誤差 \(|S - S_{ex}|\) を計算し、h の変化に対する誤差の挙動を評価します。
必要な非標準ライブラリとインストール方法
このプログラムは、以下の非標準Pythonライブラリを使用します。
numpy: 数値計算を効率的に行うためのライブラリ。matplotlib: グラフ描画のためのライブラリ。
これらのライブラリは、pip コマンドを使ってインストールできます。
pip install numpy matplotlib
必要な入力ファイル
このプログラムは入力ファイルを必要としません。すべてのパラメータはコマンドライン引数として指定されます。
生成される出力ファイル
標準出力: 選択された関数、積分範囲、各数値積分手法の計算結果が表形式で出力されます。表には、分割数 (
nh)、刻み幅 (h)、数値積分値 (S)、絶対誤差 (error) が含まれます。グラフィカルウィンドウ:
matplotlibによって2つのグラフが描画され、別ウィンドウで表示されます。上段のグラフ: 刻み幅
hに対する各手法による積分値Sの変化を示します。解析解の直線も合わせてプロットされます。下段のグラフ: 刻み幅
hに対する各手法の絶対誤差 \(|S - S_{ex}|\) の変化を示します。両軸が対数スケールで表示され、誤差の収束次数を視覚的に確認できます。
出力ファイルは直接生成されません。グラフは表示後にユーザーがEnterキーを押すことで閉じられます。
コマンドラインでの使用例 (Usage)
プログラムは以下の形式で実行します。
python integ_order_h.py xmin xmax nhscan func_type
引数の意味は以下の通りです。
xmin: 積分範囲の開始点（浮動小数点数）。xmax: 積分範囲の終了点（浮動小数点数）。nhscan: 積分範囲の分割数を制御する整数。hのスキャンは \(2^{\text{ih}}\) のステップで実行されます。より大きな値ほど、より詳細なhのステップで計算が行われます。デフォルト値は18です。func_type: 積分対象の関数タイプを指定する文字列。以下のいずれかを指定します。lorentz(デフォルト): \(f(x) = 1 / (1 + x^2)\)gauss: \(f(x) = e^{-x^2}\)exp: \(f(x) = e^x\)
コマンドラインでの具体的な使用例
ここでは、Lorentz 関数を積分範囲 \(0.0\) から \(1.0\) まで、nhscan を 5 と設定して実行する例を示します。
python integ_order_h.py 0.0 1.0 5 lorentz
実行結果の説明:
上記のコマンドを実行すると、まず以下のような内容が標準出力に表示されます。
Numerical integration using different approximations
Function: lorentz
integration range: 0.0 - 1.0
Exact: 0.78539816
Rieman
nh h S error
1 1.0 0.5 0.2854
2 0.5 0.6 0.1854
4 0.25 0.7058 0.07959
8 0.125 0.7466 0.03875
16 0.0625 0.7667 0.01869
Trapezoid
nh h S error
1 1.0 0.75 0.0354
2 0.5 0.775 0.0104
4 0.25 0.7828 0.002591
8 0.125 0.7847 0.0006478
16 0.0625 0.7852 0.0001619
Simpson
nh h S error
2 0.5 0.7833 0.002098
4 0.25 0.7852 0.0001662
8 0.125 0.7854 1.03e-05
16 0.0625 0.7854 6.43e-07
Bode
nh h S error
4 0.25 0.7854 1.419e-06
8 0.125 0.7854 2.22e-08
16 0.0625 0.7854 3.47e-10
（表示される数値は計算環境やPythonのバージョンによって若干異なる場合があります。）
この出力の後、matplotlib によって生成されたグラフウィンドウが表示されます。ウィンドウには、上部に各手法の積分値と解析解の比較、下部に各手法の絶対誤差の対数プロットが表示され、h を減少させる（分割数を増やす）ことで誤差が減少する様子や、各手法の収束の速さ（誤差の傾き）を視覚的に確認することができます。
グラフウィンドウが表示された後、標準出力には Press ENTER to exit>> と表示され、Enterキーを押すまでプログラムは待機します。Enterキーを押すと、グラフウィンドウが閉じられ、プログラムが終了します。