**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())