NMath User's Guide

TOC | Previous | Next | Index

35.3 Dormand–Prince Method (.NET, C#, CSharp, VB, Visual Basic, F#)

The Dormand–Prince (DOPRI) method is a member of the Runge–Kutta family of ODE solvers. It uses six function evaluations to calculate fourth- and fifth-order accurate solutions. Dormand–Prince is currently the default method in MATLAB's ode45 solver. NMath provides two solvers which use Dormand-Prince methods:

● Class RungeKutta45OdeSolver solves an initial value ODE using an explicit Runge-Kutta (4,5) formula known as the Dormand-Prince pair.

● Class RungeKutta5OdeSolver solves an initial value ODE using a non-adaptive explicit Runge-Kutta formula of order 5. This is a non-adaptive solver. The step sequence is determined by given vector of time values, but the derivative function is evaluated multiple times per step. The solver implements the Dormand-Prince method of order 5 in a general framework of explicit Runge-Kutta methods.

For example, the following code shows how to use the RungeKutta45OdeSolver to solve a non-stiff system of equations describing the motion of a rigid body without external forces:

y1' = y2*y3,       y1(0) = 0
y2' = -y1*y3,      y2(0) = 1
y3' = -0.51*y1*y2, y3(0) = 1

This function describes the system of differential equations:

Code Example – C# ordinary differential equations (ODE)

static DoubleVector Rigid( double t, DoubleVector y )
{
  var dy = new DoubleVector( 3 );
  dy[0] = y[1] * y[2];
  dy[1] = -y[0] * y[2];
  dy[2] = -0.51 * y[0] * y[1];
  return dy;
}

Code Example – VB ordinary differential equations (ODE)

Shared Function Rigid(T As Double, Y As DoubleVector)
  Dim DY As New DoubleVector(3)
  DY(0) = Y(1) * Y(2)
  DY(1) = -Y(0) * Y(2)
  DY(2) = -0.51 * Y(0) * Y(1)
  Return DY
End Function

Parameter t is the time parameter, and y is the state vector.

First, construct the solver:

Code Example – C# ordinary differential equations (ODE)

var solver = new RungeKutta45OdeSolver();

Code Example – VB ordinary differential equations (ODE)

Dim Solver As New RungeKutta45OdeSolver()

Next, construct the time span vector. If this vector contains exactly two points, the solver interprets these to be the initial and final time values. Step size and function output points are provided automatically by the solver. Here the initial time value t0 is 0.0 and the final time value is 12.0.

Code Example – C# ordinary differential equations (ODE)

var timeSpan = new DoubleVector( 0.0, 12.0 );

Code Example – VB ordinary differential equations (ODE)

Dim TimeSpan As New DoubleVector(0.0, 12.0)

Specify the initial y vector.

Code Example – C# ordinary differential equations (ODE)

var y0 = new DoubleVector( 0.0, 1.0, 1.0 );

Code Example – VB ordinary differential equations (ODE)

Dim Y0 As New DoubleVector(0.0, 1.0, 1.0)

Optionally, construct solver options. Here we set the absolute and relative tolerances to use.

Code Example – C# ordinary differential equations (ODE)

var solverOptions =
  new RungeKutta45OdeSolver.Options
  {
    AbsoluteTolerance = new DoubleVector( 1e-4, 1e-4, 1e-5 ),
    RelativeTolerance = 1e-4,
    Refine = 1
  };

Code Example – VB ordinary differential equations (ODE)

Dim SolverOptions As New RungeKutta45OdeSolver.Options()
SolverOptions.AbsoluteTolerance =
  New DoubleVector(0.0001, 0.0001, 0.00001)
SolverOptions.RelativeTolerance = 0.0001
SolverOptions.Refine = 1

Construct the delegate representing our system of differential equations.

Code Example – C# ordinary differential equations (ODE)

var odeFunction =
  new Func<double, DoubleVector, DoubleVector>( Rigid );

Code Example – VB ordinary differential equations (ODE)

Dim ODEFunction As New Func(Of Double, DoubleVector, 
  DoubleVector)(AddressOf Rigid)

Finally, solve the problem. The solution is returned as a key/value pair. The first element of the pair is the time span vector, the second element is the corresponding solution values. That is, if the computed solution function is y then y(soln.Key[i]) = soln.Value[i].

Code Example – C# ordinary differential equations (ODE)

RungeKutta45OdeSolver.Solution<DoubleMatrix> soln =
  solver.Solve( odeFunction, timeSpan, y0, solverOptions );
Console.WriteLine( "T = " + soln.T );
Console.WriteLine( "Y = " );
Console.WriteLine( soln.Y.ToTabDelimited() );

Code Example – VB ordinary differential equations (ODE)

Dim Soln As RungeKutta45OdeSolver.Solution(Of DoubleMatrix) = 
Solver.Solve(ODEFunction, TimeSpan, Y0, SolverOptions)
Console.WriteLine("T = {0}", Soln.T)
Console.WriteLine("Y = ")
Console.WriteLine(Soln.Y.ToTabDelimited())

Top

Top