muller Subroutine

private subroutine muller(me, ax, bx, fax, fbx, xzero, fzero, iflag)

Improved Muller method (for real roots only). Will fall back to bisection if any step fails.

Reference

  • D. E. Muller, "A Method for Solving Algebraic Equations Using an Automatic Computer", Mathematical Tables and Other Aids to Computation, 10 (1956), 208-215. link
  • Roots.jl, Julia version of standard Muller

Type Bound

muller_solver

Arguments

Type IntentOptional Attributes Name
class(muller_solver), intent(inout) :: me
real(kind=wp), intent(in) :: ax

left endpoint of initial interval
real(kind=wp), intent(in) :: bx

right endpoint of initial interval
real(kind=wp), intent(in) :: fax

f(ax)
real(kind=wp), intent(in) :: fbx

f(ax)
real(kind=wp), intent(out) :: xzero

abscissa approximating a zero of f in the interval ax,bx
real(kind=wp), intent(out) :: fzero

value of f at the root (f(xzero))
integer, intent(out) :: iflag

status flag (0=root found, -2=max iterations reached)

Calls

proc~~muller~~CallsGraph proc~muller root_module::muller_solver%muller f f proc~muller->f proc~bisect root_module::bisect proc~muller->proc~bisect proc~choose_best root_module::choose_best proc~muller->proc~choose_best proc~converged root_module::root_solver%converged proc~muller->proc~converged proc~solution root_module::root_solver%solution proc~muller->proc~solution

Source Code

    subroutine muller (me,ax,bx,fax,fbx,xzero,fzero,iflag)

    implicit none

    class(muller_solver),intent(inout) :: me
    real(wp),intent(in)  :: ax      !! left endpoint of initial interval
    real(wp),intent(in)  :: bx      !! right endpoint of initial interval
    real(wp),intent(in)  :: fax     !! `f(ax)`
    real(wp),intent(in)  :: fbx     !! `f(ax)`
    real(wp),intent(out) :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out) :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)  :: iflag   !! status flag (`0`=root found, `-2`=max iterations reached)

    real(wp) :: a,b,c,cx,fa,fb,fc,fcx,x,q,q2,q1,aa,bb,cc,delta,dp,dm,denon,bprev,fbprev
    integer  :: i    !! iteration counter
    logical  :: x_ok !! the new estimate `x` is ok to use

    iflag = 0

    ! pick a third point in the middle [this could also be an optional input]
    cx  = bisect(ax,bx)
    fcx = me%f(cx)
    if (me%solution(cx,fcx,xzero,fzero)) return

    ! [a,b,c]
    a = ax; fa = fax
    b = cx; fb = fcx
    c = bx; fc = fbx

    bprev = huge(1.0_wp)
    fbprev = huge(1.0_wp)

    do i = 1, me%maxiter

        ! muller step:
        q     = (c - b)/(b - a)
        q2    = q**2
        q1    = q + 1.0_wp
        aa    = q*fc - q*q1*fb + q2*fa
        bb    = (q1+q)*fc - q1**2*fb + q2*fa
        cc    = q1*fc
        delta = sqrt(max(0.0_wp, bb**2 - 4.0_wp * aa*cc)) ! to avoid complex roots
        dp    = bb + delta
        dm    = bb - delta
        if (abs(dp) > abs(dm)) then
            denon = dp
        else
            denon = dm
        end if
        x_ok = denon /= 0.0_wp
        if (x_ok) x = c - 2.0_wp*(c - b)*cc/denon

        ! make sure that x is ok, in the correct interval, and distinct.
        ! if not, fall back to bisection on that interval
        if (fa*fb < 0.0_wp) then  ! root in (a,b)
            if (.not. x_ok .or. x<=a .or. x>=b) x = bisect(a,b)
            c  = b
            fc = fb
            b  = x
        else  ! root in (b,c)
            if (.not. x_ok .or. x<=b .or. x>=c) x = bisect(b,c)
            a  = b
            fa = fb
            b  = x
        end if
        ! values are now [a,b,c], with b being the new estimate

        ! function evaluation for next estimate:
        fb = me%f(b)
        if (abs(fb)<=me%ftol) exit

        ! stopping criterion
        if (me%converged(a,c) .or. i == me%maxiter) then
            if ( i == me%maxiter ) iflag = -2 ! max iterations exceeded
            exit
        end if

        bprev = b
        fbprev = fb

    end do

    call choose_best(b,bprev,fb,fbprev,xzero,fzero)

    end subroutine muller