blendtf Subroutine

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

BlendTF blended method of trisection and false position methods.

Reference

  • E Badr, S Almotairi, A El Ghamry, "A Comparative Study among New Hybrid Root Finding Algorithms and Traditional Methods", Mathematics 2021, 9, 1306.

Type Bound

blendtf_solver

Arguments

Type IntentOptional Attributes Name
class(blendtf_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~~blendtf~~CallsGraph proc~blendtf root_module::blendtf_solver%blendtf f f proc~blendtf->f proc~choose_best root_module::choose_best proc~blendtf->proc~choose_best proc~converged root_module::root_solver%converged proc~blendtf->proc~converged proc~solution root_module::root_solver%solution proc~blendtf->proc~solution

Source Code

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

    implicit none

    class(blendtf_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) :: a1,a2,b1,b2,fa,fb,a,b,xt1,xt2,xf,x,fx,fxt2,&
                fxf,xprev,fxt1,fxprev,fa1,fa2,fb1,fb2
    integer :: i !! iteration counter

    iflag = 0
    a = ax; b = bx
    fa = fax; fb = fbx
    a1 = a; a2 = a
    b1 = b; b2 = b
    fa1 = fa; fa2 = fa
    fb1 = fb; fb2 = fb
    xprev = huge(1.0_wp)
    fxprev = huge(1.0_wp)

    do i = 1, me%maxiter

        if (fa == fb) then
            ! should fallback to bisection if this happens -TODO
            iflag = -3
            return
        end if

        xt1  = (b + 2.0_wp * a) / 3.0_wp
        xt2  = (2.0_wp * b + a) / 3.0_wp
        xf   = a - (fa*(b-a))/(fb-fa)
        x    = xt1
        fxt1 = me%f(xt1); if (me%solution(xt1,fxt1,xzero,fzero)) return
        fxt2 = me%f(xt2); if (me%solution(xt2,fxt2,xzero,fzero)) return
        fxf  = me%f(xf);  if (me%solution(xf,fxf,xzero,fzero))   return
        fx   = fxt1

        if (abs(fxf) < abs(fxt1)) then
            x = xf
            fx = fxf
        elseif (abs(fxt2) < abs(fxt1)) then
            x = xt2
            fx = fxt2
        end if

        ! apply the convergence tols to [a,b]
        if (me%converged(a, b) .or. i == me%maxiter) then
            call choose_best(x,xprev,fx,fxprev,xzero,fzero)
            if (i == me%maxiter) iflag = -2 ! max iterations reached
            exit
        end if
        xprev = x
        fxprev = fx

        if ((fa * fxt1) < 0.0_wp) then
            b1 = xt1; fb1 = fxt1
        else if ((fxt1 * fxt2) < 0.0_wp) then
            a1 = xt1; fa1 = fxt1
            b1 = xt2; fb1 = fxt2
        else
            a1 = xt2; fa1 = fxt2
        end if

        if (fa*fxf < 0.0_wp) then
            b2 = xf; fb2 = fxf
        else
            a2 = xf; fa2 = fxf
        end if

        if (a1>a2) then
            a = a1
            fa = fa1
        else
            a = a2
            fa = fa2
        end if
        if (b1<b2) then
            b = b1
            fb = fb1
        else
            b = b2
            fb = fb2
        end if

    end do

    end subroutine blendtf