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

プログラムの動作

integ_doubleexp.py は、倍指数関数型数値積分法（Double Exponential Integration, DE法）を用いて、与えられた関数を特定区間で数値積分するPythonスクリプトです。

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

関数定義: 積分対象となる関数 $f(x) = \sqrt{1 - x^2}$ がコード内で定義されています。
DE法による数値積分: 指定された区間 $[-1.0, 1.0]$ において、倍指数関数型変数変換を用いた台形公式により数値積分を行います。
解析解との比較: 積分対象の関数 $f(x) = \sqrt{1 - x^2}$ の区間 $[-1.0, 1.0]$ における解析解が $\pi/2$ であることを利用し、数値積分結果との誤差を評価します。
収束性の評価: 積分に用いるブロック数 (nblock) を増やしながら、数値積分結果と誤差がどのように変化するかを標準出力に表示します。

このプログラムは、特に積分区間の端点で導関数が発散するような特異性を持つ関数（本例では $x = \pm 1$ で $f'(x)$ が無限大になる）に対して、他の数値積分法よりも高い精度と収束性を持つDE法の適用例を示し、その有効性を確認することを目的としています。

原理

integ_doubleexp.py は、倍指数関数型積分法（DE法）を適用して数値積分を行います。DE法は、変数変換によって積分区間を無限区間に写像し、その上で台形公式を適用することで、端点特異性を持つ関数や高精度な積分が必要な場合に高い性能を発揮します。

本プログラムで積分する関数は $f(x) = \sqrt{1 - x^2}$ であり、積分区間は $[xmin, xmax]$（デフォルトは $[-1, 1]$）です。この積分 $I = \int_{xmin}^{xmax} f(x) dx$ は、以下の2段階の変数変換を経て計算されます。

線形変換: 最初に、積分区間 $[xmin, xmax]$ を標準的な区間 $[-1, 1]$ に線形変換します。 $x = \frac{xmax - xmin}{2} x_{norm} + \frac{xmax + xmin}{2}$ このとき、$dx = \frac{xmax - xmin}{2} dx_{norm}$ となります。積分は以下の形になります。 $$I = \int_{-1}^{1} f\left(\frac{xmax - xmin}{2} x_{norm} + \frac{xmax + xmin}{2}\right) \frac{xmax - xmin}{2} dx_{norm}$$
倍指数関数型変換: 次に、$x_{norm}$ を無限区間 $(-\infty, \infty)$ の変数 $u$ に変換する倍指数関数型変換 $\phi(u)$ を導入します。このプログラムでは以下の変換が用いられています。 $x_{norm} = \phi(u) = \tanh\left(\frac{\pi}{2} \sinh u\right)$ この変換に伴い、$dx_{norm} = \phi'(u) du$ が必要です。 $\phi'(u)$ は次のように計算されます。 $$\phi'(u) = \frac{d}{du} \tanh\left(\frac{\pi}{2} \sinh u\right) = \text{sech}^2\left(\frac{\pi}{2} \sinh u\right) \cdot \frac{\pi}{2} \cosh u = \frac{1}{\cosh^2\left(\frac{\pi}{2} \sinh u\right)} \cdot \frac{\pi}{2} \cosh u$$ これらの変換を適用すると、元の積分は以下の無限区間での積分に変換されます。 $$I = \int_{-\infty}^{\infty} f\left(\frac{xmax - xmin}{2} \phi(u) + \frac{xmax + xmin}{2}\right) \frac{xmax - xmin}{2} \cdot \frac{\pi}{2} \frac{\cosh u}{\cosh^2\left(\frac{\pi}{2} \sinh u\right)} du$$
台形公式による近似: この無限区間での積分を、有限区間 $[umin, umax]$（デフォルトは $[-2.0, 2.0]$）に制限し、台形公式で近似します。ステップ幅を $h = (umax - umin) / nblock$ とすると、積分は次のように近似されます。 $$I \approx h \sum_{i=0}^{nblock} \left[ f\left(\frac{xmax - xmin}{2} \phi(u_i) + \frac{xmax + xmin}{2}\right) \cdot \frac{xmax - xmin}{2} \cdot \frac{\pi}{2} \cdot \frac{\cosh u_i}{\cosh^2\left(\frac{\pi}{2} \sinh u_i\right)} \right]$$ ここで、$u_i = umin + i \cdot h$ です。

