dop853_class Derived Type

type, public :: dop853_class


Contents

Variables

Type-Bound Procedures

Source Code

dop853_class

Components

Type Visibility Attributes Name Initial
integer, private :: n = 0

the dimension of the system
integer, private :: nfcn = 0

number of function evaluations
integer, private :: nstep = 0

number of computed steps
integer, private :: naccpt = 0

number of accepted steps
integer, private :: nrejct = 0

number of rejected steps (due to error test), (step rejections in the first step are not counted)
integer, private :: nrdens = 0

number of components, for which dense output is required. for 0 < nrdens < n the components (for which dense output is required) have to be specified in icomp(1),...,icomp(nrdens).
real(kind=wp), private :: h = 0.0_wp

predicted step size of the last accepted step
integer, private :: iprint = output_unit

switch for printing error messages if iprint==0 no messages are being printed if iprint/=0 messages are printed with write (iprint,*) ...
integer, private :: nmax = 100000

the maximal number of allowed steps.
integer, private :: nstiff = 1000

test for stiffness is activated after step number j*nstiff (j integer), provided nstiff>0. for negative nstiff the stiffness test is never activated.
real(kind=wp), private :: hinitial = 0.0_wp

initial step size, for hinitial=0 an initial guess is computed with help of the function hinit.
real(kind=wp), private :: hmax = 0.0_wp

maximal step size, defaults to xend-x if hmax=0.
real(kind=wp), private :: safe = 0.9_wp

safety factor in step size prediction
real(kind=wp), private :: fac1 = 0.333_wp

parameter for step size selection. the new step size is chosen subject to the restriction fac1 <= hnew/hold <= fac2
real(kind=wp), private :: fac2 = 6.0_wp

parameter for step size selection. the new step size is chosen subject to the restriction fac1 <= hnew/hold <= fac2
real(kind=wp), private :: beta = 0.0_wp

is the beta for stabilized step size control (see section iv.2). positive values of beta ( <= 0.04 ) make the step size control more stable.
integer, private, dimension(:), allocatable :: icomp

dimension(nrdens) the components for which dense output is required
real(kind=wp), private, dimension(:), allocatable :: cont

dimension(8*nrdens)
integer, private :: iout = 0

copy of iout input to dop853
real(kind=wp), private :: xold = 0.0_wp
real(kind=wp), private :: hout = 0.0_wp
procedure(deriv_func), private, pointer :: fcn => null()

subroutine computing the value of f(x,y)
procedure(solout_func), private, pointer :: solout => null()

subroutine providing the numerical solution during integration. if iout>=1, it is called during integration.

Type-Bound Procedures

procedure, public :: initialize => set_parameters

initialization routine.

  • private subroutine set_parameters(me, n, fcn, solout, iprint, nstiff, nmax, hinitial, hmax, safe, fac1, fac2, beta, icomp, status_ok)

    Set the optional inputs for dop853.

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_class), intent(inout) :: me
    integer, intent(in) :: n

    the dimension of the system (size of and vectors)
    procedure(deriv_func) :: fcn

    subroutine computing the value of

    procedure(solout_func), optional :: solout

    subroutine providing the numerical solution during integration. if iout>=1, it is called during integration. supply a dummy subroutine if iout=0.
    integer, intent(in), optional :: iprint

    switch for printing error messages if iprint==0 no messages are being printed if iprint/=0 messages are printed with write (iprint,*) ...
    integer, intent(in), optional :: nstiff

    test for stiffness is activated after step number j*nstiff (j integer), provided nstiff>0. for negative nstiff the stiffness test is never activated.
    integer, intent(in), optional :: nmax

    the maximal number of allowed steps.
    real(kind=wp), intent(in), optional :: hinitial

    initial step size, for hinitial=0 an initial guess is computed with help of the function hinit.
    real(kind=wp), intent(in), optional :: hmax

    maximal step size, defaults to xend-x if hmax=0.
    real(kind=wp), intent(in), optional :: safe

    safety factor in step size prediction
    real(kind=wp), intent(in), optional :: fac1

    parameter for step size selection. the new step size is chosen subject to the restriction fac1 <= hnew/hold <= fac2
    real(kind=wp), intent(in), optional :: fac2

    parameter for step size selection. the new step size is chosen subject to the restriction fac1 <= hnew/hold <= fac2
    real(kind=wp), intent(in), optional :: beta

    is the beta for stabilized step size control (see section iv.2). positive values of beta ( <= 0.04 ) make the step size control more stable.
    integer, intent(in), optional, dimension(:) :: icomp

    the components for which dense output is required (size from 0 to n).
    logical, intent(out) :: status_ok

    will be false for invalid inputs.

