【fortran】二分法

　今回は非線形方程式

$$e^x = x^2$$

の解を数値計算で解く。今回は二分法で求める。

二分法
サンプルプログラム

二分法

　ある関数$y = f(x)$を考える。区間$[a,b]$で連続な$f(x)$が

$$f(a)f(b) < 0 \tag{1}$$

ならば、この区間に$f(x) = 0$を満たす解が少なくとも1つ存在する。今、始めに式(1)を満たす$a,b$が与えられたとき、解の推定値として

$$c = \frac{a+b}{2} \tag{2}$$

を与える。もし$f(c) = 0$であれば、$c$は$f(x)$の解である。$f(c)\neq 0 $の場合、$f(c)$と同符号の関数値を持つ$a , b$と$c$を置き換える。つまりまとめると、

$f(a)f(c) > 0$ ($f(b)f(c) < 0$)の場合、$a$を$c$で置き換える

$f(a)f(c) < 0$ ($f(b)f(c) > 0$)の場合、$b$を$c$で置き換える

この操作を繰り返し行う。繰り返し行い、$|a-b|$が十分に小さい値になると$c$は解に収束したとみなす。

　二分法はは区間$[a,b]$に解があれば必ず真の解に収束するロバストな解法であるが、収束が遅いという欠点を持つ。

サンプルプログラム

program bisec

  implicit none                            !暗黙の型宣言
  real, parameter    :: eps = 1.0e-6       !適当な小さい値
  integer, parameter :: loopmax = 100      !反復回数の上限（これを設定しないと無限に計算する可能性、、、）
  real               :: a,b,c,res
  real               :: f
  integer            :: loop

  write(6,*) 'input a and b'
  read(5,*) a,b
  if(f(a)*f(b) &gt; 0.0) then
    write(6,*) 'f(a)*f(b) &gt; 0'
    stop
  end if

  do loop = 1 , loopmax
    c = 0.5*(a+b)
    if(f(a)*f(c) &lt; 0.0) then
      b = c
    else
      a = c
    end if
    res = abs(a-b)
    if(res &lt; eps) exit
  end do

  if(res &lt; eps) then
    write(6,*) 'converge: loop =',loop
    write(6,*) 'x =',0.5*(a+b)
  else
    write(6,*) 'not converge!'
  end if

end program bisec

real function f(x)

  implicit none
  real :: x

  f = x**2-exp(x)                         !関数（exp(x)とx^2の交点を求める）

end function f

program bisec

implicit none !暗黙の型宣言

real, parameter :: eps = 1.0e-6 !適当な小さい値

integer, parameter :: loopmax = 100 !反復回数の上限（これを設定しないと無限に計算する可能性、、、）

real :: a,b,c,res

real :: f

integer :: loop

write(6,*) 'input a and b'

read(5,*) a,b

if(f(a)*f(b) > 0.0) then

write(6,*) 'f(a)*f(b) > 0'

stop

end if

do loop = 1 , loopmax

c = 0.5*(a+b)

if(f(a)*f(c) < 0.0) then

b = c

else

a = c

end if

res = abs(a-b)

if(res < eps) exit

end do

if(res < eps) then

write(6,*) 'converge: loop =',loop

write(6,*) 'x =',0.5*(a+b)

else

write(6,*) 'not converge!'

end if

end program bisec

real function f(x)

implicit none

real :: x

f = x**2-exp(x) !関数（exp(x)とx^2の交点を求める）

end function f

使い方はターミナルで

gfortran bisec.f90

1	gfortran bisec.f90

とし、実行して$a,b$を代入する。

　例えば$a=1,b=-1$と入力すると、

converge: loop =          21
x = -0.703467846

1 2	converge: loop = 21 x = -0.703467846

と出力され、解が近似的に出力される。