今回は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
となりました。