はじめに
広義積分であるガウス積分を計算する場合は、一般的に重積分を用いる。
なので、計算方法は少し複雑である。
そこで、今回は別の方法でガウス積分の近似値を数値計算を行い求められるか調査した。
具体的には、マクローリン展開を用いることでガウス積分を$x^k$の積分の無限和として表した。次に、計算機で近似値を計算することができるか調査した。
以下、結果の概要を述べる。
ガウス関数の近似
マクローリン展開を用いてガウス関数の近似を行った。
ガウス関数の計算誤差
ガウス関数の計算誤差は、$x$の値が大きくなるほど急激に増加する傾向がある。
したがって、積分範囲をうまく調整しないと値が発散してしまう。
ガウス積分の数値計算
計算機を用いてガウス積分の数値計算を行った。
このように、積分範囲をうまく設定できれば、かなり高い精度でガウス積分の数値計算を行うことができる。
また、以下にstreamlitを用いたpythonのwebアプリも作成した。
ガウス積分について
今回扱うガウス積分は以下の積分である。
S=\int_{0}^{\infty}e^{-x^2}dx=\frac{\sqrt{\pi}}{2}
このガウス積分の近似値を計算機による数値計算で求めたい。
マクローリン展開
関数$e^{x}$について、収束半径は無限大で、以下のようにマクローリン展開を行うことができる。
e^x=\sum_{n=0}^{\infty}\frac{1}{n!}x^{n}
したがって、関数$e^{-x^2}$のマクローリン展開は以下のようになる。
e^{-x^2}=\sum_{n=0}^{\infty}\frac{1}{n!}(-x^{2})^n
この関数は一様収束することを認めると、以下の式が成立する。
(シグマとインテグラルの交換則)
ガウス積分と極限
S=\int_{0}^{\infty} e^{-x^2}dx=\int_{0}^{\infty}\sum_{n=0}^{\infty}\frac{1}{n!}(-x^2)^{n}dx=\sum_{n=0}^{\infty}\frac{1}{n!}\int_{0}^{\infty}(-x^2)^{n}dx=\sum_{n=0}^{\infty}\frac{1}{n!}[(-1)^{n} \frac{1}{2n+1}x^{2n+1}]^{\infty}_0
そこで、$b\to \infty,n\to \infty$のとき、以下の式が成立する。
S=\lim_{b\to \infty}\lim_{n\to \infty}\sum_{n=0}^{n}(-1)^{n} \frac{1}{2n+1}\frac{1}{n!}[x^{2n+1}]^{b}_0=\lim_{b\to \infty}\lim_{n\to \infty}\sum_{n=0}^{n}(-1)^{n} \frac{(b^{2n+1})}{(2n+1)n!}
数値計算
Sの近似値を算出するため、以下の関数を計算機で計算する。
g(n,b)=\sum_{n=0}^{n}(-1)^{n} \frac{(b^{2n+1})}{(2n+1)n!}
ただし、$n=100,0\le b\le 4$のもとで計算した。
プログラム及び結果
ガウス関数について
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# xの分割数
m=100
x=np.linspace(0,3,m)
# テイラー展開の項数
n=100
S_ary=[]
for p in range(m):
S=0
for i in range(n+1):
S1=1
for k in range(i):
f=(-x[p]**2)/(k+1)
S1=S1*f
S+=S1
S_ary.append(S)
# ガウス関数の理論値
gauss=np.exp(-x**2)
# 誤差の計算
err=abs((np.array(S_ary)-gauss)/gauss)
plt.plot(x,err,color='red')
plt.xlabel("xの値")
plt.ylabel("ガウス関数の近似値の誤差")
plt.title("マクローリン展開を用いたガウス関数の計算誤差")
plt.grid()
plt.legend()
plt.savefig('マクローリン展開を用いたガウス関数の計算誤差.png')
plt.show()
plt.plot(x,S_ary,label='ガウス関数の近似値',color='red')
plt.plot(x,gauss,label='ガウス関数の理論値',color='blue')
plt.xlabel("xの値")
plt.ylabel("ガウス関数の値")
plt.title("マクローリン展開を用いたガウス関数の計算")
plt.grid()
plt.legend()
plt.savefig('マクローリン展開を用いたガウス関数の計算.png')
plt.show
ガウス積分について
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# テイラー展開の各項の積分値を計算する関数
def f(n,b):
B=1
A=(-1)**n*b/(2*n+1)
for k in range(1,n+1):
B=B*(b**2)/k
f=A*B
return f
# テイラー展開の各項の積分値を合計する関数
def g(n,b):
S=0
for i in range(n+1):
S+=f(i,b)
return S
# xの分割数
m=100
b_ary=np.linspace(0,4,m)
# テイラー展開の項数
n=100
S_ary=[]
for p in range(m):
b=b_ary[p]
S=g(n,b)
S_ary.append(S)
plt.plot(b_ary,[(math.pi)**0.5/2]*m,'--',label='ガウス積分の理論値の半分')
plt.plot(b_ary,S_ary,label='g(n,b)の値')
plt.xlabel("b(積分の上限値)")
plt.ylabel("g(n,b)の値")
plt.title("マクローリン展開を用いたガウス積分の数値計算")
plt.legend()
plt.grid()
plt.savefig('マクローリン展開を用いたガウス積分の数値計算.png')
plt.show()
ただし、階乗計算を行った後に除算すると値がオーバーフローする。
ゆえに、各項分数計算を行うことで対処した。
まとめ
今回は、ガウス積分の近似値を計算機を用いた数値計算で算出することに挑戦した。
具体的には、まずガウス関数をマクローリン展開をした。
次に、各項について積分値を計算した。
最後に各項を合計した。
ガウス関数及びガウス積分は、高校数学では確率統計や微積分の分野に頻出の事項である。
なので、これを機に苦手意識のある方は復習をすることを勧める。
参考文献
ガウス積分について
積分と極限（無限和）の交換について
一様収束について