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.
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 )
ODEvector
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.)
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.
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
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) constwhere
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) constA default
statechk
is defined which always returns zero
for stopnow
and toofar
.
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 ) constEach 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) constwhere
fhist
is the saved history so far.
A default histchk
is defined which always returns zero
for stopnow
and toofar
.