rk_variable_step_class

type, public, abstract, extends(rk_class) :: rk_variable_step_class

Main integration class for variable step size Runge-Kutta methods

Inherits

Help

Components

Type	Visibility	Attributes		Name		Initial
type(stepsize_class),	private		::	stepsize_method			the method for varying the step size
real(kind=wp),	private,	dimension(:), allocatable	::	rtol			relative tolerance (`size(n)`)
real(kind=wp),	private,	dimension(:), allocatable	::	atol			absolute tolerance (`size(n)`)
integer,	private		::	hinit_method	=	1	if automatically computing the inital step size, which method to use. 1 = `hstart`, 2 = `hinit`.
integer,	private		::	num_rejected_steps	=	0	number of rejected steps
real(kind=wp),	private		::	last_accepted_step_size	=	zero	the last accepted step size `dt` from the integration (positive or negative)
real(kind=wp),	private		::	dense_output_step	=	zero	if > 0, the points are reported at this fixed time interval using the dense output (see initialize_variable_step)
logical,	private		::	dense_reporting	=	.false.	if the dense output is used for reporting in the current integration
real(kind=wp),	private		::	dense_output_t0	=	zero	initial time of the current integration
integer(kind=int64),	private		::	dense_output_k	=	0	index of the last dense output point reported: `t = dense_output_t0 + k*dense_output_step` (with the sign for the direction of integration)
real(kind=wp),	private,	allocatable	::	t_saved			cached `t`
real(kind=wp),	private,	dimension(:), allocatable	::	x_saved			cached `x`
real(kind=wp),	private,	dimension(:), allocatable	::	f_saved			cached `f`

Type-Bound Procedures

procedure, public :: destroy

destructor

private subroutine destroy(me)

Destructor for rk_class.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(out)			::	me

procedure, public :: stop => rk_class_stop

user-callable method to stop the integration

private subroutine rk_class_stop(me)

User-callable method to stop the integration (e.g., from f, g, or report). No further function evaluations are done, and the integrator returns, with the status set to RKLIB_ERROR_USER_STOPPED. This is not treated as an error (i.e., it does not stop the program if stop_on_errors is true).

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(inout)			::	me

procedure, public :: status => rk_class_status

get status code and message

private subroutine rk_class_status(me, istatus, message)

Get the status of an integration.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(in)			::	me
integer,	intent(out),	optional		::	istatus	status code (`<0` means an error)
character(len=:),	intent(out),	optional,	allocatable	::	message	status message

procedure, public :: failed

private function failed(me)

Returns true if there was an error. Can use rk_class_status to get more info.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(in)			::	me

Return Value logical

procedure, public :: interpolate

dense output for the last accepted step

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).

Arguments

Type	Intent	Optional	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

procedure, public :: dense_interval

interval of the last accepted step for the dense output

private subroutine dense_interval(me, available, t0, tf, x0, xf)

Returns the interval of the last accepted step, which is where interpolate can be used. For example, this can be called from the report function to get the previous time and state, so that points in between can be computed.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(inout)			::	me
logical,	intent(out)			::	available	if dense output is available
real(kind=wp),	intent(out),	optional		::	t0	initial time of the step
real(kind=wp),	intent(out),	optional		::	tf	final time of the step
real(kind=wp),	intent(out),	optional,	dimension(me%n)	::	x0	initial state of the step
real(kind=wp),	intent(out),	optional,	dimension(me%n)	::	xf	final state of the step

procedure(properties_func), public, deferred :: properties

pure function properties_func(me) result(p) Prototype

Returns the properties of the method.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_class),	intent(in)			::	me

Return Value type(rklib_properties)

properties of the method

procedure(step_func_variable), private, deferred :: step

the step routine for the rk method

subroutine step_func_variable(me, t, x, h, xf, xerr) Prototype

rk step function for the variable-step methods.

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t	initial time
real(kind=wp),	intent(in),	dimension(me%n)	::	x	initial state vector
real(kind=wp),	intent(in)		::	h	time step $\|\Delta t\|$
real(kind=wp),	intent(out),	dimension(me%n)	::	xf	final state vector
real(kind=wp),	intent(out),	dimension(me%n)	::	xerr	truncation error estimate

procedure, public :: initialize => initialize_variable_step

initialize the class (set n,f, and report)

private subroutine initialize_variable_step(me, n, f, rtol, atol, stepsize_method, hinit_method, report, g, stop_on_errors, max_number_of_steps, report_rate, solver, dense_output, dense_output_step)

