Ordinary Differential Equations

A system of ordinary differential equations (ODE's) relates the derivatives of a set of function to the function values and a time variable.

  • y' = F(t, y)

This section deals with the integration of such systems. A stiff system is a system whose numerical solution may be unstable. Special precautions need to be taken to avoid getting completely meaningless results. Numerics.NET contains classes for dealing with stiff and non-stiff systems.

Using the ODE integrator classes

All classes that implement integration algorithms inherit from the OdeIntegrator class. This class defines methods and properties common to all implementations.

The following integrators are available:

Defining the System

A system of differential equations must be supplied as a DifferentialFunction delegate. This is a delegate that takes three arguments. The first is the time value. The second is a Vector<T> that contains the current solution, and the third is another Vector<T> that is to contain the derivative of the solution. The same value should be returned as well.

The example below defines the system for the famous Lorenz attractor, defined by:

  • y1' = σ(y2 - y1)
  • y2' = y1(ρ - y3) - y2
  • y3' = y1y2 - βy3

where σ, ρ and β are parameters.

C# 
static Vector<double> LorenzFunction(double t, Vector<double> y, Vector<double> dy)
{
    if (dy == null)
        dy = Vector.Create<double>(3);

    double sigma = 10.0;
    double beta = 8 / 3.0;
    double rho = 28.0;

    dy[0] = sigma * (y[1] - y[0]);
    dy[1] = y[0] * (rho - y[2]) - y[1];
    dy[2] = y[0] * y[1] - beta * y[2];

    return dy;
}

To set up an OdeIntegrator object to integrate a system, assign the delegate to the object's DifferentialFunction property.

To define the problem fully, the initial conditions need to be specified. The initial time is set through the InitialTime property. This is a Double value. The initial value is set through the InitialValue property. This is a Vector<T> that contains the solution at the initial time. Continuing the above example:

C# 
ClassicRungeKuttaIntegrator rk4 = new ClassicRungeKuttaIntegrator();
rk4.DifferentialFunction = LorenzFunction;
rk4.InitialTime = 0.0;
rk4.InitialValue = Vector.Create(1.0, 2.0, 3.0);

Higher order systems

Many differential equations contain expressions in terms of the second and sometimes higher derivatives. These can be converted to an equivalent first order system by introducing new variables for the first (and possibly higher) derivatives. For example, the equation for a spring:

x'' = -k/m x

can be rewritten as a system of two equations:

y' = -k/m x

x' = y

Setting integration parameters

With the system of equations defined, the next step is to define the parameters of the integration.

The AbsoluteTolerance and RelativeTolerance determine the largest integration error allowed in each step of the calculation. Methods with variable step size attempt to keep the estimated error within both tolerances. For methods with fixed step size, these properties have no effect.

The default value for AbsoluteTolerance is 10-4, and for RelativeTolerance it is 10-6.

The InitialStepsize property specifies the step size to use for the first integration step. For methods that use a fixed step size, like the classic Runge-Kutta integrator, this is also the step size that will be used throughout the integration.

Some care should be taken in choosing this value. If it is too large, the errors in the solution will build up quickly. If it is too small, the integration may take an excessive amount of time.

Variable step methods, which include both the RKF and the CVODE integrator, will adjust the step size to keep the erro within the specified tolerance. Two additional properties can affect the calculation. The MinStepSize and MaxStepSize properties specify limits on the step size during the calculation.

Running the integration

Once the system has been defined and the parameters have been set, the system can be integrated.

Differential equations are usually integrated in one of two ways. In some cases, the interest is only in the solution at a final time. In other cases, the interest is in the evolution of the solution from the initial to the final time. The integrator classes contain methods that accommodate both scenarios.

The Integrate method performs the integration up to the specified final time, which is supplied as the first and only argument. As many steps as are necessary are taken. If the time of the final step is past the final time, as is often the case, the value at the final time is interpolated using the known results.

Two other integration methods are available. Each also takes a single argument that specifies the final time. As its name implies, the IntegrateSingleStep method performs a single integration step. The step is taken in the direction of the final time. The return value is the value at the step, regardless of whether it is before or after the final time. The IntegrateMultipleSteps computes full integration steps until at least the final time. This is similar to the Integrate method without arguments, except that the return value is the value at the step. No interpolation is performed.

The first call to any of these methods uses the values of the InitialTime and InitialValue properties as starting values. Subsequent steps start from the last final time. The combination of this mechanism and the fact that results are interpolated if the grid point is past the final time makes it easy and efficient to generate intermediate results at fixed intervals.

The following example uses the CVODE integrator to compute and print 100 points of the Lorentz attractor path.

C# 
CvodeIntegrator integrator = CvodeIntegrator();
integrator.InitialTime = 0;
integrator.InitialValue = Vector.Create(1.0, 2.0, 3.0);
integrator.DifferentialFunction = LorenzFunction;
integrator.InitialStepsize = 1e-3;
integrator.AbsoluteTolerance = 1e-8;
for (int i = 1; i <= 100; i++) {
    double t = i / 10.0;
    integrator.Integrate(t);
    Vector y = integrator.CurrentValue;
    Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3,10:F6}", 
        integrator.CurrentTime, y[0], y[1], y[2]);
}

