zero Function

public function zero(obj, a, b, machep, t)

Seeks the root of a function F(X) in an interval [A,B].

The interval [A,B] must be a change of sign interval for F. That is, F(A) and F(B) must be of opposite signs. Then assuming that F is continuous implies the existence of at least one value C between A and B for which F(C) = 0.

The location of the zero is determined to within an accuracy of 6 * MACHEPS * abs ( C ) + 2 * T.

Thanks to Thomas Secretin for pointing out a transcription error in the setting of the value of P, 11 February 2013.

Modified: 11 February 2013

Author:

  • Original FORTRAN77 version by Richard Brent.
  • FORTRAN90 version by John Burkardt.

Reference:

  • Richard Brent,
  • Algorithms for Minimization Without Derivatives,
  • Dover, 2002,
  • ISBN: 0-486-41998-3,
  • LC: QA402.5.B74.

Arguments

Type IntentOptional Attributes Name
class(brent_t), intent(in) :: obj

external object with the user-supplied subroutine which calculates the functions.
real(kind=dp), intent(in) :: a

the endpoints of the change of sign interval
real(kind=dp), intent(in) :: b

the endpoints of the change of sign interval
real(kind=dp), intent(in) :: machep

an estimate for the relative machine precision.
real(kind=dp), intent(in) :: t

a positive error tolerance

Return Value real(kind=dp)

the estimated value of a zero of the function F


Calls

proc~~zero~~CallsGraph proc~zero zero f f proc~zero->f

Called by

proc~~zero~~CalledByGraph proc~zero zero proc~brent brent proc~brent->proc~zero proc~hugoniot_star_state hugoniot_star_state proc~hugoniot_star_state->proc~zero proc~droeint_drop droeint_droP proc~hugoniot_star_state->proc~droeint_drop proc~outgoing_branch_cold_circulator outgoing_branch_cold_circulator proc~outgoing_branch_cold_circulator->proc~zero proc~t_rop T_roP proc~outgoing_branch_cold_circulator->proc~t_rop proc~circulator_dynamic_parameters_df circulator_dynamic_parameters_t%circulator_dynamic_parameters_df proc~outgoing_branch_cold_circulator->proc~circulator_dynamic_parameters_df proc~state_rop state_roP proc~outgoing_branch_cold_circulator->proc~state_rop proc~circulator_dynamic_parameters_f circulator_dynamic_parameters_t%circulator_dynamic_parameters_f proc~outgoing_branch_cold_circulator->proc~circulator_dynamic_parameters_f proc~circulator_resolution_from_and_to_ports circulator_resolution_from_and_to_ports proc~circulator_resolution_from_and_to_ports->proc~outgoing_branch_cold_circulator proc~incoming_branch_cold_circulator incoming_branch_cold_circulator proc~circulator_resolution_from_and_to_ports->proc~incoming_branch_cold_circulator proc~ro_pt ro_pT proc~ro_pt->proc~brent proc~solve_boundary_pt solve_boundary_PT proc~solve_boundary_pt->proc~hugoniot_star_state proc~solve_boundary_pt->proc~ro_pt proc~derivative_function_of_pstar_pt derivative_function_of_pstar_PT proc~solve_boundary_pt->proc~derivative_function_of_pstar_pt proc~solve_boundary_pt->proc~state_rop proc~state_rop_withr state_roP_withR proc~solve_boundary_pt->proc~state_rop_withr proc~function_of_pstar_pt function_of_pstar_PT proc~solve_boundary_pt->proc~function_of_pstar_pt proc~t_roe T_roE proc~t_roe->proc~brent proc~t_rop->proc~brent proc~boundary_resolution_from_and_to_ports boundary_resolution_from_and_to_ports proc~boundary_resolution_from_and_to_ports->proc~solve_boundary_pt proc~solve_boundary_mt solve_boundary_MT proc~boundary_resolution_from_and_to_ports->proc~solve_boundary_mt proc~channel_init_part1 channel_init_part1 proc~channel_init_part1->proc~ro_pt proc~circulator_dynamic_parameters_df->proc~t_rop proc~derivative_function_of_pstar_pt->proc~t_rop proc~droeint_drop->proc~t_rop proc~ffsrclink_prelax_resolution_from_and_to_ports FFsrcLink_Prelax_resolution_from_and_to_ports proc~ffsrclink_prelax_resolution_from_and_to_ports->proc~t_rop proc~function_of_pstar_mt function_of_pstar_MT proc~function_of_pstar_mt->proc~t_rop proc~main_loop main_loop proc~main_loop->proc~circulator_resolution_from_and_to_ports proc~main_loop->proc~boundary_resolution_from_and_to_ports proc~main_loop->proc~ffsrclink_prelax_resolution_from_and_to_ports proc~friction_source_term friction_source_term proc~main_loop->proc~friction_source_term proc~from_sol_to_prim_channel from_sol_to_prim_channel proc~main_loop->proc~from_sol_to_prim_channel proc~junction_resolution_from_and_to_ports Junction_resolution_from_and_to_ports proc~main_loop->proc~junction_resolution_from_and_to_ports proc~source_term_definition_channel source_term_definition_channel proc~main_loop->proc~source_term_definition_channel proc~full_physics_definition_channel full_physics_definition_channel proc~main_loop->proc~full_physics_definition_channel proc~solve_boundary_mt->proc~t_rop proc~solve_boundary_mt->proc~function_of_pstar_mt proc~solve_boundary_mt->proc~state_rop proc~solve_junction solve_junction proc~solve_junction->proc~t_roe proc~solve_junction->proc~droeint_drop proc~solve_junction->proc~state_rop proc~state_roe state_roE proc~state_roe->proc~t_roe proc~state_roe_withr state_roE_withR proc~state_roe_withr->proc~t_rop proc~state_rop->proc~t_rop proc~state_rop_withr->proc~t_rop proc~circulator_dynamic_parameters_f->proc~state_rop proc~friction_source_term->proc~state_roe proc~friction_source_term->proc~state_roe_withr proc~from_sol_to_prim_channel->proc~state_roe proc~from_sol_to_prim_channel->proc~state_roe_withr proc~function_of_pstar_pt->proc~state_rop proc~hug_brent_f hug_brent_t%hug_brent_f proc~hug_brent_f->proc~droeint_drop proc~incoming_branch_cold_circulator->proc~state_rop proc~junction_dynamic_parameters_fcn junction_dynamic_parameters_t%junction_dynamic_parameters_fcn proc~junction_dynamic_parameters_fcn->proc~droeint_drop proc~junction_resolution_from_and_to_ports->proc~solve_junction proc~sourceterms_notfriction_channel SourceTerms_NOTfriction_channel proc~sourceterms_notfriction_channel->proc~state_roe proc~sourceterms_notfriction_channel->proc~state_roe_withr program~reims_p reims_p program~reims_p->proc~channel_init_part1 program~reims_p->proc~main_loop proc~source_term_definition_channel->proc~sourceterms_notfriction_channel proc~full_physics_definition_channel->proc~source_term_definition_channel

