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).
iout=3 for sparse dense output)| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dop853_class), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout) | :: | x |
input: initial value of independent variable.
output: |
||
| real(kind=wp), | intent(inout), | dimension(:) | :: | y |
input: initial values for output: numerical solution at |
|
| real(kind=wp), | intent(in) | :: | xend |
final x-value ( |
||
| real(kind=wp), | intent(in), | dimension(:) | :: | rtol |
relative error tolerance. |
|
| real(kind=wp), | intent(in), | dimension(:) | :: | atol |
absolute error tolerance. |
|
| integer, | intent(in) | :: | iout |
switch for calling the subroutine
|
||
| integer, | intent(out) | :: | idid |
reports on successfulness upon return:
|
subroutine dop853(me,x,y,xend,rtol,atol,iout,idid)
implicit none
class(dop853_class),intent(inout) :: me
real(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(wp),dimension(:),intent(inout) :: y !! *input:* initial values for `y`. [size n]
!!
!! *output:* numerical solution at `x`.
real(wp),intent(in) :: xend !! final x-value (`xend`-`x` may be positive or negative)
real(wp),dimension(:),intent(in) :: rtol !! relative error tolerance. `rtol` and `atol`
!! can be both scalars or else both vectors of length `n`.
real(wp),dimension(:),intent(in) :: 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`:
!!
!! * `iout=0`: subroutine is never called
!! * `iout=1`: subroutine is called after every successful step
!! * `iout=2`: dense output is performed after every successful step
!! * `iout=3`: dense output is performed in steps defined by the user
!! (see `xout` above)
integer,intent(out) :: idid !! reports on successfulness upon return:
!!
!! * `idid=1` computation successful,
!! * `idid=2` comput. successful (interrupted by [[solout]]),
!! * `idid=-1` input is not consistent,
!! * `idid=-2` larger `nmax` is needed,
!! * `idid=-3` step size becomes too small.
!! * `idid=-4` problem is probably stiff (interrupted).
real(wp) :: beta,fac1,fac2,hmax,safe
integer :: i,ieco,iprint,istore,nrdens,nstiff,nmax
logical :: arret
integer :: itol !! switch for `rtol` and `atol`:
!!
!! * `itol=0`: both `rtol` and `atol` are scalars.
!! the code keeps, roughly, the local error of
!! `y(i)` below `rtol*abs(y(i))+atol`.
!!
!! * `itol=1`: both `rtol` and `atol` are vectors.
!! the code keeps the local error of `y(i)` below
!! `rtol(i)*abs(y(i))+atol(i)`.
iprint = me%iprint
me%iout = iout
arret = .false.
!check procedures:
if (.not. associated(me%fcn)) then
if ( iprint/=0 ) &
write (iprint,*) &
'Error in dop853: procedure FCN is not associated.'
idid = -1
return
end if
if (iout/=0 .and. .not. associated(me%solout)) then
if ( iprint/=0 ) &
write (iprint,*) &
'Error in dop853: procedure SOLOUT must be associated if IOUT/=0.'
idid = -1
return
end if
!scalar or vector tolerances:
if (size(rtol)==1 .and. size(atol)==1) then
itol = 0
elseif (size(rtol)==me%n .and. size(atol)==me%n) then
itol = 1
else
if ( iprint/=0 ) &
write (iprint,*) &
'Error in dop853: improper dimensions for rtol and/or atol.'
idid = -1
return
end if
! setting the parameters
me%nfcn = 0
me%nstep = 0
me%naccpt = 0
me%nrejct = 0
nmax = me%nmax
nrdens = me%nrdens !number of dense output components
! nstiff parameter for stiffness detection
if ( me%nstiff<=0 ) then
nstiff = nmax + 10 !no stiffness check
else
nstiff = me%nstiff
end if
if ( nrdens<0 .or. me%nrdens>me%n ) then
if ( iprint/=0 ) write (iprint,*) ' curious input nrdens=' , nrdens
arret = .true.
else
if ( nrdens>0 .and. iout<2 .and. iprint/=0 ) &
write (iprint,*) ' warning: set iout=2 or iout=3 for dense output '
end if
if (size(y)/=me%n) then
if ( iprint/=0 ) &
write (iprint,*) ' error: y must have n elements: size(y)= ',size(y)
arret = .true.
end if
safe = me%safe
fac1 = me%fac1
fac2 = me%fac2
beta = me%beta
if ( me%hmax==0.0_wp ) then
hmax = xend - x
else
hmax = me%hmax
end if
me%h = me%hinitial ! initial step size
! when a fail has occurred, we return with idid=-1
if ( arret ) then
idid = -1
else
! call to core integrator
call me%dp86co(x,y,xend,hmax,me%h,rtol,atol,itol,iprint, &
iout,idid,nmax,nstiff,safe,beta,fac1,fac2, &
me%nfcn,me%nstep,me%naccpt,me%nrejct)
end if
end subroutine dop853