procedure, public :: integrate => dop853

main integration routine.

  • private subroutine dop853(me, x, y, xend, rtol, atol, iout, idid)

    Numerical solution of a system of first order ordinary differential equations . This is an explicit Runge-Kutta method of order 8(5,3) due to Dormand & Prince (with stepsize control and dense output).

    Authors

    • E. Hairer and G. Wanner Universite de Geneve, Dept. De Mathematiques ch-1211 geneve 24, switzerland e-mail: ernst.hairer@unige.ch gerhard.wanner@unige.ch
    • Version of October 11, 2009 (new option iout=3 for sparse dense output)
    • Jacob Williams, Dec 2015: significant refactoring into modern Fortran.

    Reference

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_class), intent(inout) :: me
    real(kind=wp), intent(inout) :: x

    input: initial value of independent variable. output: x for which the solution has been computed (after successful return x=xend).
    real(kind=wp), intent(inout), dimension(:) :: y

    input: initial values for y. [size n]

    Read more…
    real(kind=wp), intent(in) :: xend

    final x-value (xend-x may be positive or negative)
    real(kind=wp), intent(in), dimension(:) :: rtol

    relative error tolerance. rtol and atol can be both scalars or else both vectors of length n.
    real(kind=wp), intent(in), dimension(:) :: atol

    absolute error tolerance. rtol and atol can be both scalars or else both vectors of length n. atol should be strictly positive (possibly very small)
    integer, intent(in) :: iout

    switch for calling the subroutine solout:

    Read more…
    integer, intent(out) :: idid

    reports on successfulness upon return:

    Read more…

procedure, public :: destroy => destroy_dop853

destructor.

procedure, public :: info => get_dop853_info

to get info after a run.

  • private subroutine get_dop853_info(me, n, nfcn, nstep, naccpt, nrejct, h, iout)

    Get info from a dop853_class.

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_class), intent(in) :: me
    integer, intent(out), optional :: n

    dimension of the system
    integer, intent(out), optional :: nfcn

    number of function evaluations
    integer, intent(out), optional :: nstep

    number of computed steps
    integer, intent(out), optional :: naccpt

    number of accepted steps
    integer, intent(out), optional :: nrejct

    number of rejected steps (due to error test), (step rejections in the first step are not counted)
    real(kind=wp), intent(out), optional :: h

    predicted step size of the last accepted step
    integer, intent(out), optional :: iout

    iout flag passed into dop853, used to specify how solout is called during integration.

procedure, private :: dp86co

  • private subroutine dp86co(me, x, y, xend, hmax, h, rtol, atol, itol, iprint, iout, idid, nmax, nstiff, safe, beta, fac1, fac2, nfcn, nstep, naccpt, nrejct)

    Core integrator for dop853. parameters same as in dop853 with workspace added.

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_class), intent(inout) :: me
    real(kind=wp), intent(inout) :: x
    real(kind=wp), intent(inout), dimension(:) :: y
    real(kind=wp), intent(in) :: xend
    real(kind=wp), intent(inout) :: hmax
    real(kind=wp), intent(inout) :: h
    real(kind=wp), intent(in), dimension(:) :: rtol
    real(kind=wp), intent(in), dimension(:) :: atol
    integer, intent(in) :: itol
    integer, intent(in) :: iprint
    integer, intent(in) :: iout
    integer, intent(out) :: idid
    integer, intent(in) :: nmax
    integer, intent(in) :: nstiff
    real(kind=wp), intent(in) :: safe
    real(kind=wp), intent(in) :: beta
    real(kind=wp), intent(in) :: fac1
    real(kind=wp), intent(in) :: fac2
    integer, intent(inout) :: nfcn
    integer, intent(inout) :: nstep
    integer, intent(inout) :: naccpt
    integer, intent(inout) :: nrejct

