interpolate Subroutine

private recursive subroutine interpolate(me, t, x, istatus)

Dense output: compute the state at a time within the last accepted step, using the continuous extension of the method (if it has one, see rklib_properties dense).

This is meant to be called from the report function during an integration (to get the state at times in between the reported points), or after an integration (to get the state within the final step). t must be in the interval of the last accepted step: [t_prev, t], where t is the last reported time and t_prev is the one before it (see dense_interval).

The dense output must have been enabled when the integrator was initialized (see initialize_variable_step), otherwise an error is raised.

This must not be called from the derivative function f. If it is, an error is raised (RKLIB_ERROR_DENSE_NOT_AVAILABLE during a step, or RKLIB_ERROR_DENSE_RECURSIVE while computing the extra function evaluations needed by some methods for the dense output).

If the user has stopped the integration (see rk_class_stop), nothing is computed (no function evaluations are done), and x is returned as zero (with istatus, if present, set to RKLIB_ERROR_USER_STOPPED). This is also the case if the user stops the integration in f while the extra function evaluations for the dense output are computed.

If istatus is present, it returns the status of this call (RKLIB_ERROR_NONE if x was computed), regardless of any previous errors.

Type Bound

rk_class

Arguments

Type IntentOptional Attributes Name
class(rk_class), intent(inout) :: me
real(kind=wp), intent(in) :: t

time at which to compute the state
real(kind=wp), intent(out), dimension(me%n) :: x

state at t
integer, intent(out), optional :: istatus

status of this call

Calls

proc~~interpolate~~CallsGraph proc~interpolate rk_class%interpolate proc~dense_extra_stages rk_class%dense_extra_stages proc~interpolate->proc~dense_extra_stages proc~raise_exception rk_class%raise_exception proc~interpolate->proc~raise_exception

Called by

proc~~interpolate~~CalledByGraph proc~interpolate rk_class%interpolate proc~export_dense_output rk_variable_step_class%export_dense_output proc~export_dense_output->proc~interpolate proc~integrate_to_event_variable_step rk_variable_step_class%integrate_to_event_variable_step proc~integrate_to_event_variable_step->proc~interpolate proc~integrate_to_event_variable_step->proc~export_dense_output proc~export_step rk_variable_step_class%export_step proc~integrate_to_event_variable_step->proc~export_step proc~export_step->proc~export_dense_output proc~integrate_variable_step rk_variable_step_class%integrate_variable_step proc~integrate_variable_step->proc~export_dense_output proc~integrate_variable_step->proc~export_step program~rklib_example rklib_example program~rklib_example->proc~integrate_variable_step

Source Code

    recursive subroutine interpolate(me,t,x,istatus)
        class(rk_class),intent(inout) :: me
        real(wp),intent(in) :: t !! time at which to compute the state
        real(wp),dimension(me%n),intent(out) :: x !! state at `t`
        integer,intent(out),optional :: istatus !! status of this call

        real(wp) :: theta !! normalized time within the step [0,1]
        real(wp) :: theta_tol !! tolerance for the `theta` range check
        integer :: j !! counter
        integer :: degree !! degree of the interpolating polynomial

        if (present(istatus)) istatus = RKLIB_ERROR_NONE

        if (me%dense_computing) then
            ! called from `f` while computing the extra stages for the dense output
            call fail(RKLIB_ERROR_DENSE_RECURSIVE)
            return
        end if

        if (me%stopped) then
            call user_stopped()
            return
        end if

        if (.not. allocated(me%dense_b)) then
            call fail(RKLIB_ERROR_NO_DENSE_OUTPUT)
            return
        end if
        if (.not. me%dense_enabled) then
            call fail(RKLIB_ERROR_DENSE_NOT_ENABLED)
            return
        end if
        if (.not. me%dense_available) then
            ! no step has been taken yet (e.g., the first reported point),
            ! or called from `f` during a step
            call fail(RKLIB_ERROR_DENSE_NOT_AVAILABLE)
            return
        end if

        ! allow for roundoff in the times of the step endpoints:
        theta_tol = 10*spacing(max(abs(me%dense_t),abs(me%dense_t+me%dense_h))) / abs(me%dense_h)
        theta = (t - me%dense_t) / me%dense_h
        if (theta < -theta_tol .or. theta > 1.0_wp + theta_tol) then
            call fail(RKLIB_ERROR_DENSE_OUT_OF_RANGE)
            return
        end if

        if (.not. me%dense_computed) then
            ! compute the polynomial coefficients for this step (only done once per step)
            me%dense_computing = .true.
            call me%dense_extra_stages()
            me%dense_computing = .false.
            if (me%stopped) then
                ! the user stopped the integration in `f`
                call user_stopped()
                return
            end if
            me%dense_r = me%dense_h * matmul(me%funcs(:,1:size(me%dense_b,1)), me%dense_b)
            me%dense_computed = .true.
        end if

        ! evaluate the interpolant (from the inside out):
        degree = size(me%dense_r,2)
        x = me%dense_r(:,degree)
        select case (me%dense_form)
        case (dense_form_polynomial)
            ! x0 + theta*(r1 + theta*(r2 + theta*(r3 + ...)))
            do j = degree-1, 1, -1
                x = me%dense_r(:,j) + theta*x
            end do
        case (dense_form_dopri)
            ! x0 + theta*(r1 + theta1*(r2 + theta*(r3 + theta1*(r4 + ...))))
            do j = degree-1, 1, -1
                if (mod(j,2)==1) then
                    x = me%dense_r(:,j) + (1.0_wp-theta)*x
                else
                    x = me%dense_r(:,j) + theta*x
                end if
            end do
        end select
        x = me%dense_x + theta*x

    contains

        subroutine fail(error_code)
            !! the interpolation failed
            integer,intent(in) :: error_code
            x = zero
            call me%raise_exception(error_code)
            if (present(istatus)) istatus = error_code
        end subroutine fail

        subroutine user_stopped()
            !! the user has stopped the integration (not an error, and
            !! the status has already been set by [[rk_class_stop]])
            x = zero
            if (present(istatus)) istatus = RKLIB_ERROR_USER_STOPPED
        end subroutine user_stopped

    end subroutine interpolate