← All NMath Code Examples
using System;
using System.Collections.Generic;
using System.Diagnostics;
using CenterSpace.NMath.Core;
namespace CenterSpace.NMath.Examples.CSharp
{
/// <summary>
/// Example showing how to use the RungeKutta45OdeSolver to solve
/// a nonstiff set 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
/// </summary>
class RungeKutta45OdeSolverExample
{
static void Main( string[] args )
{
Console.WriteLine();
// Simple example solving the system of differential equations of
// motion for a rigid body with no external forces. The equations
// are defined above by the function Rigid().
Console.WriteLine( "Simple Solve -----------------------------------" );
Console.WriteLine();
SimpleSolve();
Console.WriteLine( \n);
// A mass matrix and an output function are options that may
// be specified for the solver. The output function is a callback
// that is invoked at initialization, each integration step and
// after the conclusion of the last integration step. Thus the
// solution process can be monitored, and even terminated using
// the output function.
Console.WriteLine( "Output Function and Mass Matrix Solve ---------" );
Console.WriteLine();
WithOutputFunctionAndMassMatrix();
Console.WriteLine( \n);
// The above examples use the NMath class RungeKutta45OdeSolver. This
// class uses an adaptive algorithm to figure out the optimal time step
// size at each integration step. NMath also supplies a class,
// RungeKutta5OdeSolver that does NOT use an adaptive step size
// algorithm, allowing the user to specify the time values for
// each integration step.
Console.WriteLine( "Non-adaptive Step Size -------------------------" );
Console.WriteLine();
NonAdaptiveStepSizeExample();
Console.WriteLine();
Console.WriteLine( "Press Enter Key" );
Console.Read();
}
static void SimpleSolve()
{
// Construct the solver.
var solver = new RungeKutta45OdeSolver();
// 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.
var timeSpan = new DoubleVector( 0.0, 12.0 );
// Initial y vector.
var y0 = new DoubleVector( 0.0, 1.0, 1.0 );
// Construct solver options. Here we set the absolute and relative tolerances to use.
// At the ith integration step the error, e[i] for the estimated solution
// y[i] satisfies
// e[i] <= max(RelativeTolerance * Math.Abs(y[i]), AbsoluteTolerance[i])
// The solver can increase the number of output points by a specified factor
// called Refine (useful for creating smooth plots). The default value is
// 4. Here we set the Refine value to 1, meaning we do not wish any
// additional output points to be added by the solver.
RungeKutta45OdeSolver.Options solverOptions = new RungeKutta45OdeSolver.Options
{
AbsoluteTolerance = new DoubleVector( 1e-4, 1e-4, 1e-5 ),
RelativeTolerance = 1e-4,
Refine = 1
};
// Construct the delegate representing our system of differential equations...
var odeFunction = new Func<double, DoubleVector, DoubleVector>( Rigid );
// ...and solve. The solution is returned as a key/value pair. The first Keyelement of the pair is
// the time span vector, the second Valueelement of the pair is the corresponding solution values.
// That is, if the computed solution function is y then
// y(soln.Key[i]) = soln.Value[i]
RungeKutta45OdeSolver.Solution<DoubleMatrix> soln = solver.Solve( odeFunction, timeSpan, y0, solverOptions );
Console.WriteLine( "T = " + soln.T.ToString( "G5" ) );
Console.WriteLine( "Y = " );
Console.WriteLine( soln.Y.ToTabDelimited( "G5" ) );
}
/// <summary>
/// Function describing the system of differential equations.
/// </summary>
/// <param name="t">Time parameter.</param>
/// <param name="y">State vector.</param>
/// <returns>The vector of values of the derivatives of the functions
/// at a specified time t and function values y:
/// y1= y2*y3,
/// y2= -y1*y3,
/// y3= -0.51*y1*y2
/// </returns>
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;
}
/// <summary>
/// The output function to be called by the ODE solver at each integration step.
/// If this function returns true the solver continues with the next step, if
/// it returns false the solver is halted.
/// This output function just outputs to the console upon each invocation.
/// </summary>
/// <param name="t">The time value at the current integration step. If flag is equal to Initialize,
/// this will be the initial time value t0.</param>
/// <param name="y">The calculated function value at the current integration step. If flag is equal to Initialize,
/// this will be the initial function value y0.</param>
/// <param name="flag">Flag indicating what stage the solver is at:
/// Initialize - before solving begins. y and t values are the problems initial values.
/// IntegrationStep - just finished an integration step. y and t values are the values
/// calculated at that step.
/// Done - just finished last step. y and t values are the last values in the returned
/// solution.</param>
/// <returns>true if the solver is to proceed with the calculation, false forces the solver
/// to stop.</returns>
static bool OutputFunction( DoubleVector t, DoubleMatrix y, RungeKutta45OdeSolver.OutputFunctionFlag flag )
{
switch ( flag )
{
case RungeKutta45OdeSolver.OutputFunctionFlag.Initialize:
Console.WriteLine( "Output function : Initialize" );
return true;
case RungeKutta45OdeSolver.OutputFunctionFlag.IntegrationStep:
Console.WriteLine( "Output function Integration step:" );
Console.WriteLine( "t = " + t.ToString( "G5" ) );
Console.WriteLine( "y = " );
Console.WriteLine( y.ToTabDelimited( "G5" ) );
return true;
case RungeKutta45OdeSolver.OutputFunctionFlag.Done:
Console.WriteLine( "Output function: Done" );
return true;
default:
Console.WriteLine( "Output function: Unknown flag" );
return false;
}
}
static void WithOutputFunctionAndMassMatrix()
{
// Construct the solver.
var solver = new RungeKutta45OdeSolver();
// 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.
var timeSpan = new DoubleVector( 0.0, 12.0 );
// Initial y vector.
var y0 = new DoubleVector( 0.0, 1.0, 1.0 );
// Construct solver options. Here we set the absolute and relative tolerances to use.
// At the ith integration step the error, e[i] for the estimated solution
// y[i] satisfies
// e[i] <= max(RelativeTolerance*Math.Abs(y[i]), AbsoluteTolerance[i])
RungeKutta45OdeSolver.Options solverOptions = new RungeKutta45OdeSolver.Options
{
AbsoluteTolerance = new DoubleVector( 1e-4, 1e-4, 1e-5 ),
OutputFunction = new Func<DoubleVector, DoubleMatrix, RungeKutta45OdeSolver.OutputFunctionFlag, bool>( OutputFunction ),
RelativeTolerance = 1e-4
};
// Construct the delegate representing our system of differential equations...
var odeFunction = new Func<double, DoubleVector, DoubleVector>( Rigid );
// ...and solve. The solution is returned as a key/value pair. The first Keyelement of the pair is
// the time span vector, the second Valueelement of the pair is the corresponding solution values.
// That is, if the computed solution function is y then
// y(soln.Key[i]) = soln.Value[i]
RungeKutta45OdeSolver.Solution<DoubleMatrix> soln = solver.Solve( odeFunction, timeSpan, y0, solverOptions );
// Print out a few values
Console.WriteLine( "Number of time values = " + soln.T.Length );
for ( int i = 0; i < soln.T.Length; i++ )
{
if ( i % 7 == 0 )
{
Console.WriteLine( "y({0}) = {1}", soln.T[i].ToString( "G5" ), soln.Y.Row( i ).ToString( "G5" ) );
}
}
// We can also specify a mass matrix in the solver options...
var M = new DoubleMatrix( "3x3[1 2 1 2 1 3 9 4 1]" );
solverOptions.MassMatrix = M;
soln = solver.Solve( odeFunction, timeSpan, y0, solverOptions );
// Print out a few values for the mass matrix solution.
Console.WriteLine( "\nWith constant mass matrix:" );
Console.WriteLine( "Number of time values = " + soln.T.Length );
for ( int i = 0; i < soln.T.Length; i++ )
{
if ( i % 7 == 0 )
{
Console.WriteLine( "y({0}) = {1}", soln.T[i].ToString( "G5" ), soln.Y.Row( i ).ToString( "G5" ) );
}
}
// You can also solve with default values for all the options
soln = solver.Solve( odeFunction, timeSpan, y0 );
}
static void NonAdaptiveStepSizeExample()
{
// Construct the solver.
var solver = new RungeKutta5OdeSolver();
// Construct the time span vector. Using the non-adaptive solve
// function, integration steps will be performed only at the points
// in the time span vector.
var timeSpan = new DoubleVector( 85, 0, 0.15 );
// Initial y vector.
var y0 = new DoubleVector( 0.0, 1.0, 1.0 );
// Construct the delegate representing our system of differential equations...
var odeFunction = new Func<double, DoubleVector, DoubleVector>( Rigid );
// ...and solve. Each row of the returned matrix contains the calculated function values at the
// corresponding point in the input time span vector.
var soln = solver.Solve( odeFunction, timeSpan, y0 );
// Print out a few values
Console.WriteLine( "Number of time values = " + timeSpan.Length );
for ( int i = 0; i < timeSpan.Length; i++ )
{
if ( i % 7 == 0 )
{
Console.WriteLine( "y({0}) = {1}", timeSpan[i].ToString( "G5" ), soln.Row( i ).ToString( "G5" ) );
}
}
// Solve with a constant mass matrix
var M = new DoubleMatrix( "3x3[1 2 1 2 1 3 9 4 1]" );
var solverOptions = new RungeKutta5OdeSolver.Options();
solverOptions.MassMatrix = M;
soln = solver.Solve( odeFunction, timeSpan, y0, solverOptions );
// Print out a few values for the mass matrix solution.
Console.WriteLine( "\nWith constant mass matrix:" );
Console.WriteLine( "Number of time values = " + timeSpan.Length );
for ( int i = 0; i < timeSpan.Length; i++ )
{
if ( i % 7 == 0 )
{
Console.WriteLine( "y({0}) = {1}", timeSpan[i].ToString( "G5" ), soln.Row( i ).ToString( "G5" ) );
}
}
// Solve with an output function that is a class member function.
var outputFunctionObject = new SimpleNonAdaptiveOutput();
solverOptions.OutputFunction = new Func<double, DoubleVector, RungeKutta45OdeSolver.OutputFunctionFlag, bool>( outputFunctionObject.OutputFunction );
soln = solver.Solve( odeFunction, timeSpan, y0, solverOptions );
}
}
/// <summary>
/// Class containing an example output function for the Dormand Prince ODE non-adaptive solve
/// function. The class does nothing but accumulate the solutions at each solver
/// integration stop.
/// </summary>
public class SimpleNonAdaptiveOutput
{
/// <summary>
/// Gets and sets a list of the time values.
/// </summary>
public List<double> T { get; set; }
/// <summary>
/// Gets and sets a list of the calculated function values.
/// </summary>
public List<DoubleVector> Y { get; set; }
/// <summary>
/// Gets and sets a flag indicating whether or not the output
/// function/class has been initialized.
/// </summary>
public bool IsInitialized { get; set; }
/// <summary>
/// Flag indicating the output function is done being invoked.
/// </summary>
public bool IsDone { get; set; }
/// <summary>
/// Constructs and initializes a SimpleNonAdaptiveOutput object.
/// </summary>
public SimpleNonAdaptiveOutput()
{
IsInitialized = false;
IsDone = false;
T = null;
Y = null;
}
/// <summary>
/// The output function to be called by the ODE solver at each integration step. If this function returns true
/// the solver continues with the next step, if it returns false the solver is halted.
/// </summary>
/// <param name="t">The time value at the current integration step. If flag is equal to Initialize,
/// this will be the initial time value t0.</param>
/// <param name="y">The calculated function value at the current integration step. If flag is equal to Initialize,
/// this will be the initial function value y0.</param>
/// <param name="flag">Flag indicating what stage the solver is at:
/// Initialize - before solving begins. y and t values are the problems initial values.
/// IntegrationStep - just finished an integration step. y and t values are the values
/// calculated at that step.
/// Done - just finished last step. y and t values are the last values in the returned
/// solution.</param>
/// <returns>true if the solver is to proceed with the calculation, false forces the solver
/// to stop.</returns>
public bool OutputFunction( double t, DoubleVector y, RungeKutta45OdeSolver.OutputFunctionFlag flag )
{
if ( flag == RungeKutta45OdeSolver.OutputFunctionFlag.Initialize )
{
T = new List<double>();
T.Add( t );
Y = new List<DoubleVector>();
Y.Add( y );
IsInitialized = true;
IsDone = false;
return true;
}
if ( flag == RungeKutta45OdeSolver.OutputFunctionFlag.Done )
{
IsDone = true;
T.Add( t );
Y.Add( y );
return true;
}
if ( flag == RungeKutta45OdeSolver.OutputFunctionFlag.IntegrationStep )
{
T.Add( t );
Y.Add( y );
return true;
}
throw new InvalidArgumentException( "Unknown output function flag: " + flag );
}
}
}
← All NMath Code Examples
