A C++ Integrator Class
C++ classes for solving differential equations
Darrel J. Conway
Darrel, who holds a PhD in physics, has worked in educational and industrial settings building numerical models of gravitating systems ranging from black holes to Earth-orbiting spacecraft. He can be contacted at 71203.1415@compuserve.com.
When modeling the behavior of systems, scientists and engineers must frequently solve differential equations. In this article, I'll present a C++ implementation of a class hierarchy that facilitates this process by allowing for rapid incorporation of new integration methods as they become available. As an example of this hierarchy, I'll focus on adaptive stepsize Runge-Kutta integrators. The resulting code contains two integrators: a fifth-order integrator with fourth-order error control (discussed in Numerical Recipes in C, by William Press, et al.), which I refer to as a "Runge-Kutta 4(5) integrator," and the "Runge-Kutta 7(8) integrator" derived by J.H. Verner. (Verner presents several different sets of coefficients for Runge-Kutta integrators, along with a description of the methodology used to produce them.) As you'll see, incorporation of the second integrator into the class structure is a simple matter of entering the number of stages and the coefficients for the integration algorithm.
The general problem I'll address involves a system of i variables, r, in some initial configuration that can be represented by a set of differential equations of the form in Example 1. In this case, t is the independent parameter for the system, and r(i) is the ith variable of the system. I'll solve for the values of the variables r at some other value of t. To illustrate how the integrator classes can be used to solve this problem, I'll set up and integrate the Newtonian three-body problem (see the accompanying text box).
I designed the integration-class structure to be reusable without modification by implementing the following features:
- The integrators do not contain details about the system being integrated. These details are passed to the integrator after an instance of the class is created. Among other things, this means that the integrator classes do not know the dimensions of the system being integrated until an integrator object is instantiated.
- The differential equations that provide derivative information to the integrators are specified outside of the integrator classes.
- The integrator-class structure is flexible. For instance, if someone derives a new integration implementation for one of the techniques implemented in the class hierarchy, the new derivation can be incorporated into the class structure with minimal coding.
- The integrators are reusable for many different sets of differential equations.
As Figure 1 illustrates, the integration algorithms can be divided into three levels. All integration schemes share several features, including information about the derivatives (represented generically by Example 1), the initial conditions at the start of the integration, and a function to take the desired integration step. These features are defined in the Integrator base class.
The integrators need access to information about the system being integrated. This information is provided through the Derivative class and its children; see Figure 2. Each type of integration scheme is derived from the Integrator class; Figure 1 shows two such classes: the AdaptiveRungeKutta class and PredictorCorrector class. Each class contains algorithm-specific details of the integrator; for instance, the AdaptiveRungeKutta class contains pointers to the arrays of coefficients needed for the Runge-Kutta algorithms, along with functions used to evaluate the error incurred during a step and the implementation of the integration stepping routine.
Specific integrators are derived from these intermediate classes. A typical implementation at this level contains the coefficients and descriptions needed to create an instance of the integrator. All of the code required to step the integrator, including error control, is encapsulated in the intermediate class. This formulation allows you to rapidly build multiple, yet different, integrators of a given type.
The code presented here, although written and tested using Borland C++ 4.5, is "vanilla" C++, which should port easily to other development systems and platforms. However, I haven't defined the functions for the copy constructor or the assignment operator for the classes, so you'll have to write them if you need them.
Listing One is the class definition for the Integrator class and Listing Two presents the Integrator base class. The Integrator class can be used as the skeleton for any integrators you implement. It provides pointers to the memory used by the initial state of the variables and the memory location for the integrated variables, along with functions to set these pointers. The Integrator class tracks the number of variables being integrated in the dimension variable, which is set by the constructor for the class. Finally, the Integrator class contains a pointer to Derivative, a base class that provides the framework needed to obtain the derivative information in Example 1.
Integrator is a pure virtual class. Each integrator derived from it needs the Step() function to initiate an integration. This function takes one parameter: a floating-point value for the size of the integration step. Each integration scheme performs different mechanizations to calculate the integration step. Step() is included in the base class to make the interface to the integration consistent across the integrators.
The Integrator class cannot be used without corresponding derivative information. This can be problematic if the resulting Integrator is specific to the set of equations contained in the class. Consequently, I decided to create the Derivative class. Like Integrator, the Derivative class is a skeleton used to make the interface to the derivative information consistent; see Listing Three. Derivative does not contain any data, and (aside from the empty constructor and destructor) consists of a single pure virtual function interface, GetDerivatives(), that takes as its input the state of the independent variable and a pointer to the array of dependent variables whose derivatives need to be calculated. The array of dependent variables is replaced by the corresponding derivatives in GetDerivatives(). Thus, classes derived from the Derivative class write the calculated information to a data array owned by Integrator. This violates the encapsulation of the integrator's internal data, but I decided to implement the integrator this way to improve performance.
The connection between the derivative classes and the integrators makes them logical members of a larger entity: the class that represents the physical or mathematical system being modeled. Once created, the tight coupling between the integrator and the derivatives is hidden from the rest of the program. (This is the type of entity which I'll implement for the three-body system.)
Listings Four and Five present the ThreeBody class, a sample method for the implementation of derivative information. ThreeBody is derived from the Derivative class and inherits GetDerivatives(). The code implements the derivatives needed for the three-body problem.
You still need to specify where the starting and ending data are stored at each integration step. These data are set using the initState and finalState pointers; pointers are set using InitState() and FinalState(). The memory needed for these data is controlled outside of the Integrator class hierarchy.
Once you've set up the derivative information, the dimension of the integration, and the memory locations for the integration data, the external interfaces for the integrator are set. Next, you need to create a class that captures the essence of an integration technique, but leaves out the implementation details. I've implemented the adaptive stepsize Runge-Kutta integrators in the AdaptiveRungeKutta class to illustrate this middle layer of the class structure (see Listings Six and Seven).
Runge-Kutta integrators work by taking a series of substeps (or stages) from the initial point and evaluating the derivatives of the variables at these stages. The accumulated data calculated across the stages is used to generate the final integrated data. For a total step of size h, you can write the ith stage ki(n) for the nth variable y(n) in the integration scheme as in Example 2.
The arrays ai and bij contain constants specific to a particular implementation of the Runge-Kutta algorithm. Example 2 describes how to generate the ith stage when you have the preceding i-1 stages. The first stage is generated using Example 3. After all of the stages have been generated, the integration step is given by Example 4, where s is the number of stages in the implementation and ci is a third set of constant coefficients. Example 2 through Example 4 are implemented in Listing Seven (see the Step() function).
You can control the accuracy of a Runge-Kutta implementation by performing two independent integrations of the starting state across the same interval. If you then subtract the results of these integrations, you will get a local estimate of the error incurred by the numerical integration. This estimate can be used to adapt the stepsize for the integration, and thus to control the error at each step. Numerical Recipes provides a good summary of this form of stepsize control, so I won't go into details.
The error estimate 
(n) for the nth variable can be calculated from two Runge-Kutta steps of order m and m-1. If the implementation is performed appropriately, the mth order step and the (m-1)th order step will use identical Runge-Kutta stages. The error estimate will be simple to calculate from the already determined values of the stages using Example 5, where ci* is the coefficient needed for the (m-1)th order Runge-Kutta step. GetError() uses this equation to determine the largest component of and returns the value of this component; see Listing Seven.
Suppose you want to control error to some local accuracy 
. You perform an integration across a step of size h and find that the achieved accuracy was 
. If >, you need to decrease the requested step and repeat the integration from the same initial state in order to preserve the desired accuracy of the integration. On the other hand, if <, the step was good. You may want to increase the stepsize the next time you take an integration step so that you do not accumulate error from the addition of many small steps. (You may also want to increase the stepsize to make the integration proceed more rapidly!) In either case, you need a way to adapt the step. This adaptation is performed using Example 6, where 
is a "safety" factor, used to keep the new step hnew from growing too large. A typical value for is 0.9. NewStep() (Listing Seven) performs this calculation and returns the new step.
One wrinkle to this algorithm that I've included in the implementation is the ability to take fixed-size Runge-Kutta steps with full error control. I needed a way to perform accurate integration and still produce integration steps that are evenly spaced in time so that simulations would run correctly. A typical implementation of a Runge-Kutta integrator allows the steps to adjust freely, producing integrated steps that are not evenly distributed in time. If you are interested in simulating gravitating systems, such a procedure will produce large steps when the gravitating bodies are far apart and small steps when they are close together. When the resulting points are plotted, the bodies appear to move rapidly when they are far apart and slowly when they are close together. In reality, the bodies move most rapidly when they are closest together. More raw integrated points are produced in this region because the forces are stronger, so the integrator needs to take smaller steps to preserve the accuracy of the integration. The AdaptiveRungeKutta class uses a flag named fixedStep to determine if the integrator should take fixed steps with full error control. When fixedStep is true, Step() will loop through the integration, taking multiple integration steps if necessary, and return results when the requested total step has been performed. If fixedStep is false, Step() will return the "raw" integration step taken by the integrator.
The AdaptiveRungeKutta class is pure virtual because the coefficients ai, bij, and cj are never specified at this level of the class hierarchy. SetCoefficients() is provided for this function but not implemented at this level. Specific implementations of the adaptive Runge-Kutta integrators are left to the third layer of the class hierarchy.
Each implementation of a Runge-Kutta integrator must set the powers used to adjust the stepsize (see Example 6), the value of the safety factor , and the coefficients used to perform the integration. The file rk45.cpp (available electronically; see "Availability," page 3) illustrates this procedure for the RungeKutta45 class. The constructor for the class is called with an integer specifying the dimension of the problem being integrated. The Runge-Kutta 4(5) implementation is a six-stage integrator, so this parameter is passed along with the dimension to the AdaptiveRungeKutta constructor. The RungeKutta45 constructor sets the exponents and the safety factor, and then calls SetCoefficients() to set up the coefficients ai, bij, and cj. The coefficients used by this implementation were derived by Cash and Karp; see Table 1. Most of the work needed to implement the Runge-Kutta integrator is contained in the AdaptiveRungeKutta class. The implementation-specific classes function primarily to set up the implementation details. rk45.cpp is an example of this technique. For the RungeKutta45 class, the class constructor sets the tables of coefficients needed to perform the integration.
The program starsys.cpp (available electronically) implements the Trinary class used to collect the Derivative and Integrator classes together with the arrays of initial- and final-state data into one class for integration. The Trinary class is the "entity" discussed earlier. The driving program creates an instance of this class and calls TakeAStep() with the desired stepsize in order to tell the Trinary system to perform an integration step.
The Trinary constructor initializes the initial state array with data appropriate to a system consisting of three equal-mass bodies located equally distant from each other, moving tangentially to the circle that passes through the masses. This configuration is the starting point for one of the few known analytic solutions to the three-body problem, so it is a good test case for the integrators. The masses should move in such a way as to maintain the symmetry of the initial conditions: Lines connecting the masses should always form an equilateral triangle.
The symmetric problem is a good test case for the integrators, but it can be rather boring. The Trinary class contains Star1(), Star2(), and Star3() to change the initial conditions of the problem. Each function takes a pointer to arrays of positions and velocities, along with a variable for the mass, so that you can experiment with different configurations of the system.
stars.cpp (available electronically) is a simple driver for the Trinary system. As it stands, this driver integrates the default configuration for about an orbit and a half. Figure 3 shows the results of such an integration. The conditions needed for a more-complicated three-body system are included in the code in stars.cpp, but are commented out. You can generate the data for this system by removing the comment tokens.
The file rk78.cpp incorporates the 7(8) integrator into the class structure. Since the lowest level of the integrator hierarchy consists of the implementation details, this class contains the same function calls as the RungeKutta45 class. The Trinary class can use this integrator by changing the preprocessor definition of RKTYPE in starsys3.h (available electronically) from RungeKutta45 to RungeKutta78. Table 2 shows the results of integration of the test problem through 750,000 simulated seconds for both of these integrators.
The class structure here can be adapted to many different types of differential equations. I have applied it to two- and three-body gravitating systems, and to the Lorenz equations found in Chaos theory. Any continuous, differentiable system that can be modeled using equations of the form given in Example 1 can be solved using the integrators presented here. Finally, the structure used to separate the Runge-Kutta algorithm from the implementation-specific details can be extended to other types of integrators.
Cash, R. and A.H. Karp. "A Variable Order Runge-Kutta Method for Initial Value Problems with Rapidly Varying Right-Hand Sides." ACM Transactions on Mathematical Software, 16 (1990), 201.
Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C, Second Edition. Cambridge, U.K.: Cambridge University Press, 1992.
Verner, J.H. "Explicit Runge-Kutta Methods with Estimates of the Local Truncation Error." SIAM Journal of Numerical Analysis, 15 (1978).
The Three-Body Problem
The three-body problem can be stated as follows: Given three gravitating bodies in an initial configuration in space (for instance, the configuration in Figure 4) with some specified initial velocities, find how these bodies move over time.
This problem is sufficiently complicated that no analytic solution exists for the general problem. Special cases of initial conditions have solutions that provide a good framework for testing integrators. For the three-body problem in this article, the parameter t in Example 1 is time, and the integrator is used to make the system evolve in time. The variables for the problem are the positions x, y, and z and velocities vx, vy, and vz of these bodies, represented in a Cartesian system. The state of each body can be represented by the current value of t and six components representing the position and velocity in space. Example 7 presents the differential equations for the ith body.
These equations exist for each of the three bodies in the problem. G is Newton's gravitational constant, the distance between body i and body j is given by rij, and each body i has mass mi. We need to take the equations in Example 7 along with an initial set of values for the positions and velocities (x1, y1, z1, vx1, vy1, vz1, and so on) and integrate these equations numerically to obtain the new values for the variables after a time step t. One sample integration of this system is presented in this article. A second solution is shown in Figure 5.
--D.J.C.
Figure 1: Integrator class hierarchy.
Figure 2: Derivative class structure.
Figure 3: Results of integration of the symmetric three-body system.
Figure 4: The three-body system.
Figure 5: Results of integration of an asymmetric three-body system.
Example 1: Differential equations solved by the integrator class.
Example 2: Generating the ith stage when you have the preceding i-1 stages.
Example 3: First stage generated (from Example 2).
Example 4: Integration step.
Example 5: Calculating the error estimate.
Example 6: Performing adaptation.
Example 7: Differential equations for the ith body of the three-body problem.
Table 1: Coefficients for the RungeKutta45 class.
i = 1 2 3 4 5 6
ai 0 1/5 3/10 3/5 1 7/8
bi1
bi2 1/5
bi3 3/40 4/40
bi4 3/10 -9/10 6/5
bi5 -11/54 5/2 -70/27 35/27
bi6 1631/55296 175/512 575/13824 44275/110592 253/4096
ci 37/378 0 250/621 125/594 0 512/1771
ci* 2825/27648 0 18575/48384 13525/55296 277/14336 1/4
Table 2: Comparison of final position from Runge-Kutta 4(5) and Runge-Kutta 7(8) integrations.
Runge-Kutta 4(5) Runge-Kutta 7(8) Difference
x1 4697.00964137395 4697.00964144131 -0.00000006736
y1 36.5844745040766 36.5844743862259 0.0000001178507
z1 0 0 0
x2 -2316.82173637934 -2316.82173651791 0.00000013857
y2 -4086.02190850241 -4086.02190850171 -0.0000000007
z2 0 0 0
x3 -2380.18790499465 -2380.18790492343 -0.00000007122
y3 4049.43743399842 4049.43743411535 -0.00000011693
z3 0 0 0
// integrat.h -- Integrator class definition
#ifndef INTEGRAT_H
#define INTEGRAT_H
#include <math.h>
#include <stdio.h>
#include "derivs.h"
#define FALSE 0
#define TRUE 1
class Integrator {
protected:
int dimension; // Size of the "state" vector
double *initState; // Pointer to state for integration
double *finalState;// Pointer to integrated state
// Pointer to derivative class
Derivative *ddt;
public:
Integrator(int dim);
virtual ~Integrator(void);
// Access to the integration variables
void InitState(double *istate)
{ initState = istate; }
void FinalState(double *fstate)
{ finalState = fstate; }
void Derivatives(Derivative *deriv);
virtual double Step(double h) = 0;
};
#endif
// integrat.cpp -- Integrator class functions
#include "integrat.h"
Integrator::Integrator(int dim)
{
dimension = dim;
}
Integrator::~Integrator(void)
{
return;
}
void Integrator::Derivatives(Derivative *deriv)
{
ddt = deriv;
}
// derivs.h -- Class for derivative skeleton for the integrator class
#ifndef DERIVS_H
#define DERIVS_H
#define TRUE 1
#define FALSE 0
#include <math.h>
class Derivative {
private:
// No private data or functions in the skeleton
public:
Derivative(void)
{ }
virtual ~Derivative(void)
{ }
virtual int GetDerivatives(double dt, double *where) = 0;
};
#endif
// 3body.h -- Class for gravitational forces for the example problem
#ifndef THREEBODY_H
#define THREEBODY_H
#include "derivs.h"
class ThreeBody : public Derivative {
private:
double mu1, mu2, mu3;
double r1[3], r2[3], r3[3];
double v1[3], v2[3], v3[3];
public:
ThreeBody(void);
~ThreeBody(void);
int GetDerivatives(double dt, double *where);
void Mu1(double m1) { mu1 = m1; }
void Mu2(double m2) { mu2 = m2; }
void Mu3(double m3) { mu3 = m3; }
};
#endif
// 3body.cpp -- gravitational forces for the example problem; provides
// 1st derivative information of the state. Tailor code to match your problem.
#include "3body.h"
ThreeBody::ThreeBody(void)
{
mu1 = 100;
mu2 = 100;
mu3 = 100;
}
ThreeBody::~ThreeBody(void)
{
return;
}
int ThreeBody::GetDerivatives(double dt, double *where)
// Note: dt is not used here; it may be needed in other systems
{
double d12, d23, d13; // Distances from origin
double d12cubed, d23cubed, d13cubed;
// Distances from origin cubed
// where comes in with positions and velocities
r1[0] = where[0]; r1[1] = where[1]; r1[2] = where[2];
v1[0] = where[3]; v1[1] = where[4]; v1[2] = where[5];
r2[0] = where[6]; r2[1] = where[7]; r2[2] = where[8];
v2[0] = where[9]; v2[1] = where[10]; v2[2] = where[11];
r3[0] = where[12]; r3[1] = where[13]; r3[2] = where[14];
v3[0] = where[15]; v3[1] = where[16]; v3[2] = where[17];
d12 = sqrt((r1[0] - r2[0]) * (r1[0] - r2[0]) +
(r1[1] - r2[1]) * (r1[1] - r2[1]) +
(r1[2] - r2[2]) * (r1[2] - r2[2]));
d12cubed = d12 * d12 * d12;
if (d12cubed == 0.0) // masses must be separated
return FALSE;
d23 = sqrt((r2[0] - r3[0]) * (r2[0] - r3[0]) +
(r2[1] - r3[1]) * (r2[1] - r3[1]) +
(r2[2] - r3[2]) * (r2[2] - r3[2]));
d23cubed = d23 * d23 * d23;
if (d23cubed == 0.0) // masses must be separated
return FALSE;
d13 = sqrt((r1[0] - r3[0]) * (r1[0] - r3[0]) +
(r1[1] - r3[1]) * (r1[1] - r3[1]) +
(r1[2] - r3[2]) * (r1[2] - r3[2]));
d13cubed = d13 * d13 * d13;
if (d13cubed == 0.0) // masses must be separated
return FALSE;
// and returns with velocities...
where[0] = v1[0]; where[1] = v1[1]; where[2] = v1[2];
where[6] = v2[0]; where[7] = v2[1]; where[8] = v2[2];
where[12] = v3[0]; where[13] = v3[1]; where[14] = v3[2];
// and accelerations
// m1
where[3] = - mu2 * (r1[0] - r2[0]) / d12cubed
- mu3 * (r1[0] - r3[0]) / d13cubed;
where[4] = - mu2 * (r1[1] - r2[1]) / d12cubed
- mu3 * (r1[1] - r3[1]) / d13cubed;
where[5] = - mu2 * (r1[2] - r2[2]) / d12cubed
- mu3 * (r1[2] - r3[2]) / d13cubed;
// m2
where[9] = - mu1 * (r2[0] - r1[0]) / d12cubed
- mu3 * (r2[0] - r3[0]) / d23cubed;
where[10] = - mu1 * (r2[1] - r1[1]) / d12cubed
- mu3 * (r2[1] - r3[1]) / d23cubed;
where[11] = - mu1 * (r2[2] - r1[2]) / d12cubed
- mu3 * (r2[2] - r3[2]) / d23cubed;
// m3
where[15] = - mu1 * (r3[0] - r1[0]) / d13cubed
- mu2 * (r3[0] - r2[0]) / d23cubed;
where[16] = - mu1 * (r3[1] - r1[1]) / d13cubed
- mu2 * (r3[1] - r2[1]) / d23cubed;
where[17] = - mu1 * (r3[2] - r1[2]) / d13cubed
- mu2 * (r3[2] - r2[2]) / d23cubed;
return TRUE;
}
// adapt_rk.h -- Runge-Kutta class definition
#ifndef ADAPT_RK_H
#define ADAPT_RK_H
#include "integrat.h"
class AdaptiveRungeKutta : public Integrator {
protected:
int stages; // Number deriv. evals in algorithm
double *ai; // Portion of step at ith stage
double **bij; // Portions of previous stages to use
double *cj; // For RK step taken
double *errorest; // Estimate of the error
int fixedStep; // Boolean flag for fixed step mode
double accuracy; // Maximum errorest for good step
double safety; // Safety factor for adaptive steps
double goodpower; // 1/(RK order)
double badpower; // 1/(RK order - 1)
// The following are internal workspaces for the algorithm
double *accumulated;
double *workstate;
double **stepk;
double GetError(double **stepk);
double NewStep(double maxError, double h);
double independent;// indep. variable (e.g. time)
public:
AdaptiveRungeKutta(int st, int dim);
virtual ~AdaptiveRungeKutta(void);
virtual void SetCoefficients(void) = 0;
double Step(double h);
void Accuracy(double goodness)
{ accuracy = goodness; }
void FixedStep(int tf)
{ if (tf == 0) fixedStep = FALSE; else fixedStep = TRUE; }
};
#endif
// adapt_rk.cpp -- Runge-Kutta functions
#include "adapt_rk.h"
AdaptiveRungeKutta::AdaptiveRungeKutta(int st, int dim) : Integrator(dim)
{
int i;
stages = st;
fixedStep = TRUE;
accuracy = 1e-7;
independent = 0.0;
// Set the pointers to the arrays to NULL for now
ai = new double[stages];
cj = new double[stages];
errorest = new double[stages];
bij = new double*[stages];
for (i = 0; i < stages; i++)
bij[i] = new double[stages];
accumulated = new double[dimension]; // Working arrays
workstate = new double[dimension];
stepk = new double*[stages];
for (i = 0; i < stages; i++)
stepk[i] = new double[dimension];
}
AdaptiveRungeKutta::~AdaptiveRungeKutta(void)
{
int i;
if (ai)
delete [] ai;
if (bij) {
for (i = 0; i < stages; i++)
delete [] bij[i];
delete [] bij;
}
if (cj)
delete [] cj;
if (errorest)
delete [] errorest;
if (accumulated)
delete [] accumulated;
if (workstate)
delete [] workstate;
if (stepk) {
for (i = 0; i < stages; i++)
delete[] stepk[i];
delete [] stepk;
}
}
double AdaptiveRungeKutta::Step(double h)
{
double stepsize = h; // Next stepsize to take
double desiredStep = h;
double stepSoFar = 0.0; // Step taken to date
double stepTaken; // Current step
double maxErrorFound = 0.0;
int i, j, k;
int fixedPassed = TRUE;
if (fixedStep) // If in fixed step mode, don't
fixedPassed = FALSE; // break loop until desired step
for (i = 0; i < dimension; i++)
workstate[i] = initState[i];
do {
stepTaken = stepsize;
for (i = 0; i < stages; i++) {
// Fill accumulated with info from the previous stages
for (j = 0; j < dimension; j++) {
accumulated[j] = workstate[j];
for (k = 0; k < i; k++) {
accumulated[j] += bij[i][k] * stepk[k][j];
}
}
// Calculate the data for the current stage
if (!ddt -> GetDerivatives(independent + stepTaken * ai[i],
accumulated))
return 0.0; // Bad derivative ==> No step taken
for (j = 0; j < dimension; j++)
stepk[i][j] = stepTaken * accumulated[j];
}
// Here's the step used:
for (j = 0; j < dimension; j++) {
finalState[j] = workstate[j];
for (i = 0; i < stages; i++)
finalState[j] += cj[i] * stepk[i][j];
}
// Here we test the step
maxErrorFound = GetError(stepk);
// And calculate the next step to be taken
stepsize = NewStep(maxErrorFound, stepTaken);
// Now accumulate for fixed step mode
if (fixedStep) {
// Step was good
if (maxErrorFound < accuracy) {
stepSoFar += stepTaken;
if (stepsize > desiredStep - stepSoFar)
stepsize = desiredStep - stepSoFar;
// Update the working state if not at the end
if (stepSoFar < desiredStep)
for (i = 0; i < dimension; i++)
workstate[i] = finalState[i];
// If at end, set flag
else
fixedPassed = TRUE;
}
}
// Stop if step is good, else try again with the new step.
} while ((maxErrorFound > accuracy) || !fixedPassed);
if (fixedStep)
return stepSoFar;
return stepTaken;
}
double AdaptiveRungeKutta::GetError(double **stepk)
{
int i, j;
double biggestError = 0.0, currentError;
for (i = 0; i < dimension; i++) {
currentError = 0.0;
for (j = 0; j < stages; j++)
currentError += errorest[j] * stepk[j][i];
currentError = fabs(currentError);
if (currentError > biggestError)
biggestError = currentError;
}
return biggestError;
}
double AdaptiveRungeKutta::NewStep(double maxError, double h)
{
double newstep;
if (accuracy >= maxError) // Step was good
newstep = safety * h * fabs(pow(accuracy/maxError, goodpower));
else // Step was too big
newstep = safety * h * fabs(pow(accuracy/maxError, badpower));
return newstep;
}
Copyright © 1995, Dr. Dobb's Journal