【fortran】台形公式を使ってガウス積分を求めてみた！！

今回はfortranで台形公式を使ってガウス積分を求めてみたいと思います。

台形公式は関数の作る面積（積分したもの）を近似して求めたものです。

台形公式でガウス積分を求める

台形公式でガウス積分を求める

コード

では始めにコードを書いておきます。

program sekibun
implicit none
external f
integer :: n,i
real(8) :: a,b,h,x(0:10000),s,f

a=-5.0; b=5.0;
write (*,*) 'input number'
read (*,*) n
h = (b-a)/(n+1)

s = h*(f(a)+f(b))/2;

do i=0,n-1
x(i) = a + h*i
s = s + h*f(x(i))
end do

write(*,*)'S =',s
write(*,*)'S**2=',s**2

stop
end program sekibun

!以下、積分する関数
function f(x)
implicit none
real(8) :: f,x
 f=exp(-x**2)
return
end function f

台形公式

コードの説明をするうえで欠かせないのが台形公式です。

$a\to b$の範囲で上のような関数を考えます。

この面積を求めるときに、図のような細切れの台形の面積の和を考えます。

台形の細さが細いほど、実際の値に近いものが出てくることが分かりますね。

これが台形公式です。数式で書くと

$$
\begin{eqnarray}
\int_{a}^{b} f(x)dx &=& \sum_{i=0}^{N}\int_{x_i}^{x_i+1}f(x)dx\\
&=& \frac{h}{2}\left[(f_a+f_{a+x_0})+ (f_{a+x_0}+f_{a+x_1})+\cdots +(f_{a +x_N}+f_{b})+ \right]\\
&=&\frac{h}{2}(f_{a}+f_{b})+h(f_{a+x_0}+\cdots + f_{a+x_N})\tag{1}
\end{eqnarray}
$$

となります。

コード解説

台形の細さ

まず、7~10行目のコード

 a=-5.0; b=5.0;
write (*,*) 'input number'
read (*,*) n
h = (b-a)/(n+1)

では

$a,b$の値を決める
$a~b$間を何分割するか（台形の細さ）を読み込む

といったことをしています。

特に台形の細さは「read」を使ってキーボードから決めるようにしています。

一項目

$a,b$が決まったので、(1)式の一項目

$$
\frac{h}{2}(f_{a}+f_{b})
$$

を12行目で計算しています。

s = h*(f(a)+f(b))/2;

二項目

(1)式の二項目は

do i=0,n-1
x(i) = a + h*i
s = s + h*f(x(i))
end do

とします。

ここでは配列$x(i)$を使って座標を求め、その座標を面積を求めたい関数に入れています。$f(x(i))$

それらに台形の高さ$h$をかけ、すべて足しています。

$$
h(f_{a+x_0}+\cdots + f_{a+x_N})
$$

積分する関数

26~31行目に積分する関数を書きました。

function f(x)
implicit none
real(8) :: f,x
 f=exp(-x**2)
return
end function f

関数はこのように別々書いておくと、コードが汚くならずに済みます。

今回は関数として

$$
f(x)=e^{-x^2}
$$

を使いました。これを$-\infty \to \infty$の範囲で積分すれば、答えが

$$
\sqrt{\pi}
$$

になります。

今回は範囲が$-5\to5$ですが、$e^{-5^2} \sim 10^{-11}$ほどなので、十分と考えられます。

実行してみる

このコードを$n=100$で実行すると

S =   1.7724538509002776
S**2=   3.1415926535712235

となり、10桁ほどの精度で答えを出すことができています！！