Integration of ordinary differential equations with augmented parameter array

These classes and methods are similar to the un-augmented versions, but are more convenient for applications where the end-point of integration is uncertain a priori. The augmented parameter (and function) vector includes the independent variable of integration as its first element. Thus, given a set of functions f_i(x) to integrate simultaneously, the augmented set g_i={x,f_i-1. In the derivative function, the first element returned is dx/dx=1.

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, ODEaugvector

Vector dfdx(const Vector& f) const
{
   Vector df(f.n);
   df[0] = 1.0; // =dx/dx
   // insert calculations for other components
   return df;
}  

For integration functions with an adaptive stepsize, the significance scale of each parameter can be set by including a variant of the following function:

virtual Vector sigscale(const Vector& f) const
{
   Vector sig(f.n);
   sig[0] = 0.0; // ordinate: no error test necessary,
                 // otherwise can be used to set maximum space step
   for (int i=1;i

For integration functions with flexible stop conditions,
parameter dx0 is the trial value of ordinate step.
Subclasses must define a status function to report whether
integration should proceed (0), has terminated (1), or has gone too far (-1):
int status(const Vector& f) const { return 0; }

For integration methods storing the history and using flexible stop
conditions, the status test can be made on the entire history so far:
int status(const List<Vector>& f) const { return 0; }

- this function defaults to a call to the single-value status using the
most recent value calculated.

To integrate the set of functions, the class is used as follows:

Define a subclass which inherits from this class for the particular
   ODE to be solved.
Define the gradient function and the status-testing function(s), if needed.
Define the significance function, if needed.
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
double dx0; // trial integration step, e.g. 0.1
   
   

   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;
   dx0 = 0.1;
   
   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.
Call the appropriate integration method
   (where f1 is the vector of initial function values
   and Dx the integration interval):
   

   
Runge-Kutta integration, fixed step size.
      
Vector rkint(const Vector& f1, const double Dx) const
      
   
Runge-Kutta integration, fixed step size, history.
      
Array<Vector> rkinthist(const Vector& f1, const double Dx) const
      
   
Runge-Kutta integration, fixed step size, stop from status function.
      
Vector rkint(const Vector& f1) const
      
   
Runge-Kutta integration, fixed step size, stop from status function, history.
      
Array<Vector> rkinthist(const Vector& f1) const
      
   
Runge-Kutta integration, adaptive step size.
      
Vector rkaccint(const Vector& f1, const double Dx) const
      
   
Runge-Kutta integration, adaptive step size, history.
      
Array<Vector> rkaccinthist(const Vector& f1, const double Dx) const
      
   
Runge-Kutta integration, adaptive step size, stop from status function.
      
Vector rkaccint(const Vector& f1) const
      
   
Runge-Kutta integration, adaptive step size, stop from status function, history.
      
Array<Vector> rkaccinthist(const Vector& f1) const