Initialize the rk_variable_step_class.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_variable_step_class),	intent(inout)			::	me
integer,	intent(in)			::	n	number of equations
procedure(deriv_func)				::	f	derivative function
real(kind=wp),	intent(in),	optional,	dimension(:)	::	rtol	relative tolerance (if size=1, then same tol used for all equations). If not present, a default of `100*eps` is used
real(kind=wp),	intent(in),	optional,	dimension(:)	::	atol	absolute tolerance (if size=1, then same tol used for all equations). If not present, a default of `100*eps` is used
type(stepsize_class),	intent(in),	optional		::	stepsize_method	method for varying the step size
integer,	intent(in),	optional		::	hinit_method	which method (1 or 2) to use for automatic initial step size computation. 1 = use `hstart`, 2 = use `hinit`.
procedure(report_func),		optional		::	report	for reporting the steps
procedure(event_func),		optional		::	g	for stopping at an event
logical,	intent(in),	optional		::	stop_on_errors	stop the program for any errors (default is False)
integer,	intent(in),	optional		::	max_number_of_steps	max number of steps allowed
integer,	intent(in),	optional		::	report_rate	how often to call report function: `0` : no reporting (same as not associating `report`), `1` : report every point, `2` : report every other point, etc. The first and last point are always reported.
class(root_method),	intent(in),	optional		::	solver	the root-finding method to use for even finding. if not present, then `brent_solver` is used.
logical,	intent(in),	optional		::	dense_output	to enable the dense output (default is False, unless `dense_output_step` is present and not zero). If enabled, interpolate can be used (e.g., from the `report` function), and the dense output is also used for event location in integrate_to_event_variable_step. The method must have a dense output. When not enabled, there is no overhead for the dense output.
real(kind=wp),	intent(in),	optional		::	dense_output_step	if present and not zero, the points are reported at this fixed time interval (starting from the initial time), using the dense output of the method, rather than at the steps taken by the integrator. The first and last point are also reported (so the last interval will generally be shorter). This is only used if `report` is present and `report_rate>0` (`report_rate` is otherwise ignored). The method must have a dense output. This automatically enables the dense output (it is an error if `dense_output=.false.` is also specified).

procedure, public :: integrate => integrate_variable_step

private subroutine integrate_variable_step(me, t0, x0, h, tf, xf)

Main integration routine for the rk_variable_step_class.

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t0	initial time
real(kind=wp),	intent(in),	dimension(:)	::	x0	initial state
real(kind=wp),	intent(in)		::	h	initial abs(time step)
real(kind=wp),	intent(in)		::	tf	final time
real(kind=wp),	intent(out),	dimension(:)	::	xf	final state

procedure, public :: integrate_to_event => integrate_to_event_variable_step

private subroutine integrate_to_event_variable_step(me, t0, x0, h, tmax, tol, tf, xf, gf)

Event-finding integration routine for the rk_variable_step_class. Integrates until g(t,x)=0, or until t=tf (whichever happens first).

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t0	initial time
real(kind=wp),	intent(in),	dimension(me%n)	::	x0	initial state
real(kind=wp),	intent(in)		::	h	abs(time step)
real(kind=wp),	intent(in)		::	tmax	max final time if event not located
real(kind=wp),	intent(in)		::	tol	function tolerance for root finding
real(kind=wp),	intent(out)		::	tf	actual final time reached
real(kind=wp),	intent(out),	dimension(me%n)	::	xf	final state (at tf)
real(kind=wp),	intent(out)		::	gf	g value at tf

procedure, public :: info => info_variable_step

private subroutine info_variable_step(me, num_steps, num_rejected_steps, last_accepted_step_size)

Return some info about the integration.

Arguments

Type	Intent	Optional		Name
class(rk_variable_step_class),	intent(in)		::	me
integer,	intent(out),	optional	::	num_steps	number of steps taken
integer,	intent(out),	optional	::	num_rejected_steps	number of rejected steps
real(kind=wp),	intent(out),	optional	::	last_accepted_step_size	the last accepted step size `dt` from the integration (positive or negative)

procedure, private :: hstart

for automatically computing the initial step size [this is from DDEABM]

private subroutine hstart(me, a, b, y, yprime, etol, h)

