**35.4****
****Stiff Equations** (.NET, C#, CSharp, VB, Visual Basic, F#)

A *stiff*
equation is a differential equation for which common numerical methods
for solving the equation are numerically unstable, unless the step size
is taken to be extremely small. In **NMath**,
class **VariableOrderOdeSolver**
solves stiff and non-stiff ordinary differential equations. The algorithm
uses higher order methods and smaller step size when the solution varies
rapidly.

The Solve() method solves the given initial value problem of ordinary differential equation of the form

or

for problems that involve a mass matrix *M*.

The function takes

● A delegate which evaluates the right hand side of the differential equations.

● A timespan vector specifying the interval of integration [t0, tf]. The solver imposes initial conditions at t0 and integrates from t0 to tf. If the timespan vector contains two elements, the solver returns the solution evaluated at every integration step. If the timespan vector contains more than two points, the solver returns the solution evaluated at the given points. The time values must be in order, either increasing or decreasing.

● The initial
value for problem—the value of the unknown function *y* at the initial time value.

For example, the van der Pol equations in relaxation
oscillation provide an example of a stiff system.^{1} The limit cycle has portions where the solution components change slowly
and the problem is quite stiff, alternating with regions of very sharp
change where it is not stiff.

To simulate this system, first create a function containing
the equations. The following code creates the parameterized Van der Pol
equation with parameter . Parameter *t* is the time value,
and *y* is the state value.

Code Example – C# stiff ODE solver

public static DoubleVector vdp1000( double t, DoubleVector y)

{

DoubleVector dy = new DoubleVector( 2 );

double y0 = y[0];

double y1 = y[1];

dy[0] = y1;

dy[1] = 1000 * ( 1 - (y0 * y0) ) * y1 - y0;

return dy;

}

Next, create a function that returns the Jacobian at
given *t *and *y*
values.

Code Example – C# stiff ODE solver

public static DoubleMatrix vdp1000Jac(double t, DoubleVector y)

{

DoubleMatrix J = new DoubleMatrix( 2, 2 );

J[0, 0] = 0;

J[0, 1] = 1;

J[1, 0] = -2 * 1000 * y[0] * y[1] - 1;

J[1, 1] = 1000 * ( 1.0 - y[0] * y[0] );

return J;

}

Create the initial values and time interval, and encapsulate the ODE function.

Code Example – C# stiff ODE solver

var y0 = new DoubleVector( 2.0, 0.0 );

var timespan = new DoubleVector( 0.0, 3000.0 );

var odeFunc =

new Func<double, DoubleVector, DoubleVector>( vdp1000 );

Create the solver object and set up the solver options. Here we use the default relative and absolute tolerances (1e-3 and 1e-6, respectively). Also, since we have an explicit form for the Jacobian function, we set this in the solver options too.

Code Example – C# stiff ODE solver

var ode15s = new VariableOrderOdeSolver();

var options = new VariableOrderOdeSolver.Options();

options.JacobianFunction =

new Func<double, DoubleVector, DoubleMatrix>( vdp1000Jac );

Finally, solve the equation and display the solution.

Code Example – C# stiff ODE solver

VariabeOrderOdeSolver.Solution<DoubleMatrix> soln =

ode15s.Solve( vdp1000, timespan, y0, options );

Console.WriteLine( "t = " + NMathFunctions.Round(soln.T, 4) );

Console.WriteLine();

Console.WriteLine( "y = " );

Console.WriteLine(

NMathFunctions.Round( soln.Y, 4 ).ToTabDelimited() );