Brent 法 ( Brent solver )

Brent 法．

Brent 法のアルゴリズム

Brent 法のアルゴリズム．

Brent 法の特徴

長所

短所

サンプルコード

サンプルコードを以下に示す．

main.f90
 1program main
 2  use brentSolvMod
 3  implicit none
 4  double precision, parameter :: xMin    = 0.d0
 5  double precision, parameter :: xMax    = 4.d0
 6  integer         , parameter :: iterMax = 100
 7  double precision, parameter :: crit    = 1.d-12
 8  double precision            :: x1, x2
 9  
10  ! ------------------------------------------------------ !
11  ! --- [1] solve by bisection Method                  --- !
12  ! ------------------------------------------------------ !
13  x1 = 0.5d0
14  x2 = 1.2d0
15  call brentSolver( x1, x2, iterMax, crit )
16  write(6,*) "[Main@brentSolver] x1, x2, res, sin(x)"
17  write(6,*) x1, x2, analyticFunc( x1 )
18  write(6,*) "[Main@brentSolver] Answer should be..."
19  write(6,*) sin(x1), 0.707d0
20  
21  return
22end program main
brentSolvMod.f90
  1module brentSolvMod
  2contains
  3
  4  ! ====================================================== !
  5  ! === analyticFunc                                   === !
  6  ! ====================================================== !
  7  function analyticFunc( xv )
  8    implicit none
  9    double precision, intent(in)  :: xv
 10    double precision              :: analyticFunc
 11
 12    analyticFunc = sin( xv ) - 0.707d0
 13    return
 14  end function analyticFunc
 15
 16
 17  ! ====================================================== !
 18  ! === brent Solver for Root Finding                  === !
 19  ! ====================================================== !
 20  subroutine brentSolver( x1, x2, iterMax, crit )
 21    implicit none
 22    integer         , intent(in)    :: iterMax
 23    double precision, intent(in)    :: crit
 24    double precision, intent(inout) :: x1, x2
 25    integer                         :: iter
 26    double precision                :: xa, xb, xc, xd, xe, xm
 27    double precision                :: fa, fb, fc
 28    double precision                :: sval, pval, qval, rval
 29    double precision                :: epsilon, tolerance, judge1, judge2
 30
 31    ! ------------------------------------------------------ !
 32    ! --- [1] Preparation before Main Loop               --- !
 33    ! ------------------------------------------------------ !
 34    epsilon = crit
 35    xa      = x1
 36    xb      = x2
 37    xc      = xb
 38    xe      = 0.d0
 39    fa      = analyticFunc( xa )
 40    fb      = analyticFunc( xb )
 41    fc      = fb
 42    if ( fa*fb.gt.0.d0 ) then
 43       write(6,*) '[brentSolver] ERROR!! x1 and x2 must be bracketed... '
 44       write(6,*) '[brentSolver]   i.e. ( f(x1)<0<f(x2) ) or ( f(x2)<0<f(x1) ) '
 45       stop
 46    endif
 47    ! ------------------------------------------------------ !
 48    ! --- [2] Main Loop                                  --- !
 49    ! ------------------------------------------------------ !
 50    do iter=1, iterMax
 51       if ( fb*fc.gt.0.d0 ) then
 52          xc = xa
 53          fc = fa
 54          xd = xb-xa
 55          xe = xd
 56       endif
 57       if ( abs(fc).lt.abs(fb) ) then
 58          xa = xb
 59          xb = xc
 60          xc = xa
 61          fa = fb
 62          fb = fc
 63          fc = fa
 64       endif
 65       ! -- convergence check   -- !
 66       tolerance = 2.d0*epsilon*abs( xb ) + 0.5d0*crit
 67       xm        = 0.5d0 * ( xc-xb )
 68       if ( ( abs(xm).le.tolerance ).or.( fb.eq.0.d0 ) ) exit
 69       ! -- solution candidates -- !
 70       if ( ( abs(xe).ge.tolerance ).and.( abs(fa).gt.abs(fb) ) ) then
 71          ! -- Attempts inverse uadratic interpolation -- !
 72          sval = fb / fa
 73          if ( xa.eq.xc ) then
 74             pval = 2.d0 * xm * sval
 75             qval = 1.d0 - sval
 76          else
 77             qval = fa / fc
 78             rval = fb / fc
 79             pval = sval * ( 2.d0*xm*qval*( qval-rval ) - ( xb-xa )*( rval-1.d0 ) )
 80             qval = ( qval-1.d0 )*( rval-1.d0 )*( sval-1.d0 )
 81          endif
 82          ! -- check wheter in bounds -- !
 83          if ( pval.gt.0.d0 ) qval = -qval
 84          pval   = abs(pval)
 85          judge1 = 2.d0*pval
 86          judge2 = min( 3.d0*xm*qval-abs( tolerance*qval ), abs( xe*qval ) )
 87          if ( judge1.lt.judge2 ) then
 88             ! -- Accept interpolation -- !
 89             xe = xd
 90             xd = pval / qval
 91          else
 92             ! -- Decline, use bisection -- !
 93             xd = xm
 94             xe = xd
 95          endif
 96       else
 97          ! -- Bounds decreasing too slowly, use bisection method -- !
 98          xd = xm
 99          xe = xd
100       endif
101       ! -- move last best guess to a -- !
102       xa = xb
103       fa = fb
104       if ( abs(xd).gt.tolerance ) then
105          xb = xb + xd
106       else
107          xb = xb + sign( tolerance, xm )
108       endif
109       fb = analyticFunc( xb )
110    enddo
111    ! ------------------------------------------------------ !
112    ! --- [3] return x1 and x2                           --- !
113    ! ------------------------------------------------------ !
114    x1 = xb
115    x2 = fb
116    
117    return
118  end subroutine brentSolver
119
120end module brentSolvMod

実行結果例

実行結果例を以下に示す．

kent@euler ~/fortran/brentSolver $ ./main
[Main@brentSolver] x1, x2, res, sin(x)
0.78524716339422451       -6.5658589676331758E-013  -6.5658589676331758E-013
[Main@brentSolver] Answer should be...
0.70699999999934338       0.70699999999999996