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

Newton-Raphson法では、推定値$x_{n}$の真の解$x$からのずれを $\varepsilon_{n}$とすると、(6.2)式より

$$\varepsilon_{n+1}=\varepsilon_n-{{f(x_n)}\over{f'(x_n)}}$$ (23)

だが、

$$f(x+\varepsilon)=f(x)+\varepsilon f'(x)+\varepsilon^2{{f''(x)}\over{2}}+\cdots,$$ (24)

$$f'(x+\varepsilon)=f'(x)+\varepsilon f''(x)+\cdots$$ (25)

だから、

$$\varepsilon_{n+1}=-\varepsilon_n^2{{f''(x_n)}\over{2f'(x_n)}}$$ (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