# Differential Equations QuickStart Sample (C#)

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

View this sample in: Visual Basic F# IronPython

``````using System;

using Extreme.Mathematics;
using Extreme.Mathematics.Calculus.OrdinaryDifferentialEquations;

namespace Extreme.Numerics.QuickStart.CSharp
{
/// <summary>
/// Illustrates integrating systems of ordinary differential equations
/// (ODE's) using classes in the
/// Extreme.Mathematics.Calculus.OrdinaryDifferentialEquations
/// namespace.
/// </summary>
class DifferentialEquations
{
static void Main(string[] args)
{
// The license is verified at runtime. We're using
// https://www.extremeoptimization.com/trial-key

// 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;

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: {0}", cvode.CurrentTime);
// CurrentValue returns the corresponding value:
Console.WriteLine("Value at this t: {0:F14}", cvode.CurrentValue);
// 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: {0}", 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...");
}

/// <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>
private 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 = sigma * (y - y);
dy = y * (rho - y) - y;
dy = y * y - beta * y;

return dy;
}

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

dy = y * y * (1 - y);

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>
private 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) * y;

return J;
}
}
}``````