C# Runge Kutta Solver Example

← All NMath Code Examples

 

using System;
using System.Collections.Generic;
using System.Diagnostics;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Analysis;


namespace CenterSpace.NMath.Analysis.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 'Key' element of the pair is
      // the time span vector, the second 'Value' element 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 'Key' element of the pair is
      // the time span vector, the second 'Value' element 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
Top