Computes a starting step size to be used in solving initial value problems in ordinary differential equations.

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	a	the initial point of integration.
real(kind=wp),	intent(in)		::	b	a value of the independent variable used to define the direction of integration. a reasonable choice is to set `b` to the first point at which a solution is desired. you can also use `b`, if necessary, to restrict the length of the first integration step because the algorithm will not compute a starting step length which is bigger than `abs(b-a)`, unless `b` has been chosen too close to `a`. (it is presumed that hstart has been called with `b` different from `a` on the machine being used. also see the discussion about the parameter `small`.)
real(kind=wp),	intent(in),	dimension(me%n)	::	y	the vector of initial values of the `neq` solution components at the initial point `a`.
real(kind=wp),	intent(in),	dimension(me%n)	::	yprime	the vector of derivatives of the `neq` solution components at the initial point `a`. (defined by the differential equations in subroutine `me%f`)
real(kind=wp),	intent(in),	dimension(me%n)	::	etol	the vector of error tolerances corresponding to the `neq` solution components. it is assumed that all elements are positive. following the first integration step, the tolerances are expected to be used by the integrator in an error test which roughly requires that `abs(local error) <= etol` for each vector component.
real(kind=wp),	intent(out)		::	h	appropriate starting step size to be attempted by the differential equation method.

procedure, private :: hinit

for automatically computing the initial step size [this is from DOP853]

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

computation of an initial step size guess

Arguments

Type	Intent	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)

procedure, private :: begin_integration => begin_integration_rk_variable_step_class

private subroutine begin_integration_rk_variable_step_class(me)

Begin a rk_variable_step_class integration.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_variable_step_class),	intent(inout)			::	me

procedure, private :: compute_initial_step

private function compute_initial_step(me, t0, tf, x0, h0) result(dt)

Compute the initial step size.

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t0	initial time
real(kind=wp),	intent(in)		::	tf	final time
real(kind=wp),		dimension(me%n)	::	x0	initial state
real(kind=wp),	intent(in)		::	h0	user-input initial step size (if zero, then one is computed)

Return Value real(kind=wp)

step size to use

procedure, private :: export_step

private subroutine export_step(me, t, x, ok)

Export the point(s) after a step has been accepted (other than the last step). Either the end of the step, or the dense output points within the step (see export_dense_output).

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t	time at the end of the step
real(kind=wp),	intent(in),	dimension(:)	::	x	state at the end of the step
logical,	intent(out)		::	ok	false if there was an error (the integration should not continue)

procedure, private :: export_dense_output

private subroutine export_dense_output(me, t_stop, ok)

Report the dense output points (at the fixed dense_output_step interval) within the last accepted step, that are before t_stop.

Arguments

Type	Intent		Name
class(rk_variable_step_class),	intent(inout)	::	me
real(kind=wp),	intent(in)	::	t_stop	report points up to (but not including) this time
logical,	intent(out)	::	ok	false if there was an error (the integration should not continue)

procedure, private :: order

returns p, the order of the method

private pure function order(me) result(p)

Returns the order of the RK method

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_variable_step_class),	intent(in)			::	me

Return Value integer

order of the method

procedure, public :: destroy_fsal_cache

private subroutine destroy_fsal_cache(me)

Destructor for the function evaluation cache.

Arguments

Type	Intent	Optional	Attributes		Name
class(rk_variable_step_class),	intent(inout)			::	me

procedure, public :: check_fsal_cache

private subroutine check_fsal_cache(me, t, x, f)

Compute the function at t and x, or get it from the function evaluation cache if it is there (see set_fsal_cache). This is used by all the variable-step methods for the first stage of a step. It may already be in the cache because:

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t
real(kind=wp),	intent(in),	dimension(:)	::	x
real(kind=wp),	intent(out),	dimension(:)	::	f

procedure, public :: set_fsal_cache

private subroutine set_fsal_cache(me, t, x, f)

Compute the function and add it to the function evaluation cache, so that it can be reused by the next step (see check_fsal_cache). This is used for f at the end of a step: by the first-same-as-last (FSAL) methods (the last stage of the step), and by the methods whose dense output needs it as one of its extra stages (see dense_extra_stages).

Arguments

Type	Intent	Attributes		Name
class(rk_variable_step_class),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	t
real(kind=wp),	intent(in),	dimension(:)	::	x
real(kind=wp),	intent(out),	dimension(:)	::	f

