DDEABM Name
DDEABM
Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov
UPDATED VERSIONs can be found on my WEB PAGE:
http://cow.physics.wisc.edu/~craigm/idl/idl.html
Purpose
Integrate Ordinary Differential Equation with Adams-Bashford-Moulton
Major Topics
Numerical Analysis.
Calling Sequence
DDEABM, FUNCT, T0, F0, TOUT, [ PRIVATE, FUNCTARGS=, STATE= , ]
[ /INIT, /INTERMEDIATE, TSTOP=, EPSREL=, EPSABS=, ]
[ TGRID=, YGRID=, YPGRID=, NOUTGRID=, NGRID=, NFEV=, ]
[ TIMPULSE=, YIMPULSE=, ]
[ MAX_STEPSIZE=, /CONTROL, ]
[ STATUS=, ERRMSG= ]
Description
DDEABM performs integration of a system of one or more ordinary
differential equations using a Predictor-Corrector technique. An
adaptive Adams-Bashford-Moulton method of variable order between
one and twelve, adaptive stepsize, and error control, is used to
integrate equations of the form:
DF_DT = FUNCT(T, F)
T is the independent variable, F is the (possibly vector) function
value at T, and DF_DT is the derivative of F with respect to T,
evaluated at T. FUNCT is a user function which returns the
derivative of one or more equations.
DDEABM is based on the public domain procedure DDEABM.F written by
L. F. Shampine and M. K. Gordon, and available in the DEPAC package
of solvers within SLATEC library.
DDEABM is used primarily to solve non-stiff and mildly stiff
ordinary differential equations, where evaluation of the user
function is expensive, or high precision is demanded (close to the
machine precision). A Runge-Kutta technique may be more
appropriate if lower precision is acceptable. For stiff problems
neither Adams-Bashford-Moulton nor Runge-Kutta techniques are
appropriate. DDEABM has code which detects the stiffness of the
problem and reports it.
The user can operate DDEABM in three different modes:
* the user requests one output at a specific point (the default),
and maintains the integrator state using the STATE keyword;
* the user requests a grid of points (by passing an array to
TOUT), and optionally maintains the integrator state for
subsequent integrations using the STATE keyword;
* the user requests output at the natural adaptive stepsizes
using the /INTERMEDIATE keyword.
When the user requests explicit output points using TOUT, it is
likely that the output will be interpolated between two natural
stepsize points.
If the user requests a grid of outputs, and the integration fails
before reaching the requested endpoint, then control returns
immediately to the user with the appropriate STATUS code.
The initial conditions are given by F0, when T = T0. The number of
equations is determined by the number of elements in F0.
Integration stops when the independent variable T reaches the final
value of TOUT. If the user wants to continue the integration, it
can be restarted with a new call to DDEABM, and a new value of
TOUT. User Function
The user function FUNCT must be of the following form:
FUNCTION FUNCT, T, F, PRIVATE, ... [ CONTROL=CONTROL ] ...
; ... compute derivatives ...
RETURN, DF_DT
END
The function must accept at least two, but optionally three,
parameters. The first, 'T', is the scalar independent variable
where the derivatives are to be evaluated. The second, 'F', is the
vector of function values. The function must return an array of
same size as 'F'. The third positional parameter, 'PRIVATE', is a
purely optional PRIVATE parameter. FUNCT may accept more
positional parameters, but DDEABM will not use them. Any number of
keyword parameters can be passed using the FUNCTARGS keyword.
The user function FUNCT must not modify either the independent
variable 'T' or the function values 'F'.
PASSING 'CONTROL' MESSAGES TO THE USER FUNCTION
DDEABM may pass control information to the user function, other
than requests for function evaluation. DDEABM will only do this if
the /CONTROL keyword is set when DDEABM was invoked.
When control information has been enabled, the user function *must*
accept a keyword named CONTROL. A message may be passed in this
keyword. If the user function is invoked with the CONTROL keyword
defined, the user function should not evaluate the function, but
rather it must respond to the message and return the appropriate
value.
The CONTROL message will be a structure of the form,
{ MESSAGE: 'name', ... }
where the MESSAGE field indicates a control message. Additional
fields may carry more information, depending on the message.
The following control messages are implemented:
* { MESSAGE: 'INITIALIZE' } - the integrator has been initialized
and the user function may also initialize any state variables,
load large data tables, etc.
RETURN: 0 for success, negative number to indicate failure.
If the user function does not recognize a message, a value of 0
should be returned.
RESTARTING THE INTEGRATOR
When restarting the integrator, the STATE keyword can be used to
save some computation time. Upon return from DDEABM the STATE
keyword will contain a structure which describes the state of the
integrator at the output point. The user need merely pass this
variable back into the next call to continue integration. The
value of T, and the function value F, must not be adjusted between
calls to DDEABM.
If T or F0 are to be adjusted, then the integrator must be
restarted afresh *without* the previous state. This can be
achieved either by reseting STATE to undefined, or by setting the
/INIT keyword to DDEABM.
ERROR CONTROL
Local error control is governed by two keywords, EPSABS and EPSREL.
The former governs the absolute error, while the latter governs the
relative or fractional error. The error test at each iteration is:
ABS(ERROR) LE EPSREL*ABS(F) + EPSABS
A scalar value indicates the same constraint should be applied to
every equation; a vector array indicates constraints should be
applied individually to each component.
Setting EPSABS=0.D0 results in a pure relative error test on that
component. Setting EPSREL=0. results in a pure absolute error test
on that component. A mixed test with non-zero EPSREL and EPSABS
corresponds roughly to a relative error test when the solution
component is much bigger than EPSABS and to an absolute error test
when the solution component is smaller than the threshold EPSABS.
Proper selection of the absolute error control parameters EPSABS
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 EPSREL 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.
HARD LIMIT ON T
If for some reason there is a hard limit on the independent
variable T, then the user should specify TSTOP. For efficiency
DDEABM may try to integrate *beyond* the requested point of
termination, TOUT, and then interpolate backwards to obtain the
function value at the requested point. If this is not possible
because the function because the equation changes, or if there is a
discontinuity, then users should specify a value for TSTOP; the
integrator will not go past this point.
DISCONTINUITIES
If the ODE or solution has discontinuities at known points, these
should be passed to DDEABM in order to aid the solution. The
TIMPULSE and YIMPULSE keyword variables allow the user to identify
the positions of the discontinuities and their amplitudes. As T
crosses TIMPULSE(i) the solution will change from Y to
Y+YIMPULSE(*,i) in a stepwise fashion.
Discontinuities in the function to be integrated can also be
entered in this way. Although DDEABM can adapt the integration
step size to accomodate changes in the user function, it may be
better to identify such discontinuities. In that case TIMPULSE(i)
should still identify the position of discontinuity, and
YIMPULSE(*,i) should be 0.
Currently this functionality is implemented with restarts of the
integrator at the crossing points of the discontinuities.
This technique handles only discontinuities at explicitly known
values of T. If the discontinuity condition is defined in terms of
Y (or Y and T), then the condition is implicit. DDEABM does not
handle that type of condition.
You may list the TIMPULSE points in the TOUT output grid. If an
impulse point appears once in TOUT, the corresponding function
values reported in YGRID and YPGRID will be from *before* crossing
the discontinuity. If the same TIMPULSE point appears *twice* in
TOUT, then the first and second values will correspond to before
and after crossing the discontinuity, respectively.
Inputs
FUNCT - by default, a scalar string containing the name of an IDL
function to be integrated. See above for the formal
definition of FUNCT. (No default).
T0 - scalar number, upon input the starting value of the
independent variable T. Upon output, the final value of T.
Y - vector. Upon input the starting values of the function for T =
T0. Upon output, the final (vector) value of the function.
TOUT - must be at least a scalar, but optionally a vector,
specifies the desired points of output.
If TOUT is a scalar and INTERMEDIATE is not set, then DDEABM
integrates up to TOUT. A scalar value of the function at
TOUT is returned in F.
If TOUT is a scalar and /INTERMEDIATE is set, then DDEABM
computes a grid of function values at the optimal step
points of the solution. The grid of values is returned in
TGRID, YGRID, and YPGRID. The final function value,
evaluated at the last point in TOUT, is returned in F.
If TOUT is a vector, then DDEABM computes a grid of function
values at the requested points. The grid of values is
returned in TGRID, YGRID and YPGRID. The final function
value, evaluated at the last point in TOUT, is returned in
F. If integrating forwards (TOUT greater than T0), TOUT
must be strictly increasing. Generally speaking, TOUT
output points should not repeat, except for discontinuities
as noted above.
It is possible for TOUT to be less than T0, i.e., to
integrate backwards, in which case TOUT must be strictly
decreasing instead.
PRIVATE - any optional variable to be passed on to the function to
be integrated. For functions, PRIVATE is passed as the
second positional parameter. DDEABM does not examine or
alter PRIVATE.
Keyword Parameters
CONTROL - if set, then control messages will be set to the user
function as described above. If not set, then no
control messages will be passed.
EPSABS - a scalar number, the absolute error tolerance requested
in the integral computation. If less than or equal to
zero, then the value is ignored.
Default: 0
EPSREL - a scalar number, the relative (i.e., fractional) error
tolerance in the integral computation. If less than or
equal to zero, then the value is ignored.
Default: 1e-4 for float, 1d-6 for double
ERRMSG - upon return, a descriptive error message.
FUNCTARGS - A structure which contains the parameters to be passed
to the user-supplied function specified by FUNCT via
the _EXTRA mechanism. This is the way you can pass
additional data to your user-supplied function without
using common blocks. By default, no extra parameters
are passed to the user-supplied function.
INIT - if set, indicates that the integrator is to be restarted
afresh.
INTERMEDIATE - if set, indicates that the integrator is to compute
at the natural step points.
MAX_STEPSIZE - a positive scalar value, the maximum integration
step size to take per iteration. The lesser of the
"natural" step size and MAX_STEPSIZE is used. If
MAX_STEPSIZE is not specified, there is no maximum.
NFEV - upon output, the scalar number of function evaluations.
NGRID - if /INTERMEDIATE is set, the requested number of points to
compute before returning. DDEABM uses this value to
allocate storage for TGRID, YGRID, and YPGRID. Note that
DDEABM may not actually calculate this many points. The
user must use NOUTGRID upon return to determine how many
points are valid.
Default: 1
NOUTGRID - upon output, the number of grid points computed. This
may be smaller than the requested number of grid points
(either via NGRID or TOUT) if an error occurs.
STATE - upon input and return, the integrator state. Users should
not modify this structure. If the integrator is to be
restarted afresh, then the /INIT keyword should be set, in
order to clear out the old state information.
STATUS - upon output, the integer status of the integration.
1 - A step was successfully taken in the
intermediate-output mode. The code has not yet
reached TOUT.
2 - The integration to TOUT was successfully completed
(T=TOUT) by stepping exactly to TOUT.
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 STATUS
-1 - A large amount of work has been expended. (500 steps
attempted)
-2 - The error tolerances are too stringent.
-3 - The local error test cannot be satisfied because you
specified a zero component in EPSABS and the
corresponding computed solution component is zero.
Thus, a pure relative error test is impossible for
this component.
-4 - The problem appears to be stiff.
*** Task Terminated ***
-33 - The code has encountered trouble from which it
cannot recover. A error message is printed
explaining the trouble and control is returned to
the calling program. For example, this occurs when
invalid input is detected.
-16 - The user function returned a non-finite
TGRID - upon output, the grid of values of T for which the
integration is provided. Upon return, only values
TGRID(0:NOUTGRID-1) are valid. The remaining values are
undefined.
TIMPULSE - array of values of T where discontinuities occur. The
array should be in ascending order. TIMPULSE must
match YIMPULSE.
TSTOP - a scalar, specifies a hard limit for T, beyond which the
integration must not proceed.
Default: none, i.e., integration is allowed to
"overshoot"
YGRID - upon output, the grid of function values for which the
integration is provided. Upon return, only values
YGRID(*,0:NOUTGRID-1) are valid. The remaining values are
undefined.
YIMPULSE - array of discontinuity offset values, of the form
DBLARR(NEQ,NIMPULSE), where NEQ is the size of Y and
NIMPULSE is the size of TIMPULSE. YIMPULSE(*,I) is the
amount to *add* to Y when T crosses TIMPULSE(I) going
in the positive direction.
YPGRID - upon ouptut, the grid of function derivative values at
the points where the integration is provided. Upon
return, only values YPGRID(*,0:NOUTGRID-1) are valid.
The remaining values are undefined.
Examples
This is a simple example showing how to computes orbits of a
satellite around the earth using Newton's law of gravity. The
earth is assumed to be a central point mass, modeled by the
NEWTON_G function which follows. We assume that the satellite is
orbiting at a radius of 7000 km. The state vector F has six
elements consisting of the position and velocity of the satellite.
POSITION VELOCITY
F = [ X, Y, Z, VX, VY, VZ]
The function NEWTON_G below computes the derivative of F, that is,
VELOCITY ACCELERATION
dF_dt = [ VX, VY, VZ, AX, AY, AZ]
Where the acceleration vector [AX,AY,AZ] is computed using Newton's
laws.
GM = 3.986005d14 ; [MKS] - gravitational constant for earth
a0 = 7000d3 ; [m] - initial radius
v0 = sqrt(GM/a0) ; [m/s] - initial circular velocity
; POSITION VELOCITY
f0 = [a0,0,0, 0,-v0,0] ; initial state vector
t0 = 100d ; [s] Initial time value, meaningless in this case
; Initial output time grid (10000 seconds)
tout = dindgen(10000) + t0
; Integrate equations of motion, starting at T0, and proceeding to
; the maximum time of TOUT. Here the variable GM is passed using
; the PRIVATE mechanism.
f = f0 & t = t0
ddeabm, 'newton_g', t, f, tout, GM, $
tgrid=tgrid, ygrid=ygrid, ypgrid=ypgrid, noutgrid=noutgrid, $
status=status, errmsg=errmsg
Now YGRID(0:2,*) contains the 3D position of the satellite
YGRID(3:5,*) contains the 3D velocity of the satellite
YPGRID(3:5,*) contains the 3D acceleration of the satellite
An alternate way to call DDEABM is to use its natural gridpoints
rather than requesting explicit gridpoints. In that case, we need
to specify the maximum time value we are expecing with TOUT, and
a maximum number of output grid values using NGRID.
f = f0 & t = t0
tout = 10000d ;; Maximum requested time
ddeabm, 'newton_g', t, f, tout, GM, $
ngrid=3000, noutgrid=noutgrid, $
tgrid=tgrid, ygrid=ygrid, ypgrid=ypgrid, noutgrid=noutgrid, $
status=status, errmsg=errmsg
NOUTGRID contains the actual number of grid values returned by
DDEABM. If NOUTGRID is less than NGRID, then the remaining values
are to be ignored.
TGRID = TGRID(0:NOUTGRID-1)
YGRID = YGRID(*,0:NOUTGRID-1)
YPGRID = YPGRID(*,0:NOUTGRID-1)
The user can then plot these values are use them as desired. The
result should be a circular orbit at radius 7000000d meters, with
constant speed given by V0.
; WORK FUNCTION ----
; The acceleration of Newtonian gravity by a central body
; of mass M.
; T - time (not used)
; f - state vector
; f(0:2) = position vector
; f(3:5) = velocity vector
; GM - central body newtonian constant
FUNCTION NEWTON_G, t, f, GM
r = f(0:2) ; Position vector
v = f(3:5) ; Velocity vector
rsq = total(r^2,1) ;; central body distance, squared
rr = sqrt(rsq) ;; central body distance
;; Newtonian gravitational acceleration, three components
a = - GM/rsq * r/rr
;; Assemble final differential vector
df_dt = [v, a]
return, df_dt
end
References
SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE
Solvers.
"Solving Ordinary Differential Equations with ODE, STEP, and INTRP",
by L. F. Shampine and M. K. Gordon, SLA-73-1060.
SLATEC Common Mathematical Library, Version 4.1, July 1993
a comprehensive software library containing over
1400 general purpose mathematical and statistical routines
written in Fortran 77. (http://www.netlib.org/slatec/)
Modification History
Fix bug in TSTOP keyword, 09 May 2002, CM
Fix behavior of KSTEPS, which caused premature termination, 26
May 2002, CM
Fix two errors in the DDEABM_DINTP interpolation step, 04 Jul 2002, CM
Handle case of IMPULSES more correctly, 25 Sep 2002, CM
Handle case when INIT is not set (default to 1); detect
non-finite user function values and error out with STATUS code
-16; promote integer values to LONG integers; some internal
function renaming, 28 Jun 2003, CM
Fixed bug in handling of DOIMPULSE and INTERMEDIATE, 08 Mar 2004, CM
Corrected interface error in usage of NGRID. Now NGRID is
actually the number of INTERMEDIATE points to compute (and is
input only). NOUTGRID is a new keyword, which provides the
number of output grid points upon return. 08 Mar 2004, CM
Early termination is possible for INTERMEDIATE case. Handle it
properly , 08 Mar 2004, CM
Fix a bug in the handling of INIT (strangely the internal
code keeps two different INIT variables!); this really only
had an effect when continuing a previous integration; handle
impulses properly when integrating in the negative direction;
document the TIMPULSE/YIMPULSE keyword parameters; some other
small code cleanups; 16 Jul 2008, CM
Handle the case when TOUT EQ TIMPULSE, 05 Sep 2008, CM
Further work on TOUT EQ TIMPULSE, also allowing reporting of
function values on either side of a discontinuity, 07 Sep 2008, CM
Add the MAX_STEPSIZE keyword, 01 Oct 2008, CM
Make sure new impulse checks work when integrating in reverse
direction, 09 Oct 2008, CM
New interface requirement: user function must be able to handle
control messages from DDEABM; first one is INITIALIZE,
20 Oct 2008, CM
Change the control message interface so that it is
backward-compatible; the user must now set the /CONTROL keyword
to enable control messages; they are passed to the user
function via the CONTROL keyword, 08 Nov 2008, CM
Update the documentation; the largest change is the inclusion of
a new example, 16 Jan 2010, CM