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.

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:

   int verbose; // 1 or higher for diagnostic output
   double acc; // accuracy, e.g. 1.0e-8
   int nsugsteps; // suggested integration steps, e.g. 10
   int maxsteps; // max integration steps; 0 => no limit
   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:

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

   and methods Read_parameters(istream &) and Write_parameters(ostream &) to read and write the parameters (in the same order) to/from an ASCII stream.

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

4. Call the appropriate integration method:
   • Runge-Kutta integration, fixed step size.

     Vector rkint(const double x1, const double x2, const Vector& f1) const
           

   • Runge-Kutta integration, fixed step size, history.

     Array<Vector> rkinthist(const double x1, const double x2, const Vector& f1) const
           

   • 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
           

   • Previous function with slicker interface:

     Vector rkaccint(
        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
     ) 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.

   • Previous function with slicker interfsce:

     Array<Vector> rkaccinthist(
        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
     ) const