C# Runge Kutta Solver Example

using System;
using System.Collections.Generic;

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

namespace RungeKuttaExample
{
  class Program
  {
    /// <summary>
    /// Example illustrating the use of the RungeKuttaSolver class for solving
    /// first order initial value differential equations.
    /// </summary>
    /// <param name="args">Program arguments.</param>
    static void Main( string[] args )
    {
      // Consider the first order initial value problem
      // y' = (3 - 4*y)/2*x, y(1) = -4
      // which has solution
      // y = 0.75 - 19.0/(4*x^2)
      // To solve this differential equation using the 
      // Runge Kutta method, we first construct a FirstOrderInitialValueProblem
      // object which contains delegate for computing y' from x and y, and
      // the initial values.
      NMathFunctions.DoubleBinaryFunction f = delegate( double x, double y )
      {
        return ( 3.0 - 4.0 * y ) / ( 2.0 * x );
      };
      double x0 = 1.0;
      double y0 = -4.0;
      FirstOrderInitialValueProblem prob = new FirstOrderInitialValueProblem { F = f, Y0 = y0, X0 = x0 };

      // The Runge Kutta method returns the solution function, f, as a set of tabulated values
      // xi and f(xi). The initial x-value, x0 is the x0 of the initial condition. Number of
      // specified points and delta determine the remaining xi as
      // xi = x0 + i*delta
      int numPoints = 3000;
      double delta = .0005;

      // Construct a first order Runge Kutta solver (also known as Eulers method) and solve
      // our first order initial values differential equation.
      RungeKuttaSolver solver = new RungeKuttaSolver( numPoints, delta, RungeKuttaSolver.SolverOrder.First );
      KeyValuePair<double[], double[]> tabulatedValues = solver.Solve( prob );

      // Construct a delegate for the actual solution function for comparison.
      NMathFunctions.DoubleUnaryFunction actualSolutionFunction = delegate( double x )
      {
        return 0.75 - ( 19.0 / ( 4.0 * x * x ) );
      };

      // The solution to the differential equation is returned as a key-value pair. The
      // key is an array of doubles containing the x-values and the value is an array
      // containing the y = f(x) values. Below we compare the computed function values
      // to the know function values by computing the maximum absolute value of their
      // differences (l-infinity norm of the difference).
      DoubleVector xValues = new DoubleVector( tabulatedValues.Key );
      DoubleVector calculatedSolutionYvalues = new DoubleVector( tabulatedValues.Value );
      DoubleVector actualSolutionYvalues = xValues.Apply( actualSolutionFunction );
      double maxDiff = ( calculatedSolutionYvalues - actualSolutionYvalues ).InfinityNorm();
      Console.WriteLine( "Maximum error using first order Runge Kutta = " + maxDiff );

      // Solve the same problem with a second order Runge Kutta method and compute
      // the error. Note the error decreases.
      solver.Order = RungeKuttaSolver.SolverOrder.Second;
      tabulatedValues = solver.Solve( prob );
      calculatedSolutionYvalues = new DoubleVector( tabulatedValues.Value );
      maxDiff = ( calculatedSolutionYvalues - actualSolutionYvalues ).InfinityNorm();
      Console.WriteLine( "Maximum error using second order Runge Kutta = " + maxDiff );

      // You can also specify a TabulatedFunction object into which to place the
      // solution. We specify a NaturalCubicSpline function and use a fourth order
      // Runge Kutta to solve the same problem. Note the solution improves.
      TabulatedFunction solutionFunction = new NaturalCubicSpline();
      solver.Order = RungeKuttaSolver.SolverOrder.Fourth;
      solver.Solve( prob, ref solutionFunction );
      maxDiff = 0;
      for ( int i = 0; i < solutionFunction.NumberOfTabulatedValues; ++i )
      {
        double absDiff = Math.Abs( solutionFunction.GetY( i ) - actualSolutionFunction( solutionFunction.GetX( i ) ) );
        maxDiff = Math.Max( maxDiff, absDiff );
      }
      Console.WriteLine( "Maximum error using fourth order Runge Kutta = " + maxDiff );

      Console.WriteLine();
      Console.WriteLine("Press Enter Key");
      Console.Read();
    }
  }
}