procedure, private :: hinit

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

    computation of an initial step size guess

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_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
    real(kind=wp), intent(in), dimension(:) :: f0

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

    Return Value real(kind=wp)

procedure, public :: contd8

can be called in user's solout_func for dense output.

  • private function contd8(me, ii, x) result(y)

    this function can be used for continuous output in connection with the output-subroutine for dop853. it provides an approximation to the ii-th component of the solution at x.

    Arguments

    Type IntentOptional Attributes Name
    class(dop853_class), intent(in) :: me
    integer, intent(in) :: ii
    real(kind=wp), intent(in) :: x

    Return Value real(kind=wp)

Source Code

    type,public :: dop853_class

        private

        !internal variables:
        integer :: n      = 0   !! the dimension of the system
        integer :: nfcn   = 0   !! number of function evaluations
        integer :: nstep  = 0   !! number of computed steps
        integer :: naccpt = 0   !! number of accepted steps
        integer :: nrejct = 0   !! number of rejected steps (due to error test),
                                !! (step rejections in the first step are not counted)
        integer :: nrdens = 0   !! number of components, for which dense output
                                !! is required. for `0 < nrdens < n` the components
                                !! (for which dense output is required) have to be
                                !! specified in `icomp(1),...,icomp(nrdens)`.
        real(wp) :: h = 0.0_wp  !! predicted step size of the last accepted step

        !input paramters:
        !  these parameters allow
        !  to adapt the code to the problem and to the needs of
        !  the user. set them on class initialization.

        integer :: iprint = output_unit   !! switch for printing error messages
                                          !! if `iprint==0` no messages are being printed
                                          !! if `iprint/=0` messages are printed with
                                          !! `write (iprint,*)` ...

        integer :: nmax = 100000       !! the maximal number of allowed steps.
        integer :: nstiff = 1000       !! test for stiffness is activated after step number
                                       !! `j*nstiff` (`j` integer), provided `nstiff>0`.
                                       !! for negative `nstiff` the stiffness test is
                                       !! never activated.
        real(wp) :: hinitial = 0.0_wp  !! initial step size, for `hinitial=0` an initial guess
                                       !! is computed with help of the function [[hinit]].
        real(wp) :: hmax = 0.0_wp      !! maximal step size, defaults to `xend-x` if `hmax=0`.
        real(wp) :: safe = 0.9_wp      !! safety factor in step size prediction
        real(wp) :: fac1 = 0.333_wp    !! parameter for step size selection.
                                       !! the new step size is chosen subject to the restriction
                                       !! `fac1 <= hnew/hold <= fac2`
        real(wp) :: fac2 = 6.0_wp      !! parameter for step size selection.
                                       !! the new step size is chosen subject to the restriction
                                       !! `fac1 <= hnew/hold <= fac2`
        real(wp) :: beta = 0.0_wp      !! is the `beta` for stabilized step size control
                                       !! (see section iv.2). positive values of beta ( <= 0.04 )
                                       !! make the step size control more stable.

        integer,dimension(:),allocatable  :: icomp      !! `dimension(nrdens)`
                                                        !! the components for which dense output is required
        real(wp),dimension(:),allocatable :: cont       !! `dimension(8*nrdens)`

        integer :: iout = 0 !! copy of `iout` input to [[dop853]]

        !formerly in the condo8 common block:
        real(wp) :: xold = 0.0_wp
        real(wp) :: hout = 0.0_wp

        !user-defined procedures:
        procedure(deriv_func),pointer  :: fcn    => null() !! subroutine computing the value of `f(x,y)`
        procedure(solout_func),pointer :: solout => null() !! subroutine providing the
                                                           !! numerical solution during integration.
                                                           !! if `iout>=1`, it is called during integration.

        contains

        private

        procedure,public :: initialize => set_parameters    !! initialization routine.
        procedure,public :: integrate  => dop853            !! main integration routine.
        procedure,public :: destroy    => destroy_dop853    !! destructor.
        procedure,public :: info       => get_dop853_info   !! to get info after a run.

        procedure :: dp86co
        procedure :: hinit
        procedure,public :: contd8  !! can be called in user's [[solout_func]] for dense output.

    end type dop853_class