rbp Subroutine

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

Regula Falsi-Bisection-Parabolic Algorithm (RBP)

See also

  • Alojz Suhadolnik, "Combined bracketing methods for solving nonlinear equations", Applied Mathematics Letters, Volume 25, Issue 11, November 2012, Pages 1755-1760. https://www.sciencedirect.com/science/article/pii/S0893965912000778

Note

This differs from the reference in that it will fallback to bisection if the parabolic interpolation fails.

Type Bound

rbp_solver

Arguments

Type IntentOptional Attributes Name
class(rbp_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~~rbp~~CallsGraph proc~rbp rbp_solver%rbp f f proc~rbp->f proc~bisect bisect proc~rbp->proc~bisect proc~converged root_solver%converged proc~rbp->proc~converged proc~parabolic parabolic proc~rbp->proc~parabolic proc~regula_falsi_step regula_falsi_step proc~rbp->proc~regula_falsi_step proc~solution root_solver%solution proc~rbp->proc~solution proc~parabolic->proc~bisect proc~regula_falsi_step->proc~bisect

Source Code

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

    implicit none

    class(rbp_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) :: xa,xb,xc,x0p,fa,fb,fc,xp,fp,df,dx
    integer :: i  !! iteration counter
    logical :: root_found

    iflag = 0
    xa = ax
    xb = bx
    xc = bisect(xa,xb)
    x0p = xc

    ! function values:
    fa = fax
    fb = fbx
    fc = me%f(xc)
    if (me%solution(xc,fc,xzero,fzero)) return

    do i = 1, me%maxiter

        ! parabolic interpolation:
        xp = parabolic(xa,xb,xc,fa,fb,fc)
        fp = me%f(xp)
        if (me%solution(xp,fp,xzero,fzero)) return

        ! we can check for convergence here.
        ! since everything below is used for the next iteration
        root_found = me%converged(x0p,xp) .and. i>1
        if (root_found .or. i==me%maxiter) then
            xzero = xp
            fzero = fp
            if (.not. root_found) iflag = -2  ! max iterations reached
            exit
        end if

        ! new endpoints of the interval:
        if (fa*fp < 0.0_wp) then
            xb = xp
            fb = fp
            if (fa*fc > 0.0_wp) then
                xa = xc
                fa = fc
            end if
        else
            xa = xp
            fa = fp
            if (fb*fc > 0.0_wp) then
                xb = xc
                fb = fc
            end if
        end if

        ! switching mechanism between bisection and regula falsi:
        df = abs(fa-fb)
        dx = abs(xa-xb)
        if (df>10.0_wp*dx .or. df<0.1_wp*dx) then
            xc = bisect(xa,xb)
        else
            xc = regula_falsi_step(xa,xb,fa,fb,ax,bx)
        end if
        fc = me%f(xc)
        if (me%solution(xc,fc,xzero,fzero)) return

        x0p = xp

    end do

    end subroutine rbp