zeroin Subroutine

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

Uses

  • proc~~zeroin~~UsesGraph proc~zeroin brent_module::brent_class%zeroin iso_fortran_env iso_fortran_env proc~zeroin->iso_fortran_env

Find a zero of the function in the given interval to within a tolerance , where is the relative machine precision defined as the smallest representable number such that .

It is assumed that and have opposite signs.

References

See also

  1. zeroin.f from Netlib

Type Bound

brent_class

Arguments

Type IntentOptional Attributes Name
class(brent_class), 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) :: tol

desired length of the interval of uncertainty of the final result (>=0)
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 (-1=error, 0=root found)
real(kind=wp), intent(in), optional :: fax

if f(ax) is already known, it can be input here
real(kind=wp), intent(in), optional :: fbx

if f(ax) is already known, it can be input here

Called by

proc~~zeroin~~CalledByGraph proc~zeroin brent_module::brent_class%zeroin proc~brent_test brent_module::brent_test proc~brent_test->proc~zeroin proc~integrate_to_event rk_module_variable_step::rk_variable_step_class%integrate_to_event proc~integrate_to_event->proc~zeroin proc~integrate_to_event~2 rk_module::rk_class%integrate_to_event proc~integrate_to_event~2->proc~zeroin proc~halo_to_rv_diffcorr halo_orbit_module::halo_to_rv_diffcorr proc~halo_to_rv_diffcorr->proc~integrate_to_event~2 proc~rk_test rk_module::rk_test proc~rk_test->proc~integrate_to_event~2 proc~rk_test_variable_step rk_module_variable_step::rk_test_variable_step proc~rk_test_variable_step->proc~integrate_to_event

Source Code

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

    use iso_fortran_env, only: error_unit

    implicit none

    class(brent_class),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)              :: tol     !! desired length of the interval of uncertainty of the final result (>=0)
    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 (`-1`=error, `0`=root found)
    real(wp),intent(in),optional     :: fax     !! if `f(ax)` is already known, it can be input here
    real(wp),intent(in),optional     :: fbx     !! if `f(ax)` is already known, it can be input here

    real(wp),parameter :: eps   = epsilon(one)  !! original code had d1mach(4)
    real(wp) :: a,b,c,d,e,fa,fb,fc,tol1,xm,p,q,r,s

    tol1 = eps+one

    a=ax
    b=bx

    if (present(fax)) then
        fa = fax
    else
        fa=me%f(a)
    end if
    if (present(fbx)) then
        fb = fbx
    else
        fb=me%f(b)
    end if

    !check trivial cases first:
    if (fa==zero) then

        iflag = 0
        xzero = a
        fzero = fa

    elseif (fb==zero) then

        iflag = 0
        xzero = b
        fzero = fb

    elseif (fa*(fb/abs(fb))<zero) then  ! check that f(ax) and f(bx) have different signs

        c=a
        fc=fa
        d=b-a
        e=d

        do

            if (abs(fc)<abs(fb)) then
                a=b
                b=c
                c=a
                fa=fb
                fb=fc
                fc=fa
            end if

            tol1=two*eps*abs(b)+0.5_wp*tol
            xm = 0.5_wp*(c-b)
            if ((abs(xm)<=tol1).or.(fb==zero)) exit

            ! see if a bisection is forced
            if ((abs(e)>=tol1).and.(abs(fa)>abs(fb))) then
                s=fb/fa
                if (a/=c) then
                    ! inverse quadratic interpolation
                    q=fa/fc
                    r=fb/fc
                    p=s*(two*xm*q*(q-r)-(b-a)*(r-one))
                    q=(q-one)*(r-one)*(s-one)
                else
                    ! linear interpolation
                    p=two*xm*s
                    q=one-s
                end if
                if (p<=zero) then
                    p=-p
                else
                    q=-q
                end if
                s=e
                e=d
                if (((two*p)>=(three*xm*q-abs(tol1*q))) .or. &
                    (p>=abs(0.5_wp*s*q))) then
                    d=xm
                    e=d
                else
                    d=p/q
                end if
            else
                d=xm
                e=d
            end if

            a=b
            fa=fb
            if (abs(d)<=tol1) then
                if (xm<=zero) then
                    b=b-tol1
                else
                    b=b+tol1
                end if
            else
                b=b+d
            end if
            fb=me%f(b)
            if ((fb*(fc/abs(fc)))>zero) then
                c=a
                fc=fa
                d=b-a
                e=d
            end if

        end do

        iflag = 0
        xzero = b
        fzero = fb

    else

        iflag = -1
        write(error_unit,'(A)')&
            'Error in zeroin: f(ax) and f(bx) do not have different signs.'

    end if

    end subroutine zeroin