Integrate a Differential Equation

How do I integrate an ordinary differential equation?

Here is code to numerically integrate the differential equation y' = - x y:

using System;
using Meta.Numerics;
using Meta.Numerics.Analysis;

Func<double, double, double> rhs = (x, y) => - x * y;
OdeResult sln = FunctionMath.IntegrateOde(rhs, 0.0, 1.0, 2.0);
Console.WriteLine($"Numeric solution y({sln.X}) = {sln.Y}.");
Console.WriteLine($"Required {sln.EvaluationCount} evaluations.");
Console.WriteLine($"Analytic solution y({sln.X}) = {Math.Exp(-MoreMath.Sqr(sln.X) / 2.0)}");

It's easy to derive an analytic solution to this equation, so we show the analytic result for comparison.

How do I integrate coupled differential equations?

Let's solve the Lotka-Volterra equations, which model interdependent populations of predators and prey:

using System.Collections.Generic;

double A = 0.1;
double B = 0.02;
double C = 0.4;
double D = 0.02;
Func<double, IReadOnlyList<double>, IReadOnlyList<double>> lkRhs = (t, y) => {
    return(new double[] {
        A * y[0] - B * y[0] * y[1], D * y[0] * y[1] - C * y[1]
    });
};
MultiOdeSettings lkSettings = new MultiOdeSettings() {
    Listener = r =>  {Console.WriteLine($"t={r.X} rabbits={r.Y[0]}, foxes={r.Y[1]}"); }
};
MultiFunctionMath.IntegrateOde(lkRhs, 0.0, new double[] {20.0, 10.0}, 50.0, lkSettings);

Notice that we have used a settings object to add a listener that records values at intermediate times.

What about higher order equations?

The traditional way to handle a higher order differential equation is to convert it to a set of coupled first order differential equations. Consider the second order wave equation y'' = -y. Define y0 = y and y1 = y'. The original second order wave equation is then equivalent to the set of coupled first order equations y0' = y1, y1' = - y0. Here's some code that solves it to find sin(500):

Func<double, IReadOnlyList<double>, IReadOnlyList<double>> rhs1 = (x, u) => {
    return new double[] { u[1], -u[0] };
};
MultiOdeSettings settings1 = new MultiOdeSettings() { EvaluationBudget = 100000 };
MultiOdeResult result1 = MultiFunctionMath.IntegrateOde(
    rhs1, 0.0, new double[] { 0.0, 1.0}, 500.0, settings1
);
double s1 = MoreMath.Sqr(result1.Y[0]) + MoreMath.Sqr(result1.Y[1]);
Console.WriteLine($"y({result1.X}) = {result1.Y[0]}, (y)^2 + (y')^2 = {s1}");
Console.WriteLine($"Required {result1.EvaluationCount} evaluations.");

Notice that we have computed sin^2 + cos^2, which should always be 1, as a check on our result.

What if my equation is conservative?

Conservative differential equations are differential equations of the form y'' = f(x, y), where the right-hand-side does not have any explicit dependence on y'. The wave equation from the example above is conservative. Conservative differential equations can be solved more efficiently and with better respect for conserved quantities using a specialized algorithm that Meta.Numerics implements. Here is code that uses our conservative equation API to solve the wave equation problem from the example above:

Func<double, double, double> rhs2 = (x, y) => -y;
OdeSettings settings2 = new OdeSettings() { EvaluationBudget = 100000 };
OdeResult result2 = FunctionMath.IntegrateConservativeOde(
    rhs2, 0.0, 0.0, 1.0, 500.0, settings2
);
double s2 = MoreMath.Sqr(result2.Y) + MoreMath.Sqr(result2.YPrime);
Console.WriteLine($"y({result2.X}) = {result2.Y}, (y)^2 + (y')^2 = {s2}");
Console.WriteLine($"Required {result2.EvaluationCount} evaluations");

Notice that this approach requires significantly fewer right-hand-side evaluations and does a better job respecting sin^2 + cos^2 = 1.