zhang Subroutine

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

Zhang's method (with corrections from Stage).

Reference

  • A. Zhang, "An Improvement to the Brent's Method", International Journal of Experimental Algorithms (IJEA), Volume (2) : Issue (1) : 2011. https://www.cscjournals.org/download/issuearchive/IJEA/Volume2/IJEA_V2_I1.pdf
  • S. A. Stage, "Comments on An Improvement to the Brent's Method", International Journal of Experimental Algorithms (IJEA), Volume (4) : Issue (1) : 2013. https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.740.923&rep=rep1&type=pdf

Type Bound

zhang_solver

Arguments

Type IntentOptional Attributes Name
class(zhang_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~~zhang~~CallsGraph proc~zhang root_module::zhang_solver%zhang f f proc~zhang->f proc~bisect root_module::bisect proc~zhang->proc~bisect proc~choose_best root_module::choose_best proc~zhang->proc~choose_best proc~converged root_module::root_solver%converged proc~zhang->proc~converged proc~secant root_module::secant proc~zhang->proc~secant proc~solution root_module::root_solver%solution proc~zhang->proc~solution proc~swap root_module::swap proc~zhang->proc~swap proc~secant->proc~bisect

Source Code

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

    implicit none

    class(zhang_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,fa,fb,fc,s,fs
    integer :: i !! iteration counter

    iflag = 0
    a  = ax
    b  = bx
    fa = fax
    fb = fbx

    do i = 1, me%maxiter

        c = bisect(a,b)
        fc = me%f(c)
        if (me%solution(c,fc,xzero,fzero)) return

        if (fa/=fc .and. fb/=fc) then
            ! inverse quadratic interpolation
            s = a*fb*fc/((fa-fb)*(fa-fc)) + &
                b*fa*fc/((fb-fa)*(fb-fc)) + &
                c*fa*fb/((fc-fa)*(fc-fb))
            if (a<s .and. s<b) then
                fs = me%f(s)
                if (abs(fs)<=me%ftol) then
                    xzero = s
                    fzero = fs
                    return
                end if
            else
                ! s is not in (a,b)
                s = c ! just use this (there are 3 options in the reference)
                fs = fc
            end if
        else
            ! secant
            if (fa*fc<0.0_wp) then ! root in [a,c]
                s = secant(a,c,fa,fc,ax,bx)
            else ! root in [c,b]
                s = secant(c,b,fc,fb,ax,bx)
            end if
            fs = me%f(s)
            if (me%solution(s,fs,xzero,fzero)) return
        end if

        if (c>s) then
            ! ensures a <= c <= s <= b
            call swap(s,c)
            call swap(fs,fc)
        end if

        if (fc*fs<0.0_wp) then       ! root on [c,s]
            a = c
            b = s
            fa = fc
            fb = fs
        else if (fa*fc<0.0_wp) then  ! root on [a,c]
            b = c
            fb = fc
        else                         ! root on [s,b]
            a = s
            fa = fs
        end if

        if (me%converged(a,b)) exit
        if (i == me%maxiter) iflag = -2 ! max iterations reached

    end do

    ! pick the one closest to the root:
    call choose_best(a,b,fa,fb,xzero,fzero)

    end subroutine zhang