【fortran】if文を使って積分の収束値を自動で決める!!

プログラミング

今回は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

が気になると思います。

これは、

  1. 初期値nに2をかける
  2. そのn(台形の細さ)で\(f(x) = x \sin x\)の積分を行う
  3. nでやったときの積分値とn-1でやったときの積分値 (一つ前のループで出た値) の比を取る
  4. これらの比が\(10^{-6}\)以下(収束しているとみなす)だったらその値を収束値として出す
  5. \(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 

となりました。

タイトルとURLをコピーしました