Evaluate an Integral

Here is some code to numerically evaluate a simple, one-dimensional integral:

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

// This is the area of the upper half-circle, so the value should be pi/2.
IntegrationResult i = FunctionMath.Integrate(x => Math.Sqrt(1.0 - x * x), -1.0, +1.0);
Console.WriteLine($"i = {i.Estimate} after {i.EvaluationCount} evaluations.");

Can I do multi-dimensional integrals?

Yes! Here is some code to evaluate one of the famous 3-dimensional Watson integrals:

Interval zeroToPi = Interval.FromEndpoints(0.0, Math.PI);
IntegrationResult watson = MultiFunctionMath.Integrate(
   x => 1.0 / (1.0 - Math.Cos(x[0]) * Math.Cos(x[1]) * Math.Cos(x[2])),
   new Interval[] { zeroToPi, zeroToPi, zeroToPi }
);

To do multidimensional integrals: use the MultiFunctionMath class instead of the FunctionMath class, provide a function that accepts an IReadOnlyList instead of just one double, and provide an IReadOnlyList to define the integration limits. The default precision targets for multi-dimensional integrals are lower than for one-dimensional integrals, and the default evaluation count limits are higher, but you can adjust them to your liking using IntegrationSettings.

How can I change the target precision?

Use the overloads that accept IntegrationSettings parameters and pass them a settings object that specifies your desired precision. Here is our half-circle integral evaluated to just 0.01% accuracy, which will require many fewer function evaluations:

IntegrationResult i = FunctionMath.Integrate(
    x => Math.Sqrt(1.0 - x * x), -1.0, +1.0,
    new IntegrationSettings() { RelativePrecision = 1.0E-4 }
);

You can specify the target precision either in absolute terms by setting the AbsolutePrecision property (e.g. AbsolutePrecision = 0.01 indicates that the uncertainty should be 0.01 or less) or in relative terms by settings the RelativePrecision property (e.g. RelativePrecision = 0.01 indicates that the uncertainty should be 1% or less).

If both AbsolutePrecision and RelativePrecision are specified, evaluation will terminate as soon as one of the criteria is satisfied. You can use this behavior to require a purely absolute or purely relative termination criteria by setting the other property value to 0.0, so that it will never be satisfied.

You can also change the maximum allowed number of evaluations by setting the EvaluationBudget property. If the required precision is not achieved after this number of function evaluations, the method will throw a NonconvergenceException.

Can you handle singular integrals?

Yes! Assuming, of course, that the singularity still integrates to a finite value. Soft singularities like logs are no problem:

// This integrand has a log singularity at x = 0.
// The value of this integral is the Catalan constant.
IntegrationResult soft = FunctionMath.Integrate(
    x => -Math.Log(x) / (1.0 + x * x), 0.0, 1.0
);

Harder singularities (power laws close to the 1/x limit of integrability) can also be integrated, but only as long as precision requirements aren't too high:

// This integral has a power law singularity at x = 0.
// The value of this integral is 4.
IntegrationResult hard = FunctionMath.Integrate(
    x => Math.Pow(x, -3.0 / 4.0), 0.0, 1.0,
    new IntegrationSettings() { RelativePrecision = 1.0E-6 }
);

Can you handle integrals over infinite ranges?

Yes! Just specify Double.PositiveInfinity (and/or Double.NegativeInfinity) as an integration limit.

// The value of this infinite integral is sqrt(pi)
IntegrationResult infinite = FunctionMath.Integrate(
    x => Math.Exp(-x* x), Double.NegativeInfinity, Double.PositiveInfinity
);

Can I monitor the progress of integration?

Yes! In the IntegrationSettings object, set the Listener property to an Action delegate (a function that accepts one IntegrationResult argument and has no return value). Your listener will be called periodically as the integral is evaluated, allowing you to monitor the number of evaluations and precision achieved before the termination criteria are met.

Here is a monster-example which pulls it all together to monitor progress toward the numerical evaluation of a Box integral: the six-dimensional integral that measures the average distance between two points in the unit cube:

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

Func<IReadOnlyList<double>,double> distance = z => {
    ColumnVector x = new ColumnVector(z[0], z[1], z[2]);
    ColumnVector y = new ColumnVector(z[3], z[4], z[5]);
    ColumnVector d = x - y;
    return(d.Norm());
};

Interval oneBox = Interval.FromEndpoints(0.0, 1.0);
Interval[] sixBox = new Interval[] { oneBox, oneBox, oneBox, oneBox, oneBox, oneBox };

IntegrationSettings settings= new IntegrationSettings() {
    RelativePrecision = 1.0E-4,
    AbsolutePrecision = 0.0,
    Listener = r => {
        Console.WriteLine($"Estimate {r.Estimate} after {r.EvaluationCount} evaluations.");
    }
};

IntegrationResult numeric = MultiFunctionMath.Integrate(distance, sixBox, settings);
Console.WriteLine($"The numeric result is {numeric.Estimate}.");

double analytic = 4.0 / 105.0 + 17.0 / 105.0 * Math.Sqrt(2.0) - 2.0 / 35.0 * Math.Sqrt(3.0)
    + Math.Log(1.0 + Math.Sqrt(2.0)) / 5.0 + 2.0 / 5.0 * Math.Log(2.0 + Math.Sqrt(3.0))
    - Math.PI / 15.0;
Console.WriteLine($"The analytic result is {analytic}.");