    Next: 7 常微分方程式の数値解法：初期値問題 Up: 6 方程式の数値解法 Previous: 2 二分法

# 3 Newton-Raphson法と二分法のハイブリッド法

Newton-Raphson法では、推定値 の真の解 からのずれを とすると、(6.2)式より (23)

だが、 (24) (25)

だから、 (26)

となり、良い推定値さえ得られていれば、収束が速いことが判る。しかし、 解を捜す区間に変曲点があると、近傍の解には収束しないか、収束しても 一般には収束は速くない。 そこで、二分法と組み合わせて、確実に、しかも収束の速いプログラムを用意すると 良い。以下のプログラムは、そんな例である function rtsafe(funcd,x1,x2,xacc)
implicit NONE
Integer MAXIT
Real rtsafe,x1,x2,xacc
external funcd
Parameter (MAXIT=100)   !Maximum allowed number of iterations.
!	Using a combination of Newton-Raphson and bisection, find the root of a
!	function bracketed between x1 and x2. The root, returned as the
!	function value 'rtsafe', will be refined until its accuracy is known
!	within +-'xacc'. 'funcd' is a user-supplied subroutine which returns
!	both the function value and the first derivative of the function.
Integer j
Real df,dx,dxold,f,fh,fl,temp,xh,xl
call funcd(x1,fl,df)
call funcd(x2,fh,df)
if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
*	pause 'root must be bracketed in rtsafe'
if(fl.eq.0.)then
rtsafe=x1
return
else if(fh.eq.0.)then
rtsafe=x2
return
else if(fl.lt.0.)then   !Orient the search so that f(x1)<0.
xl=x1
xh=x2
else
xh=x1
xl=x2
endif
rtsafe=.5*(x1+x2)       !Initialize the guess for root,
dxold=abs(x2-x1)        !the "stepsize before last,"
dx=dxold                !and the last step.
call funcd(rtsafe,f,df)
do j=1,MAXIT            !Loop over allowed iterations.
if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0. !Bisect if Newton out of
*     .or. abs(2.*f).gt.abs(dxold*df) ) then      !range, or not decreasing
dxold=dx                                     !fast enough.
dx=0.5*(xh-xl)
rtsafe=xl+dx
if(xl.eq.rtsafe)return        !Change in root is negligible.
else                            !Newton step acceptable. Take it.
dxold=dx
dx=f/df
temp=rtsafe
rtsafe=rtsafe-dx
if(temp.eq.rtsafe)return
endif
if(abs(dx).lt.xacc) return      !Convergence criterion.
call funcd(rtsafe,f,df)         !The one new function evaluation per
if(f.lt.0.) then                !iteration.
xl=rtsafe                     !Maintain the bracket on the root.
else
xh=rtsafe
endif
end do
pause 'rtsafe exceeding maximum iterations'
return
END


Jun Makino