To restart the calculation from new starting values, use the Restart. Without arguments, the integration is restarted using the InitialTime and InitialValue properties. It is also possible to pass in a new start time and vector of initial values.

Integration Algorithms

The Classic 4th Order Runge-Kutta Method

The best known method for integrating systems of ODE's is probably the 4th order Runge-Kutta method. This is due, no doubt, to its simplicity. But that simplicity comes at a cost. The stepsize is fixed, so there is no way to manage the error.

The ClassicRungeKuttaIntegrator class implements this algorithm.

The 4th-5th Order Runge-Kutta-Fehlberg Method

The Runge-Kutta-Fehlberg method uses a combination of a 4th and 5th order method to estimate the integration error in each step. A new step size is computed based on this information. The ability to dynamically change the step size makes this method suitable for most non-stiff problems.

The RungeKuttaFehlbergIntegrator class implements this algorithm.

The CVODE Integrator

The CVODE integrator is part of the SUNDIALS suite of nonlinear and differential/algebraic equation solvers developed at the Lawrence Livermore National Laboratory. CVODE can be used to solve both stiff and non-stiff systems of ordinary differential equations. It is the integrator of choice when high precision is required, or when the system of ODE's is stiff.

The CvodeIntegrator class is based on the original CVODE code.

Stiff Problems

A system of differential equations is stiff when the numerical integration may be unstable. As a result, the step size needed to produce accurate results is very small and the computation is very slow.

Techniques exist to handle stiff systems. The CVODE integrator implements such a technique. We will use the Van Der Pol equation to illustrate. This is a second order equation:

  • y'' = μ(1 - y2)y' - y

This is a second order equation. The first step is to rewrite it as a system of first order equations. We do this by setting y1 = y and y2 = y'. The above equation then becomes:

  • y1' = y2
  • y2' = μ(1 - y12)y2 - y1
C# 
static Vector<double> VanDerPolFunction(double t, Vector<double> y, Vector<double> dy)
{
    if (dy == null)
        dy = Vector.Create<double>(2);

    double mu = 1000.0;

    dy[0] = y[1];
    dy[1] = mu * (1 - y[0] * y[0]) * y[1] - y[0];

    return dy;
}

We can integrate this system using a Runge-Kutta-Fehlberg variable step method:

C# 
RungeKuttaFehlbergIntegrator nonstiff = new RungeKuttaFehlbergIntegrator();
nonstiff.InitialTime = 0;
nonstiff.InitialValue = Vector.Create(2.0, 0.0);
nonstiff.DifferentialFunction = VanDerPolFunction;
nonstiff.InitialStepsize = 1e-3;
for (int i = 1; i <= 10; i++)
{
    double t = i * 25.0;
    nonstiff.Integrate(t);
    Vector<double> y = nonstiff.CurrentValue;
    Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3}",
        nonstiff.CurrentTime, y[0], y[1], nonstiff.IterationsNeeded);
}

When we run this code, we find that thousands of steps are necessary in each iteration. Fortunately, the CVODE integrator has the option of using a method specifically designed to handle stiff problems. One slight complication with this method is that it requires the Jacobian of the differential function. The Jacobian is a square matrix whose elements are the partial derivatives of the differential function with respect to each component of the solution. The Jacobian must be supplied as a DifferentialJacobianFunction delegate. For the Van Der Pol equation, the Jacobian is:

C# 
static Matrix<double> VanDerPolJacobian(double t, Vector<double> y, Vector<double> dy, Matrix<double> J)
{
    if (J == null)
        J = Matrix.Create<double>(2, 2);

    double mu = 1000.0;

    J[0, 0] = 0.0;
    J[1, 0] = -2 * mu * y[0] * y[1] - 1.0;

    J[0, 1] = 1.0;
    J[1, 1] = mu * (1 - y[0] * y[0]);

    return J;
}

Armed with the Jacobian, we can create a more efficient integrator for the Van Der Pol system. Two modifications are necessary. First, we must call the constructor of the CvodeIntegrator integrator with a parameter that specifies that the system we want to integrate is stiff. The parameter is a OdeKind. Second, we must pass the Jacobian to the integrator by calling the SetDenseJacobian method. This gives us:

C# 
CvodeIntegrator stiff = new CvodeIntegrator(OdeKind.Stiff);
stiff.InitialTime = 0;
stiff.InitialValue = Vector.Create(2.0, 0.0);
stiff.DifferentialFunction = VanDerPolFunction;
stiff.SetDenseJacobian(VanDerPolJacobian);
for (int i = 1; i <= 10; i++)
{
    double t = i * 25.0;
    stiff.Integrate(t);
    Vector<double> y = nonstiff.CurrentValue;
    Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3}",
            stiff.CurrentTime, y[0], y[1], stiff.IterationsNeeded);
}

This integrator takes only a handful of steps in each iteration, and uses only a few hundred steps to compute the solution at the most difficult point. The total running time was cut down by a factor of around 100 or more.