ridders Subroutine

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

Ridders method to find a root of f(x).

See also

  • Ridders, C., "A new algorithm for computing a single root of a real continuous function", IEEE Trans. on Circuits and Systems, Vol 26, Issue 11, Nov 1979.

Type Bound

ridders_solver

Arguments

Type IntentOptional Attributes Name
class(ridders_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, -3=singularity in the algorithm)

Calls

proc~~ridders~~CallsGraph proc~ridders root_module::ridders_solver%ridders f f proc~ridders->f proc~bisect root_module::bisect proc~ridders->proc~bisect proc~converged root_module::root_solver%converged proc~ridders->proc~converged proc~solution root_module::root_solver%solution proc~ridders->proc~solution

Source Code

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

    implicit none

    class(ridders_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, `-3`=singularity in the algorithm)

    integer  :: i !! counter
    real(wp) :: fh,fl,fm,fnew,denom,xh,xl,xm,xnew

    ! initialize:
    iflag = 0
    fl    = fax
    fh    = fbx
    xl    = ax
    xh    = bx
    xzero = huge(1.0_wp)

    do i = 1, me%maxiter

        xm = bisect(xl,xh)
        fm = me%f(xm)
        if (me%solution(xm,fm,xzero,fzero)) return

        denom = sqrt(fm**2-fl*fh)
        if (denom == 0.0_wp) then
            xzero = xm
            fzero = fm
            iflag = -3        ! can't proceed: denominator is zero [TODO: add a bisection if this happens]
            exit
        end if

        xnew = xm + (xm-xl)*(sign(1.0_wp,fl-fh)*fm/denom)
        if (me%converged(xzero,xnew)) then  ! relative convergence in x
            ! additional check to prevent false convergence
            if (me%converged(xl,xm) .or. me%converged(xm,xh)) exit
        end if

        xzero = xnew
        fnew  = me%f(xzero)
        fzero = fnew
        if (abs(fnew) <= me%ftol) exit    ! abs convergence in f

        ! to keep the root bracketed:
        if (sign(fm,fnew) /= fm) then
            xl = xm
            fl = fm
            xh = xzero
            fh = fnew
        else if (sign(fl,fnew) /= fl) then
            xh = xzero
            fh = fnew
        else if (sign(fh,fnew) /= fh) then
            xl = xzero
            fl = fnew
        end if

        if (me%converged(xl,xh)) exit    ! relative convergence in x
        if (i == me%maxiter) iflag = -2  ! max iterations exceeded

    end do

    end subroutine ridders