Differential Equations in C# QuickStart Sample

Illustrates integrating systems of ordinary differential equations (ODE’s) in C#.

View this sample in: Visual Basic F# IronPython

using System;

using Numerics.NET;
using Numerics.NET.Calculus.OrdinaryDifferentialEquations;

namespace Numerics.NET.QuickStart.CSharp
{
    /// <summary>
    /// Illustrates integrating systems of ordinary differential equations
    /// (ODE's) using classes in the 
    /// Numerics.NET.Calculus.OrdinaryDifferentialEquations
    /// namespace.
    /// </summary>
    class DifferentialEquations
    {
        static void Main(string[] args)
        {
            // The license is verified at runtime. We're using
            // a 30 day trial key here. For more information, see
            //     https://numerics.net/trial-key
            Numerics.NET.License.Verify("64542-18980-57619-62268");

            // The ClassicRungeKuttaIntegrator class implements the
            // well-known 4th order fixed step Runge-Kutta method.
            ClassicRungeKuttaIntegrator rk4 = new ClassicRungeKuttaIntegrator();

            // The differential equation is expressed in terms of a 
            // DifferentialFunction delegate. This is a function that
            // takes a double (time value) and two Vectors (y value and
            // return value)  as arguments.
            //
            // The Lorentz function below defines the differential function
            // for the Lorentz attractor.
            rk4.DifferentialFunction = Lorentz;

            // To perform the computations, we need to set the initial time...
            rk4.InitialTime = 0.0;
            // and the initial value.
            rk4.InitialValue = Vector.Create(1.0, 0.0, 0.0);
            // The Runge-Kutta integrator also requires a step size:
            rk4.InitialStepsize = 0.1;

            Console.WriteLine("Classic 4th order Runge-Kutta");
            for (int i = 0; i <= 5; i++) {
                double t = 0.2 * i;
                // The Integrate method performs the integration.
                // It takes as many steps as necessary up to
                // the specified time and returns the result:
                var y = rk4.Integrate(t);
                // The IterationsNeeded always shows the number of steps
                // needed to arrive at the final time.
                Console.WriteLine("{0:F2}: {1,20:F14} ({2} steps)", t, y, rk4.IterationsNeeded);
            }

            // The RungeKuttaFehlbergIntegrator class implements a variable-step
            // Runge-Kutta method. This method chooses the stepsize, and so
            // is generally more reliable.
            RungeKuttaFehlbergIntegrator rkf45 = new RungeKuttaFehlbergIntegrator();

            rkf45.InitialTime = 0.0;
            rkf45.InitialValue = Vector.Create(1.0, 0.0, 0.0);
            rkf45.DifferentialFunction = Lorentz;
            rkf45.AbsoluteTolerance = 1e-8;

            Console.WriteLine("Classic 4/5th order Runge-Kutta-Fehlberg");
            for (int i = 0; i <= 5; i++) {
                double t = 0.2 * i;
                var y = rkf45.Integrate(t);
                Console.WriteLine("{0:F2}: {1,20:F14} ({2} steps)", t, y, rkf45.IterationsNeeded);
            }

            // The CVODE integrator, part of the SUNDIALS suite of ODE solvers,
            // is the most advanced of the ODE integrators.
            CvodeIntegrator cvode = new CvodeIntegrator();

            cvode.InitialTime = 0.0;
            cvode.InitialValue = Vector.Create(1.0, 0.0, 0.0);
            cvode.DifferentialFunction = Lorentz;
            cvode.AbsoluteTolerance = 1e-8;

            Console.WriteLine("CVODE (multistep Adams-Moulton)");
            for (int i = 0; i <= 5; i++) {
                double t = 0.2 * i;
                var y = cvode.Integrate(t);
                Console.WriteLine("{0:F2}: {1,20:F14} ({2} steps)", t, y, cvode.IterationsNeeded);
            }

            //
            // Other properties and methods
            //

            // The IntegrateSingleStep method takes just a single step
            // in the direction of the specified final time:
            cvode.IntegrateSingleStep(2.0);
            // The CurrentTime property returns the corresponding time value.
            Console.WriteLine($"t after single step: {cvode.CurrentTime}");
            // CurrentValue returns the corresponding value:
            Console.WriteLine($"Value at this t: {cvode.CurrentValue:F14}");
            // The IntegrateMultipleSteps method performs the integration
            // until at least the final time, and returns the last
            // value that was computed:
            cvode.IntegrateMultipleSteps(2.0);
            Console.WriteLine($"t after multiple steps: {cvode.CurrentTime}");

            //
            // Stiff systems
            //

            Console.WriteLine("\nStiff systems");

            // The CVODE integrator is the only ODE integrator that can
            // handle stiff problems well. The following example uses
            // an equation for the size of a flame. The smaller 
            // the initial size, the more stiff the equation is.

            // We compare the performance of the variable step Runge-Kutta method
            // and the CVODE integrator:

            double delta = 0.0001;
            double tFinal = 2 / delta;

            rkf45 = new RungeKuttaFehlbergIntegrator();
            rkf45.InitialTime = 0.0;
            rkf45.InitialValue = Vector.Create(delta);
            rkf45.DifferentialFunction = Flame;

            Console.WriteLine("Classic 4/5th order Runge-Kutta-Fehlberg");
            for (int i = 0; i <= 10; i++) {
                double t = i * 0.1 * tFinal;
                var y = rkf45.Integrate(t);
                Console.WriteLine("{0:F2}: {1,20:F14} ({2} steps)", t, y, rkf45.IterationsNeeded);
            }

            // The CVODE integrator will use a special (implicit) method
            // for stiff problems. To select this method, pass 
            // EdeKind.Stiff as an argument to the constructor.
            cvode = new CvodeIntegrator(OdeKind.Stiff);
            cvode.InitialTime = 0.0;
            cvode.InitialValue = Vector.Create(delta);
            cvode.DifferentialFunction = Flame;
            // When solving stiff problems, a Jacobian for the system
            // must be supplied. See below for the definition.
            cvode.SetDenseJacobian(FlameJacobian);

            // Notice how the CVODE integrator takes a lot fewer steps,
            // and is also more accurate than the variable-step
            // Runge-Kutta method.
            Console.WriteLine("CVODE (Stiff - Backward Differentiation Formula)");
            for (int i = 0; i <= 10; i++) {
                double t = i * 0.1 * tFinal;
                var y = cvode.Integrate(t);
                Console.WriteLine("{0:F2}: {1,20:F14} ({2} steps)", t, y, cvode.IterationsNeeded);
            }

            Console.Write("Press Enter key to exit...");
            Console.ReadLine();
        }
      
