Integration of ordinary differential equations

The classes and methods defined here allow the integration of a set of simultaneous ordinary differential equations (ODEs) of a real variable. Two approaches are used to allow integrations to be performed. Methods are provided for which a function name can be passed as an argument, along with control information to define the integration. Classes are provided from which a user-defined function can inherit, and methods in the class can be used to perform the integration.

Methods to perform integration

Runge-Kutta integration with accuracy checking

Vector rkaccint(
   // Data in:
   const double x1, // start of interval (initial function values known)
   const double x2, // end of interval (function values to be calculated)
   const double acc, // desired accuracy parameter (e.g. 1e-4)
   const int nsugsteps, // suggested number of steps over interval
   const int maxsteps, // max steps allowed (0 => no limit)
   Vector (*dfdx)(const double, const Vector&), // derivatives routine
   const int verbose, // true for extra diagnostic output
   const double eps, // numerical bandwidth (e.g. 1e-6) used in scaling
      //  and to detect stepsize underflow
   const Vector& f1, // initial function values
   // Data out:
   int& nsteps, // number of integration steps needed
   int& hitmax, // true if reached max iterations
   int& underflow // true if space step became too small
)

Classes to define and integrate functions

ODE integrator class, ODEvector

1. Define a subclass which inherits from this class for the particular ODE to be solved.
2. Define a method in the class to evaluate the gradients given the ordinate and a vector of parameters:

      Vector dfdx(const double x, const Vector& f) const
      

   (This is the most simple form. Some of the integration methods can use more information, as described below.)
3. Set control parameters for the integration. The control parameters (constituents of the ODEvector class) are:

   double acc; // accuracy, e.g. 1.0e-8
   int nsugsteps; // suggested integration steps, e.g. 10
   int maxsteps; // max integration steps; 0 => no limit
   int verbose; // 1 for diagnostic output
   double eps; // bandwidth, e.g. 1.0e-10
   double sfac; // step change safety factor, e.g. 0.9
      

   There is a call to set defaults, the unimaginatively-named

      set_defaults();
      

   The defaults are:

      acc = 1.0e-8;
      nsugsteps = 10;
      maxsteps = 0;
      verbose = 0;
      eps = 1.0e-10;
      sfac = 0.9;
      

   The control parameters are publically-aceessible, and can also be set individually.

4. Call the appropriate integration method:
   • Runge-Kutta integration with accuracy checking. Adapts stepsize to ensure accuracy. Compares step with 2 halfsteps to determine error. Uses 5th order scalings to predict improved stepsizes. Additional error reduction implemented.

     Vector rkaccint(
        // Data in:
        const double x1, // start of interval (initial function values known)
        const double x2, // end of interval (function values to be calculated)
        const Vector& f1, // initial function values
        // Data out:
        int& nsteps, // number of integration steps needed
        int& hitmax, // true if reached max iterations
        int& underflow // true if space step became too small
     ) const
           

   • Runge-Kutta integration with accuracy checking and flexible stop conditions. Adapts stepsize to ensure accuracy. Compares step with 2 halfsteps to determine error. Uses 5th order scalings to predict improved stepsizes. Additional error reduction implemented. Accepts `gone too far' and `stop now' flags from gradient routine to allow greater control e.g. for stopping when dependent values (or a function of them) reach certain values. Integrates as usual, but stops at end of interval or when `stop now' flag set. If `gone too far' is set, will reduce stepsize and try again.

     Vector rksaccint(
        // Data in:
        const double x1, // start of interval (initial function values known)
        const double x2, // end of interval (function values to be calculated)
        const Vector& f1, // initial function values
        // Data out:
        int& nsteps, // number of integration steps needed
        int& hitmax, // true if reached max iterations
        int& underflow, // true if space step became too small
        int& earlystop, // true if stopped on `stop now' condition
        double& xearly // position of early stop
     ) const
           

     Uses a gradient function of the following form:

        Vector dfdx(const double x, const Vector& f, int& stopnow, int& toofar) const
        

     where stopnow is returned as 1 to stop the integration, and toofar is returned as 1 if the integration has proceeded too far.

     A default gradient function of this form is defined, calling the simple gradient function described above and checking the state with a call to a function

     void statechk(const double x, const Vector& f, int& stopnow, int& toofar) const
        

     A default statechk is defined which always returns zero for stopnow and toofar.
   • Runge-Kutta integration with accuracy checking and flexible stop conditions; also saves history of intermediate states from the integration. Adapts stepsize to ensure accuracy. Compares step with 2 halfsteps to determine error. Uses 5th order scalings to predict improved stepsizes. Additional error reduction implemented. Accepts `gone too far' and `stop now' flags from gradient routine to allow greater control e.g. for stopping when dependent values (or a function of them) reach certain values. Integrates as usual, but stops at end of interval or when `stop now' flag set. If `gone too far' is set, will reduce stepsize and try again. Records x,f values at end of each step in history arrays. These are passed to a user analysis routine at the end of each successful integration step, allowing the integration to be stopped early.

     Array<Vector> rkaccinthist(
        // Data in:
        const double x1, // start of interval (initial function values known)
        const double x2, // end of interval (function values to be calculated)
        const Vector& f1, // initial function values
        // Data out:
        int& nsteps, // number of integration steps needed
        int& hitmax, // true if reached max iterations
        int& underflow, // true if space step became too small
        int& earlystop // true if stopped on `stop now' condition
     ) const
        

     Each element of the array of returned vectors is a vector whose first element is the ordinate x at that point and the remaining elements are the function values.

     The flexible stop condition is as above, with the addition of a test based on the state history:

     void histchk(const List<Vector>& fhist, int& stopnow, int& toofar) const
        

     where fhist is the saved history so far.

     A default histchk is defined which always returns zero for stopnow and toofar.