solve an initial value problem in ordinary differential equations using an adams-bashforth method.
subroutine ddeabm uses the adams-bashforth-moulton
predictor-corrector formulas of orders one through twelve to
integrate a system of neq first order ordinary differential
equations of the form du/dx = df(x,u)
when the vector y(*) of initial values for u(*) at x=t is given.
the subroutine integrates from t to tout. it is easy to continue the
integration to get results at additional tout. this is the interval
mode of operation. it is also easy for the routine to return with
the solution at each intermediate step on the way to tout. this is
the intermediate-output mode of operation.
ddeabm is primarily designed to solve non-stiff and mildly stiff differential equations when derivative evaluations are expensive, high accuracy results are needed or answers at many specific points are required. ddeabm attempts to discover when it is not suitable for the task posed.
quantities which are used as input items are:
neq, t, y(*), tout, info(*), rtol, atol
quantities which may be altered by the code are:
t, y(*), info(1), rtol, atol, idid
the parameters are
neq -- this is the number of (first order) differential
equations to be integrated.
t -- this is a double precision value of the independent
variable.
y(*) -- this double precision array contains the solution
components at t.
tout -- this is a double precision point at which a solution is
desired.
info(*) -- the basic task of the code is to integrate the
differential equations from t to tout and return an
answer at tout. info(*) is an integer array which is used
to communicate exactly how you want this task to be
carried out.
rtol, atol -- these double precision quantities represent
relative and absolute error tolerances which you
provide to indicate how accurately you wish the
solution to be computed. you may choose them to be
both scalars or else both vectors.
idid -- this scalar quantity is an indicator reporting what
the code did. you must monitor this integer variable to
decide what action to take next.
the first call of the code is defined to be the start of each new problem. read through the descriptions of all the following items, provide sufficient storage space for designated arrays, set appropriate variables for the initialization of the problem, and give information about how you want the problem to be solved.
neq -- set it to the number of differential equations.
(neq >= 1)
t -- set it to the initial point of the integration.
you must use a program variable for t because the code
changes its value.
y(*) -- set this vector to the initial values of the neq solution
components at the initial point. you must dimension y at
least neq in your calling program.
tout -- set it to the first point at which a solution
is desired. you can take tout = t, in which case the code
will evaluate the derivative of the solution at t and
return. integration either forward in t (tout > t) or
backward in t (tout < t) is permitted.
the code advances the solution from t to tout using
step sizes which are automatically selected so as to
achieve the desired accuracy. if you wish, the code will
return with the solution and its derivative following
each intermediate step (intermediate-output mode) so that
you can monitor them, but you still must provide tout in
accord with the basic aim of the code.
the first step taken by the code is a critical one
because it must reflect how fast the solution changes near
the initial point. the code automatically selects an
initial step size which is practically always suitable for
the problem. by using the fact that the code will not step
past tout in the first step, you could, if necessary,
restrict the length of the initial step size.
for some problems it may not be permissible to integrate
past a point tstop because a discontinuity occurs there
or the solution or its derivative is not defined beyond
tstop. when you have declared a tstop point (see info(4)),
you have told the code not to integrate past tstop.
in this case any tout beyond tstop is invalid input.
info(*) -- use the info array to give the code more details about
how you want your problem solved. this array should be
dimensioned of length 4. you must respond to all of
the following items which are arranged as questions. the
simplest use of the code corresponds to answering all
questions as yes ,i.e. setting all entries of info to 0.
info(1) -- this parameter enables the code to initialize
itself. you must set it to indicate the start of every
new problem.
**** is this the first call for this problem ...
yes -- set info(1) = 0
no -- not applicable here.
see below for continuation calls. ****
info(2) -- how much accuracy you want of your solution
is specified by the error tolerances rtol and atol.
the simplest use is to take them both to be scalars.
to obtain more flexibility, they can both be vectors.
the code must be told your choice.
**** are both error tolerances rtol, atol scalars ...
yes -- set info(2) = 0
and input scalars for both rtol and atol
no -- set info(2) = 1
and input arrays for both rtol and atol ****
info(3) -- the code integrates from t in the direction
of tout by steps. if you wish, it will return the
computed solution and derivative at the next
intermediate step (the intermediate-output mode) or
tout, whichever comes first. this is a good way to
proceed if you want to see the behavior of the solution.
if you must have solutions at a great many specific
tout points, this code will compute them efficiently.
**** do you want the solution only at
tout (and not at the next intermediate step) ...
yes -- set info(3) = 0
no -- set info(3) = 1 ****
info(4) -- to handle solutions at a great many specific
values tout efficiently, this code may integrate past
tout and interpolate to obtain the result at tout.
sometimes it is not possible to integrate beyond some
point tstop because the equation changes there or it is
not defined past tstop. then you must tell the code
not to go past.
**** can the integration be carried out without any
restrictions on the independent variable t ...
yes -- set info(4)=0
no -- set info(4)=1
and define the stopping point me%tstop ****
rtol, atol -- you must assign relative (rtol) and absolute (atol)
error tolerances to tell the code how accurately you want
the solution to be computed. they must be defined as
program variables because the code may change them. you
have two choices --
both rtol and atol are scalars. (info(2)=0)
both rtol and atol are vectors. (info(2)=1)
in either case all components must be non-negative.
the tolerances are used by the code in a local error test
at each step which requires roughly that
abs(local error) <= rtol*abs(y)+atol
for each vector component.
(more specifically, a euclidean norm is used to measure
the size of vectors, and the error test uses the magnitude
of the solution at the beginning of the step.)
the true (global) error is the difference between the true
solution of the initial value problem and the computed
approximation. practically all present day codes,
including this one, control the local error at each step
and do not even attempt to control the global error
directly. roughly speaking, they produce a solution y(t)
which satisfies the differential equations with a
residual r(t), dy(t)/dt = df(t,y(t)) + r(t),
and, almost always, r(t) is bounded by the error
tolerances. usually, but not always, the true accuracy of
the computed y is comparable to the error tolerances. this
code will usually, but not always, deliver a more accurate
solution if you reduce the tolerances and integrate again.
by comparing two such solutions you can get a fairly
reliable idea of the true error in the solution at the
bigger tolerances.
setting atol=0.0_wp results in a pure relative error test on
that component. setting rtol=0. results in a pure absolute
error test on that component. a mixed test with non-zero
rtol and atol corresponds roughly to a relative error
test when the solution component is much bigger than atol
and to an absolute error test when the solution component
is smaller than the threshold atol.
proper selection of the absolute error control parameters
atol requires you to have some idea of the scale of the
solution components. to acquire this information may mean
that you will have to solve the problem more than once. in
the absence of scale information, you should ask for some
relative accuracy in all the components (by setting rtol
values non-zero) and perhaps impose extremely small
absolute error tolerances to protect against the danger of
a solution component becoming zero.
the code will not attempt to compute a solution at an
accuracy unreasonable for the machine being used. it will
advise you if you ask for too much accuracy and inform
you as to the maximum accuracy it believes possible.
the principal aim of the code is to return a computed solution at tout, although it is also possible to obtain intermediate results along the way. to find out whether the code achieved its goal or if the integration process was interrupted before the task was completed, you must check the idid parameter.
t -- the solution was successfully advanced to the
output value of t.
y(*) -- contains the computed solution approximation at t.
you may also be interested in the approximate derivative
of the solution at t. it is contained in me%ypout.
idid -- reports what the code did
*** task completed ***
reported by positive values of idid
idid = 1 -- a step was successfully taken in the
intermediate-output mode. the code has not
yet reached tout.
idid = 2 -- the integration to tout was successfully
completed (t=tout) by stepping exactly to tout.
idid = 3 -- the integration to tout was successfully
completed (t=tout) by stepping past tout.
y(*) is obtained by interpolation.
*** task interrupted ***
reported by negative values of idid
idid = -1 -- a large amount of work has been expended.
(maxnum steps attempted)
idid = -2 -- the error tolerances are too stringent.
idid = -3 -- the local error test cannot be satisfied
because you specified a zero component in atol
and the corresponding computed solution
component is zero. thus, a pure relative error
test is impossible for this component.
idid = -4 -- the problem appears to be stiff.
idid = -5,-6,-7,..,-32 -- not applicable for this code
but used by other members of depac or possible
future extensions.
*** task terminated ***
reported by the value of idid=-33
idid = -33 -- the code has encountered trouble from which
it cannot recover. a message is printed
explaining the trouble and control is returned
to the calling program. for example, this occurs
when invalid input is detected.
rtol, atol -- these quantities remain unchanged except when
idid = -2. in this case, the error tolerances have been
increased by the code to values which are estimated to be
appropriate for continuing the integration. however, the
reported solution at t was obtained using the input values
of rtol and atol.
this code is organized so that subsequent calls to continue the
integration involve little (if any) additional effort on your
part. you must monitor the idid parameter in order to determine
what to do next.
recalling that the principal task of the code is to integrate
from t to tout (the interval mode), usually all you will need
to do is specify a new tout upon reaching the current tout.
do not alter any quantity not specifically permitted below,
in particular do not alter neq, t, y(*), the internal work variables
or the differential equation in subroutine df. any such alteration
constitutes a new problem and must be treated as such, i.e.
you must start afresh.
you cannot change from vector to scalar error control or vice
versa (info(2)) but you can change the size of the entries of
rtol, atol. increasing a tolerance makes the equation easier
to integrate. decreasing a tolerance will make the equation
harder to integrate and should generally be avoided.
you can switch from the intermediate-output mode to the
interval mode (info(3)) or vice versa at any time.
if it has been necessary to prevent the integration from going
past a point tstop (info(4)), keep in mind that the
code will not integrate to any tout beyond the currently
specified tstop. once tstop has been reached you must change
the value of tstop or set info(4)=0. you may change info(4)
or tstop at any time but you must supply the value of tstop
whenever you set info(4)=1.
the parameter info(1) is used by the code to indicate the
beginning of a new problem and to indicate whether integration
is to be continued. you must input the value info(1)=0
when starting a new problem. you must input the value
info(1)=1 if you wish to continue after an interrupted task.
do not set info(1)=0 on a continuation call unless you
want the code to restart at the current t.
*** following a completed task ***
if
idid = 1, call the code again to continue the integration
another step in the direction of tout.
idid = 2 or 3, define a new tout and call the code again.
tout must be different from t. you cannot change
the direction of integration without restarting.
*** following an interrupted task ***
to show the code that you realize the task was
interrupted and that you want to continue, you
must take appropriate action and reset info(1) = 1
if
idid = -1, the code has attempted maxnum steps.
if you want to continue, set info(1) = 1 and
call the code again. an additional maxnum steps
will be allowed.
idid = -2, the error tolerances rtol, atol have been
increased to values the code estimates appropriate
for continuing. you may want to change them
yourself. if you are sure you want to continue
with relaxed error tolerances, set info(1)=1 and
call the code again.
idid = -3, a solution component is zero and you set the
corresponding component of atol to zero. if you
are sure you want to continue, you must first
alter the error criterion to use positive values
for those components of atol corresponding to zero
solution components, then set info(1)=1 and call
the code again.
idid = -4, the problem appears to be stiff. it is very
inefficient to solve such problems with ddeabm.
the code ddebdf in depac handles this task
efficiently. if you are absolutely sure you want
to continue with ddeabm, set info(1)=1 and call
the code again.
*** following a terminated task ***
if
idid = -33, you cannot continue the solution of this
problem. an attempt to do so will result in your
run being terminated.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ddeabm_class), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | neq |
the number of (first order) differential equations to be integrated. (neq >= 1) |
||
| real(kind=wp), | intent(inout) | :: | t |
value of the independent variable. on input, set it to the initial point of the integration. the code changes its value. on output, the solution was successfully advanced to the output value of t. |
||
| real(kind=wp), | intent(inout), | dimension(neq) | :: | y | ||
| real(kind=wp), | intent(in) | :: | tout | |||
| integer, | intent(inout), | dimension(4) | :: | info | ||
| real(kind=wp), | intent(inout), | dimension(:) | :: | rtol | ||
| real(kind=wp), | intent(inout), | dimension(:) | :: | atol | ||
| integer, | intent(out) | :: | idid |
subroutine ddeabm(me, neq, t, y, tout, info, rtol, atol, idid) implicit none class(ddeabm_class), intent(inout) :: me integer, intent(in) :: neq !! the number of (first order) differential !! equations to be integrated. (neq >= 1) real(wp), intent(inout) :: t !! value of the independent variable. !! on input, set it to the initial point of !! the integration. the code changes its value. !! on output, the solution was successfully !! advanced to the output value of t. real(wp), dimension(neq), intent(inout) :: y real(wp), intent(in) :: tout integer, dimension(4), intent(inout) :: info real(wp), dimension(:), intent(inout) :: rtol real(wp), dimension(:), intent(inout) :: atol integer, intent(out) :: idid character(len=16) :: xern3 if (me%error) return ! if user-triggered error ! check for an apparent infinite loop if (info(1) == 0) me%icount = 0 if (me%icount >= 5) then if (t == me%tprev) then write (xern3, '(1pe15.6)') t call report_error('ddeabm', & 'an apparent infinite loop has been detected. '// & 'you have made repeated calls at t = '//xern3// & ' and the integration has not advanced. check the '// & 'way you have set parameters for the call to the '// & 'code, particularly info(1).', 13, 2) idid = -33 return end if end if !make sure the deriv function was set: if (.not. associated(me%df)) then call report_error('ddeabm', 'the derivative function DF '// & ' has not been associated.', 0, 0) idid = -33 return end if me%tprev = t call me%ddes(neq, t, y, tout, info, rtol, atol, idid, & me%ypout, me%yp, me%yy, me%wt, me%p, me%phi, me%alpha, me%beta, & me%psi, me%v, me%w, me%sig, me%g, me%gi, & me%h, me%eps, me%x, me%xold, me%hold, me%told, me%delsgn, & me%tstop, me%twou, me%fouru, me%start, & me%phase1, me%nornd, me%stiff, me%intout, me%ns, me%kord, & me%kold, me%init, me%ksteps, & me%kle4, me%iquit, me%kprev, me%ivc, me%iv, me%kgi) if (me%error) return ! user-triggered error if (idid /= (-2)) me%icount = me%icount + 1 if (t /= me%tprev) me%icount = 0 end subroutine ddeabm