Source Code

function zero(obj, a, b, machep, t)
    !! Seeks the root of a function F(X) in an interval [A,B].
    !!
    !! The interval [A,B] must be a change of sign interval for F.
    !! That is, F(A) and F(B) must be of opposite signs.  Then
    !! assuming that F is continuous implies the existence of at least
    !! one value C between A and B for which F(C) = 0.
    !!
    !! The location of the zero is determined to within an accuracy
    !! of 6 * MACHEPS * abs ( C ) + 2 * T.
    !!
    !! Thanks to Thomas Secretin for pointing out a transcription error in the
    !! setting of the value of P, 11 February 2013.
    !!
    !! Modified: 11 February 2013
    !!
    !!  Author:
    !!
    !!    - Original FORTRAN77 version by Richard Brent.
    !!    - FORTRAN90 version by John Burkardt.
    !!
    !!  Reference:
    !!
    !!    - Richard Brent,
    !!    - Algorithms for Minimization Without Derivatives,
    !!    - Dover, 2002,
    !!    - ISBN: 0-486-41998-3,
    !!    - LC: QA402.5.B74.
    class(brent_t), intent(in) :: obj !! external object with the user-supplied 
                                   !! subroutine which calculates the functions.
    real(dp), intent(in) :: a,b    !! the endpoints of the change of sign interval
    real(dp), intent(in) :: machep !! an estimate for the relative machine precision.
    real(dp), intent(in) :: t      !! a positive error tolerance
    real(dp) :: zero   !! the estimated value of a zero of the function F

    real(dp) :: c
    real(dp) :: d
    real(dp) :: e
    real(dp) :: fa
    real(dp) :: fb
    real(dp) :: fc
    real(dp) :: m
    real(dp) :: p
    real(dp) :: q
    real(dp) :: r
    real(dp) :: s
    real(dp) :: sa
    real(dp) :: sb
    real(dp) :: tol

    ! Make local copies of A and B.
    sa = a
    sb = b
    fa = obj%f(sa)
    fb = obj%f(sb)
  
    c = sa
    fc = fa
    e = sb - sa
    d = e
    do    
        if (abs(fc) < abs(fb)) then    
            sa = sb
            sb = c
            c = sa
            fa = fb
            fb = fc
            fc = fa    
        endif
  
        tol = 2.0_dp * machep * abs(sb) + t
        m = 0.5_dp * (c - sb)
        if (abs(m) <= tol .or. fb == 0.0_dp) then
            exit
        endif
  
        if (abs(e) < tol .or. abs (fa) <= abs(fb)) then    
            e = m
            d = e    
        else    
            s = fb / fa    
            if (sa == c) then    
                p = 2.0_dp * m * s
                q = 1.0_dp - s    
            else    
                q = fa / fc
                r = fb / fc
                p = s * (2.0_dp * m * q * (q - r) - (sb - sa) * (r - 1.0_dp))
                q = (q - 1.0_dp) * (r - 1.0_dp) * (s - 1.0_dp)    
            endif
  
            if (0.0_dp < p) then
                q = - q
            else
                p = - p
            endif
  
            s = e
            e = d
            if (2.0_dp * p < 3.0_dp * m * q - abs(tol * q) .and. &
                                    p < abs (0.5_dp * s * q)) then
                d = p / q
            else
                e = m
                d = e
            endif
        endif
  
        sa = sb
        fa = fb
        if (tol < abs(d)) then
            sb = sb + d
        else if (0.0_dp < m) then
            sb = sb + tol
        else
            sb = sb - tol
        endif
  
        fb = obj%f(sb)
        if ((0.0_dp < fb .and. 0.0_dp < fc ) .or. &
                (fb <= 0.0_dp .and. fc <= 0.0_dp)) then
            c = sa
            fc = fa
            e = sb - sa
            d = e
        endif
    enddo
    zero = sb
end function zero