integ_scipy_quad.py 技術ドキュメント

プログラムの動作

integ_scipy_quad.py は、Pythonの数値計算ライブラリ SciPy およびグラフ描画ライブラリ Matplotlib を用いて、様々な数値積分手法をデモンストレーションするプログラムです。

主な機能は以下の通りです。

定積分:
- 関数 \(f(x) = \frac{4}{1+x^2}\) を \(0\) から \(1\) まで積分し、円周率 \(\pi\) の近似値を計算します。これには scipy.integrate.quad (QUADPACKルーチン) と scipy.integrate.romberg (ロンバーグ積分) の2つの手法を使用します。
広義積分:
- 関数 \(f(x) = x^5 e^{-x} \sin(x)\) を \(0\) から無限大 (\(\infty\)) まで積分します。これも scipy.integrate.quad を使用します。
離散データの累積積分:
- 関数 \(f(x) = x^4\) の離散データに対して、scipy.integrate.cumulative_trapezoid を用いて累積積分を計算します。
- 計算結果は解析解 \(\frac{x^5}{5}\) と比較され、Matplotlib を使用してグラフとして可視化されます。

このプログラムは、SciPy が提供する強力な数値積分ツール群の利用方法を具体的に示すことで、科学技術計算における積分の課題解決に貢献します。

原理

このプログラムでは、以下の数値積分手法とその数学的原理が利用されています。

1. QUADPACKによる定積分 (`scipy.integrate.quad`)

scipy.integrate.quad は、適応型求積法（Adaptive Quadrature）に基づいています。これは、積分区間を動的に分割し、各部分区間でガウス求積法などの高精度な求積規則を適用することで、高い精度と効率を実現します。ユーザーが指定した被積分関数、積分区間、およびオプションのパラメータ（許容誤差など）に基づいて、積分値とその推定誤差を返します。円周率 \(\pi\) の計算には、以下の積分式を利用します。

\[\pi = \int_0^1 \frac{4}{1+x^2} dx\]

この積分は、逆正接関数 \(\arctan(x)\) の定義から導かれます。

2. ロンバーグ積分 (`scipy.integrate.romberg`)

ロンバーグ積分は、台形公式に基づき、リチャードソン補外を繰り返し適用することで、積分の精度を高める数値積分法です。異なるステップサイズで計算された台形公式の結果から、より高次の誤差項をキャンセルすることで、収束を加速させます。これにより、少ない関数評価回数で高い精度を得られる場合があります。

3. 台形公式による累積積分 (`scipy.integrate.cumulative_trapezoid`)

台形公式は、数値積分の中でも最も基本的な手法の一つです。区間 \([a, b]\) における関数 \(f(x)\) の積分を、この区間を複数の小区間に分割し、各小区間において関数を直線で近似（台形）することで計算します。隣接する2点 \(x_i, x_{i+1}\) における台形公式の面積 \(A_i\) は以下で与えられます。

\[A_i = \frac{f(x_i) + f(x_{i+1})}{2} (x_{i+1} - x_i)\]

cumulative_trapezoid は、この台形公式を連続的に適用し、各点までの積分の累積値を計算します。プログラムでは、関数 \(f(x) = x^4\) の累積積分を計算し、その解析解と比較します。\(x^4\) の不定積分は以下の通りです。

\[\int x^4 dx = \frac{x^5}{5} + C\]

ここで \(C\) は積分定数です。initial=0 オプションにより、\(x=0\) での積分値が \(0\) に設定されるため、\(C=0\) となり、解析解は \(\frac{x^5}{5}\) となります。

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

本プログラムの実行には、以下のPython非標準ライブラリが必要です。

NumPy: 数値計算を効率的に行うための基盤ライブラリ。
SciPy: 科学技術計算のツールを提供するライブラリ。特に integrate モジュールを使用します。
Matplotlib: グラフ描画を行うためのライブラリ。

これらのライブラリは、pip コマンドを使用してインストールできます。

pip install numpy scipy matplotlib

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。すべてのデータ（被積分関数、積分区間、離散データ）はプログラムのソースコード内で直接定義され、生成されます。

生成される出力ファイル

本プログラムは、実行時に直接ファイルを生成して保存することはありません。出力は主に以下の2種類です。

標準出力 (コンソール):
- integrate.quad および integrate.romberg による各積分結果が、(積分値, 推定誤差) のタプル形式でコンソールに出力されます。
グラフィカル出力 (Matplotlibウィンドウ):
- 離散データに対する cumulative_trapezoid による累積積分結果と、その解析解である \(x^5/5\) の比較グラフが Matplotlib のウィンドウに表示されます。このウィンドウはユーザーが手動で閉じるまで表示され続けます。

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

本プログラムは引数を取らずに直接実行されます。

python integ_scipy_quad.py

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

以下のコマンドを実行することで、プログラムが動作し、数値積分結果とグラフが出力されます。

python integ_scipy_quad.py

実行すると、まずコンソールに以下の出力が表示されます。

(res, error)= (3.141592653589793, 3.487868498007664e-14)
(res, error)= 3.141592653589793
(res, error)= (7.986663249053229, 3.844391295325852e-13)

それぞれの出力の意味は以下の通りです。

1行目: integrate.quad を用いた円周率 \(\pi\) の計算結果（積分値と推定誤差）。
2行目: integrate.romberg を用いた円周率 \(\pi\) の計算結果（積分値のみ）。
3行目: integrate.quad を用いた広義積分 \(\int_0^\infty x^5 e^{-x} \sin(x) dx\) の計算結果（積分値と推定誤差）。

コンソール出力に続き、Matplotlib によって「Comparison of Cumulative Integration」というタイトルのグラフウィンドウが表示されます。このグラフは、関数 \(f(x)=x^4\) の離散データに対するcumulative_trapezoidによる累積積分結果（赤い丸）と、解析解 \(x^5/5\)（青い線）を比較したものです。両者がほぼ完全に一致していることを確認できます。