pegasus Subroutine

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

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

See also

  • G.E. Mullges & F. Uhlig, "Numerical Algorithms with Fortran", Springer, 1996. Section 2.8.2, p 35.

Type Bound

pegasus_solver

Arguments

Type IntentOptional Attributes Name
class(pegasus_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~~pegasus~~CallsGraph proc~pegasus root_module::pegasus_solver%pegasus f f proc~pegasus->f proc~choose_best root_module::choose_best proc~pegasus->proc~choose_best proc~converged root_module::root_solver%converged proc~pegasus->proc~converged proc~secant root_module::secant proc~pegasus->proc~secant proc~solution root_module::root_solver%solution proc~pegasus->proc~solution proc~bisect root_module::bisect proc~secant->proc~bisect

Source Code

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

    implicit none

    class(pegasus_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)

    integer :: i !! counter
    real(wp) :: x1,x2,x3,f1,f2,f3,f1tmp,denom

    ! initialize:
    iflag = 0
    x1    = ax
    x2    = bx
    f1    = fax
    f2    = fbx

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

        ! secant step
        x3 = secant(x1,x2,f1,f2,ax,bx)

        f3  = me%f(x3)  ! calculate f3
        if (me%solution(x3,f3,xzero,fzero)) return

        ! determine a new inclusion interval:
        if (f2*f3<0.0_wp) then  ! root on (x2,x3)
            x1 = x2
            f1 = f2
            f1tmp = f1
        else  ! root on (x1,x3)
            f1tmp = f1
            denom = f2 + f3
            if (denom /= 0.0_wp) then
                ! proceed as normal
                f1 = f1 * f2 / denom
            else
                ! can't proceed, keep as is.
                ! [need a find a test case where this happens -TODO]
            end if
        end if

        x2 = x3
        f2 = f3

        call choose_best(x1,x2,f1tmp,f2,xzero,fzero)

        if (me%converged(x1,x2)) exit   ! check for convergence
        if (i == me%maxiter) iflag = -2 ! max iterations exceeded

    end do

    end subroutine pegasus