bdqrf Subroutine

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

Bisected Direct Quadratic Regula Falsi (BDQRF) root solver method to find the root of a 1D function.

See also

  • R. G. Gottlieb, B. F. Thompson, "Bisected Direct Quadratic Regula Falsi", Applied Mathematical Sciences, Vol. 4, 2010, no. 15, 709-718.

Type Bound

bdqrf_solver

Arguments

Type IntentOptional Attributes Name
class(bdqrf_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~~bdqrf~~CallsGraph proc~bdqrf root_module::bdqrf_solver%bdqrf f f proc~bdqrf->f proc~bisect root_module::bisect proc~bdqrf->proc~bisect proc~choose_best root_module::choose_best proc~bdqrf->proc~choose_best proc~converged root_module::root_solver%converged proc~bdqrf->proc~converged proc~solution root_module::root_solver%solution proc~bdqrf->proc~solution

Source Code

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

    implicit none

    class(bdqrf_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) :: xdn,ydn,xup,yup,d,xm,ym,a,b,y2
    integer :: i !! counter

    ! initialize:
    iflag = 0
    xzero = ax
    fzero = fax
    y2    = fbx

    if (fzero<0.0_wp) then
        xdn = ax
        ydn = fzero
        xup = bx
        yup = y2
    else
        xup = ax
        yup = fzero
        xdn = bx
        ydn = y2
    end if

    ! main loop:
    do i = 1, me%maxiter

        xm = bisect(xup,xdn)
        ym = me%f(xm)
        if (me%solution(xm,ym,xzero,fzero)) return ! Convergence

        d = (xup - xdn) / 2.0_wp
        a = (yup + ydn - 2.0_wp*ym)/(2.0_wp*d**2)
        b = (yup - ydn)/(2.0_wp*d)

        xzero = xm - 2.0_wp*ym / (b * (1.0_wp + sqrt(1.0_wp - 4.0_wp*a*ym/b**2)))
        fzero = me%f(xzero)
        if (me%solution(xzero,fzero,xzero,fzero)) return ! Convergence

        if (fzero>0.0_wp) then
            yup = fzero
            xup = xzero
            if (ym<0.0_wp) then
                ydn = ym
                xdn = xm
            end if
        else
            ydn = fzero
            xdn = xzero
            if (ym>0.0_wp) then
                yup = ym
                xup = xm
            end if
        end if

        if (me%converged(xdn,xup) .or. i==me%maxiter) then
            call choose_best(xdn,xup,ydn,yup,xzero,fzero)
            if (i==me%maxiter) iflag = -2 ! maximum number of iterations
            exit
        end if

    end do

    end subroutine bdqrf