区間が無限の積分を数値計算で行う(ガウス積分)

積分を台形公式などを用いて数値計算するには、積分区間を細かく分割することが重要となるが、積分範囲が広くなるにつれて計算量が増大する。特に積分区間が無限である場合(広義積分)はそのまま数値計算を行うことは困難である。

例を挙げると、$\int_{-\infty}^{\infty} e^{-x^2}dx=\sqrt{\pi}$(ガウス積分)は自然科学で頻出の積分だが、そのままの形では台形公式を用いて計算することはできない。
これを計算するには以下のような工夫で無限区間を有限にする必要がある。

###1.積分区間を途中で打ち切る
これは指数関数のように急激に減衰する場合に有効である。ガウス積分もこの方法を用いると[-3:3]の区間で数値積分するだけで解析解との誤差は$10^{-5}$のオーダーに収まる。しかし$\int_{-\infty}^{\infty}\frac{dx}{1+x^2}$では被積分関数の減衰が緩やかであるためにこの方法で計算することは難しい。

###2.置換積分を用いる
定義域が無限区間でも値域が有限となるような関数を用いて置換積分を行えば積分区間を狭くすることができる。例えば$y=\tanh(x)$や$y=\frac{2\arctan(x)}{\pi}$の値域は[-1:1]となり積分区間の大幅な短縮ができる。

ガウス積分をこの方法を用いて数値計算してみよう。

\int_{-\infty}^{\infty} e^{-x^2}dxについてy=\frac{2\arctan x}{\pi}とする。\\
x:-\infty\rightarrow\infty \Leftrightarrow y:-1\rightarrow 1\\
x=\tan(\frac{\pi y}{2})\\
dx=\frac{\pi}{2\cos^2(\frac{\pi y}{2})}dyより\\
\int_{-\infty}^{\infty} e^{-x^2}dx=\int_{-1}^{1}\frac{\pi e^{-\tan^2({\frac{\pi y}{2}})}}{2\cos^2(\frac{\pi y}{2})}dy

被積分関数をプロットすると図のようになり、[-1:1]で無限に発散する点はないので、これに対して台形公式を適用すればよい精度でガウス積分が計算できると期待される。

c言語での台形公式のコードは以下の通り

#include <stdio.h>
#include <math.h>

double f(double x){
double y;

y = M_PI*exp(-pow(tan(M_PI*x/2),2))/(2*pow(cos(M_PI*x/2),2));

return y;
}

int main(void){
int i,N=50000;
double S,dx;

dx=1./(double)N;
S=0.;

for(i=0;i<=N-1;i++){
S += (f(i*dx)+f((i+1)*dx))/2.*dx;
}

printf("計算結果: %e 解析解: %e",2.*S,sqrt(M_PI));

return 0;
}

計算結果は以下の通り

計算結果: 1.772454e+000 解析解: 1.772454e+000