今回はfortranを使って、積分の収束値を求めていきたいと思います。
使うのは「if」文です。
条件分岐をさせて、収束する値をコンピューターに自動で決めてもらおうと思います。
具体的にはこのようなコードになります。
program trapezoidal
implicit none
external f
integer :: n,i,j
real(8) :: a,b,h, x(0:1000000),sint,s(0:100),f,delta
a=-1.0d0; b=1.0d0;
s(0)=0
n=10
do j=1,100
n = n*2
h=(b-a)/(n+1)
sint= h*(f(a) + f(b))/2
do i=0,n-1
x(i)=a+i*h
sint=sint + h*f(x(i))
end do
s(j) = sint
if (s(j-1)==0) then
write(*,*)"1 loop"
else
write(*,*) ' integral over ', a,'~',b
write(*,*) ' S = ', s(j)
write(*,*) j,"loops"
delta = 10.0*2.0**(-j+1)*(1-s(j)/s(j-1))
if(delta < 1.0d-6) exit
end if
end do
stop
end program trapezoidal
function f(x)
implicit none
real(8) :: f,x
f=x*sin(x)
return
end function f
上の例では
$$
f(x) = x\sin x
$$
を-1~1の範囲で積分しています。
では、具体的に説明していきます。
こちらのコードでは 「【fortran】台形公式を使ってガウス積分を求めてみた！！」 で使った台形公式を使っているのでぜひ参考にしてください。
if文で収束値を自動で決める
doループ
まず、10~33行目のdoループ
s(0)=0
n=10
do j=1,100
n = n*2
h=(b-a)/(n+1)
sint= h*(f(a) + f(b))/2
do i=0,n-1
x(i)=a+i*h
sint=sint + h*f(x(i))
end do
s(j) = sint
if (s(j-1)==0) then
write(*,*)"1 loop"
else
write(*,*) ' integral over ', a,'~',b
write(*,*) ' S = ', s(j)
write(*,*) j,"loops"
delta = 10.0*2.0**(-j+1)*(1-s(j)/s(j-1))
if(delta < 1.0d-6) exit
end if
end do
が気になると思います。
これは、
- 初期値nに2をかける
- そのn（台形の細さ）で\(f(x) = x \sin x\)の積分を行う
- nでやったときの積分値とn-1でやったときの積分値 （一つ前のループで出た値） の比を取る
- これらの比が\(10^{-6}\)以下（収束しているとみなす）だったらその値を収束値として出す
- \(10^{-6}\)以上なら1に戻る
ということをしています。
4.5でif文を使っています。
ではこのループ処理の中身を解説していきたいと思います。
1.初期値nに2をかける
10,11行目は初期値\(n=10,s(0) =0\)を決めています。
\(s(0)\)は積分値の初期値でとりあえず0としています。
また、\(n=10\)に2を次々とかけていきます。
2.そのn（台形の細さ）で\(f(x) = x \sin x\)の積分を行う
h=(b-a)/(n+1) sint= h*(f(a) + f(b))/2 do i=0,n-1 x(i)=a+i*h sint=sint + h*f(x(i)) end do s(j) = sint
「【fortran】台形公式を使ってガウス積分を求めてみた！！」 でやったように積分を行っていきます。
\(j\)回目のループで出た積分値と\(j+1\)回目に出た積分値を比較したいので、\(s(j)\)に積分値を代入しています。
3．nでやったときの積分値とn-1でやったときの積分値 （一つ前のループで出た値） の比を取る
if (s(j-1)==0) then
write(*,*)"1 loop"
else
write(*,*) ' integral over ', a,'~',b
write(*,*) ' S = ', s(j)
write(*,*) j,"loops"
delta = 10.0*2.0**(-j+1)*(1-s(j)/s(j-1))
if(delta < 1.0d-6) exit
end if
積分値の比較には
$$
\Delta_n =\Delta\times 2^{-n}\left|1-\frac{I(\Delta_{n+1})}{I(\Delta_n)}\right|\tag{1}
$$
を使います。
積分の幅を\(\Delta\)とした場合の積分値を\(I(\Delta)\)とします。このとき、初期値\(\Delta\)を\(2^{n}\)倍したときの\(\Delta_n\)が十分小さい（\(10^{-6}\)くらい）なら、収束しているといえます。
(1)式に
$$
\begin{cases}
\Delta = 10(初期値)\\
I(\Delta_{n+1})=s(j)\\
I(\Delta_{n})=s(j-1)
\end{cases}
$$
を代入したものが
delta = 10.0*2.0**(-j+1)*(1-s(j)/s(j-1))
となります。
この「delta」が\(10^{-6}\)より小さいときループを抜けます。
if(delta < 1.0d-6) exit
4.比が\(10^{-6}\)以下だったらその値を収束値として出す
ループから外れ、最後に残った値が収束値となります。
\(f(x) =x\sin x\)
のときは
S = 0.60233834098286954
となりました。
