【二分法(Bisection method)】Fortranによる二分法アルゴリズム

はじめに

二分法とは

数値解析における**二分法（Bisection method）**は、解を含む区間の中間点をを求める操作を繰り返すことによって方程式を解く求根アルゴリズムです。その動作が示すように反復法の一種です。

また、$f(x) = 0$ の解を少なくとも一つ求める際、ここでは割愛しますが、ニュートン法と違って、安全確実に解を求めることができます。
しかし、ニュートン法と比較して収束は遅いです。

目標

目標は、中間値の定理が成立する条件のもとで、$f(x) = 0$ となる解を一つ求めることです。

中間値の定理

方程式 $f(x)$が閉区間 $[x_1, x_2]$ で連続であり、$f(x_1)f(x_2) < 0$ であれば、$x_1 < x < x_2$ で $f(x) = 0$ となる $x$ が存在する。

つまり、方程式が連続であり、なおかつ関数値の符号が異なる初期条件を与えることができれば必ず収束します。その際、関数が単調増加あるいは単調減少であれば、区間上限を十分に大きく、区間下限を十分に小さくすることで適切な初期条件となります。

また、**有効桁数は繰り返しの回数に比例して増える（線形収束）**ことから、あらかじめ解の精度を予測することも可能です。

二分法アルゴリズム

1. 前提：$f(x_1)$と$f(x_2)$で符号が異なる区間下限 $x_1$ と区間上限 $x_2$を初期値として定めます。
2. $f(x_1) > 0$ なら $x_p = x_1$ 、 $x_m = x_2$ とおき、そうでなければ、 $x_p = x_2$ 、 $x_m = x_1$ とおきます。
3. 中点 $x_{mid} = \frac{x_p + x_m}{2}$ を計算します。
4. $f(x_{mid})$ の符号を調べます。 $f(x_{mid}) > 0$ なら $x_p = x_{mid}$ とし、$f(x_{mid}) < 0$ なら $x_m = x_{mid}$ とします。
5. 3に戻って操作を繰り返します。

コード

Fortran90の基礎に関してはこちらの記事にまとめています。
【Fortran超入門】プログラム全体の流れ・文法の基礎[FORTRAN]

main.f90

program bisection_method
  implicit none
  real x1, x2, acc, root
  integer eflag
! エラー検出用にeflagに初期値0を設定
  eflag = 0
! 精度(accuracy)値を設定
  acc = 1.0e-7
  x1 = 1.0
  x2 = 2.0
! 上記でx1, x2に初期値を設定しているが、下記read分をコメントアウトすることで対話的に入力可となる
! read *, x1, x2
! 二分法
  call bisection(x1, x2, acc, root, eflag)
! エラーフラグ値が0であれば解（根）を表示
  if (eflag == 0) then
    print *, 'Root=', root
  else
!-- エラーが検出されればエラー値（1or2）を表示
    print *, 'ERROR=', eflag
  endif
end program bisection_method

bisection.f90

subroutine bisection(x1, x2, acc, root, eflag)
  implicit none
  real x1, x2, acc, root
  integer eflag

  integer, parameter::IMAX=30
  real f1, f2, xp, xm, xmid, fmid
  integer i
! y値を計算
  call func(x1, f1)
  call func(x2, f2)
  if (f1*f2 >= 0.0) then
!-- エラー値を１に設定（解がある保証がない）
    eflag = 1
!-- subroutineからメインプログラムへ返る
    return
  endif

  if (f1 > 0.0) then
    xp = x1
    xm = x2
  else
    xp = x2
    xm = x1
  endif

  do i=1, IMAX
    xmid = (xp + xm)/2.0
    print *, xmid, xp, xm, ABS(xp-xm)
!-- 精度が十分であればrootに値を代入してメインプログラムに返る
    if (ABS(xp-xm) < acc*2.0) then
      root = xmid
      return
    endif
    call func(xmid, fmid)
    print *, fmid
    if (fmid > 0.0) xp = xmid
    if (fmid < 0.0) xm = xmid
    if (fmid == 0.0) then
      root = xmid
      return
    endif
  enddo
!-- エラー値を2に設定（精度を満たすことに失敗）
  eflag = 2
!-- 精度を満たせなかったため、中間値を解とする
  root = (xp + xm)/2.0
end subroutine bisection

subroutine func(x, f)
  implicit none
  real x, f
  f = x**2 - 2.0
end subroutine func

コンパイル/実行

ファイルを分割した場合、gfortranコマンドに続けて列挙します。

$ gfortran main.f90 bisection.f90

$ ./a.out
1
3
   2.00000000       3.00000000       1.00000000       2.00000000
   2.00000000
.
.
.
  -1.19209290E-07
   1.41421366       1.41421366       1.41421354       1.19209290E-07
 Root=   1.41421366

参考

二分法 from Wekipedia

おわりに

こちらの記事が役に立ったという方はいいね、お願いします(^^)