hinit Function

private function hinit(me, x, y, posneg, f0, hmax, atol, rtol)

computation of an initial step size guess

Note

This routine is from dop853. It was modified for this module.

Type Bound

rk_variable_step_class

Arguments

Type IntentOptional Attributes Name
class(rk_variable_step_class), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in), dimension(:) :: y

dimension(n)
real(kind=wp), intent(in) :: posneg

posneg = sign(1.0_wp,xend-x)
real(kind=wp), intent(in), dimension(:) :: f0

dimension(n)
real(kind=wp), intent(in) :: hmax
real(kind=wp), intent(in), dimension(:) :: atol
real(kind=wp), intent(in), dimension(:) :: rtol

Return Value real(kind=wp)


Calls

proc~~hinit~~CallsGraph proc~hinit rklib_module::rk_variable_step_class%hinit f f proc~hinit->f proc~order rklib_module::rk_variable_step_class%order proc~hinit->proc~order properties properties proc~order->properties

Called by

proc~~hinit~~CalledByGraph proc~hinit rklib_module::rk_variable_step_class%hinit proc~compute_initial_step rklib_module::rk_variable_step_class%compute_initial_step proc~compute_initial_step->proc~hinit proc~integrate_to_event_variable_step rklib_module::rk_variable_step_class%integrate_to_event_variable_step proc~integrate_to_event_variable_step->proc~compute_initial_step proc~integrate_variable_step rklib_module::rk_variable_step_class%integrate_variable_step proc~integrate_variable_step->proc~compute_initial_step program~rklib_example rklib_example program~rklib_example->proc~integrate_variable_step

Source Code

    function hinit(me,x,y,posneg,f0,hmax,atol,rtol)

    implicit none

    class(rk_variable_step_class),intent(inout) :: me
    real(wp),intent(in)               :: x
    real(wp),dimension(:),intent(in)  :: y       !! dimension(n)
    real(wp),intent(in)               :: posneg  !! posneg = sign(1.0_wp,xend-x)
    real(wp),dimension(:),intent(in)  :: f0      !! dimension(n)
    real(wp),intent(in)               :: hmax
    real(wp),dimension(:),intent(in)  :: atol
    real(wp),dimension(:),intent(in)  :: rtol

    real(wp) :: der12,der2,dnf,dny,h,h1,hinit,sk
    integer :: i
    integer :: iord  !! order of the method
    real(wp),dimension(me%n) :: f1,y1

    iord = me%order()

    ! compute a first guess for explicit euler as
    !   h = 0.01 * norm (y0) / norm (f0)
    ! the increment for explicit euler is small
    ! compared to the solution
    dnf = zero
    dny = zero
    do i = 1 , me%n
        sk = atol(i) + rtol(i)*abs(y(i))
        dnf = dnf + (f0(i)/sk)**2
        dny = dny + (y(i)/sk)**2
    end do
    if ( dnf<=1.0e-10_wp .or. dny<=1.0e-10_wp ) then
        h = 1.0e-6_wp
    else
        h = sqrt(dny/dnf)*0.01_wp
    end if
    h = min(h,hmax)
    h = sign(h,posneg)
    ! perform an explicit euler step
    do i = 1 , me%n
        y1(i) = y(i) + h*f0(i)
    end do
    call me%f(x+h,y1,f1)
    ! estimate the second derivative of the solution
    der2 = zero
    do i = 1 , me%n
        sk = atol(i) + rtol(i)*abs(y(i))
        der2 = der2 + ((f1(i)-f0(i))/sk)**2
    end do
    der2 = sqrt(der2)/h
    ! step size is computed such that
    !  h**iord * max ( norm (f0), norm (der2)) = 0.01
    der12 = max(abs(der2),sqrt(dnf))
    if ( der12<=1.0e-15_wp ) then
        h1 = max(1.0e-6_wp,abs(h)*1.0e-3_wp)
    else
        h1 = (0.01_wp/der12)**(1.0_wp/iord)
    end if

    h = min(100.0_wp*abs(h),h1,hmax)
    hinit = sign(h,posneg)

    end function hinit