>  Docs Center  >  Libraries  >  Markwardt  >  DDEABM






  Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
  UPDATED VERSIONs can be found on my WEB PAGE:


  Integrate Ordinary Differential Equation with Adams-Bashford-Moulton

Major Topics

  Numerical Analysis.

Calling Sequence

    [ STATUS=, ERRMSG= ]


  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

User Function

  The user function FUNCT must be of the following form:
    ; ... compute derivatives ...
  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'.
  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
  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.
  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.
  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:
  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.
  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.
  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.


  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
  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
                      *** Task Interrupted ***
                Reported by negative values of STATUS
          -1 - A large amount of work has been expended. (500 steps
          -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
  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
  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
  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.


  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.
    F = [ X, Y, Z, VX, VY, VZ]
  The function NEWTON_G below computes the derivative of F, that is,
    dF_dt = [ VX, VY, VZ, AX, AY, AZ]
  Where the acceleration vector [AX,AY,AZ] is computed using Newton's
  GM = 3.986005d14 ; [MKS] - gravitational constant for earth
  a0 = 7000d3 ; [m] - initial radius
  v0 = sqrt(GM/a0) ; [m/s] - initial circular 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.
  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.
  ; 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
    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


  SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE
  "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

© 2020 Harris Geospatial Solutions, Inc. |  Legal
My Account    |    Store    |    Contact Us