コード内の計算は、この数式に対応しています。

phiui = tanh(pi / 2.0 * sinh(ui)) は $x_{norm} = \phi(u_i)$ を計算します。
xi = xmin + (phiui + 1.0) / 2.0 * (xmax - xmin) は、$\phi(u_i)$ から元の区間での $x_i$ を計算します。これは $x = \frac{xmax - xmin}{2} x_{norm} + \frac{xmax + xmin}{2}$ の式で $x_{norm} = \phi(u_i)$ を代入したものに他なりません。
yi = func(xi) は、$f(x_i)$ を計算します。
phi1 = coshnh / (coshsinh * coshsinh) は、$\phi'(u)$ の係数部分 $\frac{\cosh u_i}{\cosh^2\left(\frac{\pi}{2} \sinh u_i\right)}$ を計算します。
ループ内の S += yi * phi1 は、台形公式の和の部分 $\sum f(x_i) \cdot \frac{\cosh u_i}{\cosh^2\left(\frac{\pi}{2} \sinh u_i\right)}$ を累積します。
ループ後の S *= pi / 2.0 * h * (xmax - xmin) / 2.0 は、残りの係数 $\frac{\pi}{2} \cdot h \cdot \frac{xmax - xmin}{2}$ を乗算し、最終的な積分値 $I$ を求めます。

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

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

numpy: 高度な数値計算（特に配列操作）に使用されますが、このスクリプトでは np としてインポートされているものの、直接的なNumPy配列操作は行われていません。将来的な拡張のために含まれている可能性があります。

numpy のインストールは、pip コマンドを使用して行えます。

pip install numpy

math モジュールはPythonの標準ライブラリであるため、別途インストールは不要です。コード内でコメントアウトされている from integ_gauss_legendre import integ_GL は、ガウス・ルジャンドル求積法の実装モジュールであり、このプログラムの実行には直接必要ありません。

必要な入力ファイル

integ_doubleexp.py は、外部の入力ファイルを必要としません。積分対象の関数 func(x) および積分区間 xmin_def, xmax_def、変換区間 umin_def, umax_def、ブロック数 nblock_max_def は、すべてプログラムのソースコード内に直接定義されています。

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

integ_doubleexp.py は、いかなるファイルも生成しません。すべての実行結果は、標準出力（コンソール）にテキスト形式で表示されます。出力には、数値積分結果、解析解との誤差、および解析解自体が含まれます。

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

integ_doubleexp.py はコマンドライン引数を受け付けません。プログラムを実行するには、Pythonインタプリタを使用してスクリプトファイルを直接指定します。

python integ_doubleexp.py

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

以下のコマンドを実行することで、プログラムが開始され、設定された積分区間と関数に対する倍指数関数型積分法の結果がコンソールに表示されます。

python integ_doubleexp.py

実行結果の例:

Numerical integration using by double exponential integra for [-1, 1]

Analytical values:
  f(-1.0)=0.0
  integ(f)[-1.0, 1.0]=1.5707963267948966

nBlock  S               Error
2       1.5459368565150867      -0.02485947027980991
3       1.5714399767512762      0.000643649956379574
4       1.5707963267948966      -2.220446049250313e-16
5       1.5707963267948966      -2.220446049250313e-16
...
32      1.5707963267948966      -2.220446049250313e-16

出力の説明:

Numerical integration using by double exponential integra for [-1, 1]：プログラムの目的を示します。
Analytical values:：解析解の値が表示されます。
- f(-1.0)=0.0：積分区間の下限における関数値。
- integ(f)[-1.0, 1.0]=1.5707963267948966：関数 $f(x)=\sqrt{1-x^2}$ を区間 $[-1, 1]$ で積分した解析解（$\pi/2$）です。
nBlock S Error：以降のテーブルの見出しです。
- nBlock：倍指数関数型積分の際に使用される、変換後の区間におけるブロック数（分割数）を示します。
- S：nBlock の値に対応する数値積分結果です。
- Error：数値積分結果 S と解析解 Fx との差 (S - Fx) を示します。この値が0に近づくほど、数値積分の精度が高いことを意味します。

この出力例から、nBlock が増加するにつれて S の値が解析解に急速に収束し、非常に小さな誤差で積分が実行されていることが確認できます。特に nBlock=4 あたりから、誤差が非常に小さい値になっているのがわかります。