Source Code

    type,extends(rk_class),abstract,public :: rk_variable_step_class

        !! Main integration class for variable step size Runge-Kutta methods

        private

        type(stepsize_class) :: stepsize_method  !! the method for varying the step size

        real(wp),dimension(:),allocatable :: rtol  !! relative tolerance (`size(n)`)
        real(wp),dimension(:),allocatable :: atol  !! absolute tolerance (`size(n)`)

        integer :: hinit_method = 1 !! if automatically computing the inital step size, which
                                    !! method to use. 1 = `hstart`, 2 = `hinit`.

        integer :: num_rejected_steps = 0 !! number of rejected steps
        real(wp) :: last_accepted_step_size = zero !! the last accepted step size `dt` from the integration
                                                   !! (positive or negative)

        ! for reporting at a fixed step size using the dense output:
        real(wp) :: dense_output_step = zero !! if > 0, the points are reported at this fixed
                                             !! time interval using the dense output (see [[initialize_variable_step]])
        logical :: dense_reporting = .false. !! if the dense output is used for reporting in the current integration
        real(wp) :: dense_output_t0 = zero !! initial time of the current integration
        integer(int64) :: dense_output_k = 0 !! index of the last dense output point reported:
                                             !! `t = dense_output_t0 + k*dense_output_step` (with the
                                             !! sign for the direction of integration)

        ! function evaluation cache (see [[set_fsal_cache]]):
        real(wp),allocatable :: t_saved !! cached `t`
        real(wp),dimension(:),allocatable :: x_saved  !! cached `x`
        real(wp),dimension(:),allocatable :: f_saved !! cached `f`

        contains

        private

        procedure(step_func_variable),deferred :: step !! the step routine for the rk method
        procedure,public :: initialize => initialize_variable_step  !! initialize the class (set n,f, and report)
        procedure,public :: integrate => integrate_variable_step
        procedure,public :: integrate_to_event => integrate_to_event_variable_step
        procedure,public :: info => info_variable_step

        procedure :: hstart  !! for automatically computing the initial step size [this is from DDEABM]
        procedure :: hinit   !! for automatically computing the initial step size [this is from DOP853]
        procedure :: begin_integration => begin_integration_rk_variable_step_class
        procedure :: compute_initial_step
        procedure :: export_step
        procedure :: export_dense_output
        procedure :: order !! returns `p`, the order of the method
        procedure,public :: destroy_fsal_cache
        procedure,public :: check_fsal_cache
        procedure,public :: set_fsal_cache

    end type rk_variable_step_class

rk_variable_step_class Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public, abstract, extends(rk_class) :: rk_variable_step_class

Inherits

Inherited by

Components

Type-Bound Procedures

procedure, public :: destroy

private subroutine destroy(me)

Arguments

procedure, public :: stop => rk_class_stop

private subroutine rk_class_stop(me)

Arguments

procedure, public :: status => rk_class_status

private subroutine rk_class_status(me, istatus, message)

Arguments

procedure, public :: failed

private function failed(me)

Arguments

Return Value logical

procedure, public :: interpolate

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

Arguments

procedure, public :: dense_interval

private subroutine dense_interval(me, available, t0, tf, x0, xf)

Arguments

procedure(properties_func), public, deferred :: properties

pure function properties_func(me) result(p) Prototype

Arguments

Return Value type(rklib_properties)

procedure(step_func_variable), private, deferred :: step

subroutine step_func_variable(me, t, x, h, xf, xerr) Prototype

Arguments

procedure, public :: initialize => initialize_variable_step

private subroutine initialize_variable_step(me, n, f, rtol, atol, stepsize_method, hinit_method, report, g, stop_on_errors, max_number_of_steps, report_rate, solver, dense_output, dense_output_step)

Arguments

procedure, public :: integrate => integrate_variable_step

private subroutine integrate_variable_step(me, t0, x0, h, tf, xf)

Arguments

procedure, public :: integrate_to_event => integrate_to_event_variable_step

private subroutine integrate_to_event_variable_step(me, t0, x0, h, tmax, tol, tf, xf, gf)

Arguments

procedure, public :: info => info_variable_step

private subroutine info_variable_step(me, num_steps, num_rejected_steps, last_accepted_step_size)

Arguments

procedure, private :: hstart

private subroutine hstart(me, a, b, y, yprime, etol, h)

Arguments

procedure, private :: hinit

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

Arguments

Return Value real(kind=wp)

procedure, private :: begin_integration => begin_integration_rk_variable_step_class

private subroutine begin_integration_rk_variable_step_class(me)

Arguments

procedure, private :: compute_initial_step

private function compute_initial_step(me, t0, tf, x0, h0) result(dt)

Arguments

Return Value real(kind=wp)

procedure, private :: export_step

private subroutine export_step(me, t, x, ok)

Arguments

procedure, private :: export_dense_output

private subroutine export_dense_output(me, t_stop, ok)

Arguments

procedure, private :: order

private pure function order(me) result(p)

Arguments

Return Value integer

procedure, public :: destroy_fsal_cache

private subroutine destroy_fsal_cache(me)

Arguments

procedure, public :: check_fsal_cache

private subroutine check_fsal_cache(me, t, x, f)

Arguments

procedure, public :: set_fsal_cache

private subroutine set_fsal_cache(me, t, x, f)