A system of ordinary differential equations (ODE's) relates the
derivatives of a set of function to the function values and
a time variable.
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.
The Extreme Optimization Numerical Libraries for .NET
contains classes for dealing with stiff and nonstiff 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:
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 VectorT
that contains the current solution, and the third is another
VectorT
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:

y_{1}' = σ(y_{2}  y_{1})

y_{2}' = y_{1}(ρ  y_{3})  y_{2}

y_{3}' = y_{1}y_{2}  βy_{3}
where σ, ρ and β are parameters.
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;
}
Shared Function LorenzFunction(t As Double, y As Vector(Of Double), dy As Vector(Of Double)) As Vector(Of Double)
If (dy Is Nothing) Then
dy = Vector.Create(Of Double)(3)
End If
Dim sigma = 10.0
Dim beta = 8 / 3.0
Dim 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
End Function
No code example is currently available or this language may not be supported.
let LorenzFunction t (y : Vector<double>) (dy : Vector<double>) =
let sigma = 10.0
let beta = 8.0 / 3.0
let 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]
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 VectorT that contains
the solution at the initial time. Continuing the above example:
ClassicRungeKuttaIntegrator rk4 = new ClassicRungeKuttaIntegrator();
rk4.DifferentialFunction = LorenzFunction;
rk4.InitialTime = 0.0;
rk4.InitialValue = Vector.Create(1.0, 2.0, 3.0);
Dim rk4 As ClassicRungeKuttaIntegrator = New ClassicRungeKuttaIntegrator()
rk4.DifferentialFunction = AddressOf LorenzFunction
rk4.InitialTime = 0.0
rk4.InitialValue = Vector.Create(1.0, 2.0, 3.0)
No code example is currently available or this language may not be supported.
let rk4 = new ClassicRungeKuttaIntegrator()
rk4.DifferentialFunction < 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 RungeKutta 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.
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.
CvodeIntegrator integrator = CvodeIntegrator();
integrator.InitialTime = 0;
integrator.InitialValue = Vector.Create(1.0, 2.0, 3.0);
integrator.DifferentialFunction = LorenzFunction;
integrator.InitialStepsize = 1e3;
integrator.AbsoluteTolerance = 1e8;
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]);
}
Dim integrator As New CvodeIntegrator()
integrator.InitialTime = 0
integrator.InitialValue = Vector.Create(1.0, 2.0, 3.0)
integrator.DifferentialFunction = LorenzFunction
integrator.InitialStepsize = 1e3
integrator.AbsoluteTolerance = 1e8
For i As Integer = 1 To 100
Dim t As Double = i / 10.0
integrator.Integrate(t)
Dim y As Vector = integrator.CurrentValue
Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3,10:F6}", _
integrator.CurrentTime, y(0), y(1), y(2))
Next
No code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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.
The Classic 4th Order RungeKutta Method
The best known method for integrating systems of ODE's is
probably the 4th order RungeKutta 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 4th5th Order RungeKuttaFehlberg Method
The RungeKuttaFehlberg 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 nonstiff problems.
The RungeKuttaFehlbergIntegrator
class implements this algorithm.
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 nonstiff 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.
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:
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
y_{1} = y and y_{2} = y'.
The above equation then becomes:

y_{1}' = y_{2}

y_{2}' = μ(1  y_{1}^{2})y_{2}  y_{1}
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;
}
Shared Function VanDerPolFunction(t As Double, y As Vector(Of Double), dy As Vector(Of Double)) As Vector(Of Double)
If (dy Is Nothing) Then dy = Vector.Create(Of Double)(2)
Dim mu As Double = 1000.0
dy(0) = y(1)
dy(1) = mu * (1  y(0) * y(0)) * y(1)  y(0)
Return dy
End Function
No code example is currently available or this language may not be supported.
let VanDerPolFunction t (y : Vector<double>) (dy : Vector<double>) =
let mu = 1000.0
dy.[0] < y.[1]
dy.[1] < mu * (1.0  y.[0] * y.[0]) * y.[1]  y.[0]
dy
We can integrate this system using a RungeKuttaFehlberg variable step method:
RungeKuttaFehlbergIntegrator nonstiff = new RungeKuttaFehlbergIntegrator();
nonstiff.InitialTime = 0;
nonstiff.InitialValue = Vector.Create(2.0, 0.0);
nonstiff.DifferentialFunction = VanDerPolFunction;
nonstiff.InitialStepsize = 1e3;
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);
}
Dim nonstiff As RungeKuttaFehlbergIntegrator = New RungeKuttaFehlbergIntegrator()
nonstiff.InitialTime = 0
nonstiff.InitialValue = Vector.Create(2.0, 0.0)
nonstiff.DifferentialFunction = AddressOf VanDerPolFunction
nonstiff.InitialStepsize = 0.001
For i = 1 To 10
Dim t As Double = i * 25.0
nonstiff.Integrate(t)
Dim y As Vector(Of Double) = nonstiff.CurrentValue
Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3}",
nonstiff.CurrentTime, y(0), y(1), nonstiff.IterationsNeeded)
Next
No code example is currently available or this language may not be supported.
let nonstiff = new RungeKuttaFehlbergIntegrator()
nonstiff.InitialTime < 0.0
nonstiff.InitialValue < Vector.Create(2.0, 0.0)
nonstiff.DifferentialFunction < DifferentialFunction(VanDerPolFunction)
nonstiff.InitialStepsize < 1e3
for i in 1..10 do
let t = (float i) * 25.0
nonstiff.Integrate(t) > ignore
let 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:
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;
}
Shared Function VanDerPolJacobian(t As Double, y As Vector(Of Double), dy As Vector(Of Double), J As Matrix(Of Double)) As Matrix(Of Double)
If (J Is Nothing) Then J = Matrix.Create(Of Double)(2, 2)
Dim mu As Double = 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
End Function
No code example is currently available or this language may not be supported.
let VanDerPolJacobian t (y : Vector<double>) (dy : Vector<double>) (J : Matrix<double>) =
let mu = 1000.0
J.[0, 0] < 0.0
J.[1, 0] < 2.0 * mu * y.[0] * y.[1]  1.0
J.[0, 1] < 1.0
J.[1, 1] < mu * (1.0  y.[0] * y.[0])
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:
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);
}
Dim stiff As CvodeIntegrator = New CvodeIntegrator(OdeKind.Stiff)
stiff.InitialTime = 0
stiff.InitialValue = Vector.Create(2.0, 0.0)
stiff.DifferentialFunction = AddressOf VanDerPolFunction
stiff.SetDenseJacobian(AddressOf VanDerPolJacobian)
For i = 1 To 10
Dim t As Double = i * 25.0
stiff.Integrate(t)
Dim y As Vector(Of Double) = nonstiff.CurrentValue
Console.WriteLine("{0:F3} {1,10:F6} {2,10:F6} {3}",
stiff.CurrentTime, y(0), y(1), stiff.IterationsNeeded)
Next
No code example is currently available or this language may not be supported.
let stiff = new CvodeIntegrator(OdeKind.Stiff)
stiff.InitialTime < 0.0
stiff.InitialValue < Vector.Create(2.0, 0.0)
stiff.DifferentialFunction < DifferentialFunction(VanDerPolFunction)
stiff.SetDenseJacobian(DifferentialJacobianFunction(VanDerPolJacobian))
for i in 1..10 do
let t = (float i) * 25.0
stiff.Integrate(t) > ignore
let 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.