## 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