        /// <summary>
        /// Represents the differential function for the Lorentz attractor.
        /// </summary>
        /// <param name="t">The time value.</param>
        /// <param name="y">The current value.</param>
        /// <param name="dy">On output, the first derivatives.</param>
        /// <returns>A reference to <paramref name="dy"/>.</returns>
        /// <remarks><paramref name="dy"/> may be <see langword="null"/>
        /// on input.</remarks>
        static Vector<double> Lorentz(double t, Vector<double> y, Vector<double> dy) 
        {
            if (dy == null)
                dy = Vector.Create<double>(3);

            double sigma = 10.0;
            double beta = 8.0 / 3.0;
            double rho = 28.0;

            dy[0] = sigma * (y[1] - y[0]);
            dy[1] = y[0] * (rho - y[2]) - y[1];
            dy[2] = y[0] * y[1] - beta * y[2];

            return dy;
        }

        /// <summary>
        /// Represents the differential function for the flame expansion
        /// problem.
        /// </summary>
        static Vector<double> Flame(double t, Vector<double> y, Vector<double> dy) 
        {
            if (dy == null)
                dy = Vector.Create<double>(1);

            dy[0] = y[0] * y[0] * (1 - y[0]);

            return dy;
        }
       
        /// <summary>
        /// Represents the Jacobian of the differential function 
        /// for the flame expansion problem.
        /// </summary>
        /// <param name="t">The time value.</param>
        /// <param name="y">The current value.</param>
        /// <param name="dy">The current value of the first derivatives.</param>
        /// <param name="J">A Matrix that, on output, contains the value
        /// of the Jacobian.</param>
        /// <returns>A reference to <paramref name="J"/>.</returns>
        /// <remarks>The Jacobian is the matrix of partial derivatives of each
        /// equation in the system with respect to each variable in the system.
        /// <paramref name="J"/> may be <see langword="null"/>
        /// on input.</remarks>
        static Matrix<double> FlameJacobian(double t, Vector<double> y, Vector<double> dy, Matrix<double> J) {
            
            if (J == null)
                J = Matrix.Create<double>(1, 1);

            J[0, 0] = (2 - 3 * y[0]) * y[0];

            return J;
        }
    }
}