Numerical Integration in more than Two Dimensions

The NumericalIntegratorND class is the abstract base class for all classes that implement numerical integration in two dimensions. It inherits from IterativeAlgorithm. The AbsoluteTolerance and RelativeTolerance properties set the desired precision as specified by the ConvergenceCriterion property. The default value for both tolerances is SqrtEpsilon (roughly 10^-8). MaxIterations sets the maximum number of iterations. The default value for this property depends on the algorithm used. IterationsNeeded returns the actual number of iterations performed after the algorithm has completed.

Specific to this class are the Order and EvaluationsNeeded properties, and the overloaded Integrate method.

The Order property gives the order of the integration algorithm. The order of an integration algorithm is the highest degree of a general polynomial (or multinomial in this case) whose integral is calculated exactly by the algorithm. A method of order three integrates cubic polynomials exactly. Many methods have a fixed order. For some algorithms, the order depends on the input values.

The EvaluationsNeeded property returns the total number of times the target function was evaluated while approximating the integral. This property is a more reliable indication of the efficiency of an algorithm than IterationsNeeded. For some algorithms, the number of function evaluations grows exponentially with each iteration, while for others it is a simple multiple. Even though higher order methods are slower, they usually require less subdivisions of the integration interval, which makes them more desirable for smooth target functions. For target functions with integrable singularities, a low order method is usually preferred.

The integrand must be specified as a FuncT, TResult that maps a VectorT to a real number. It can be set through the Integrand property. The boundaries of the integration region are determined by the LowerBounds and UpperBounds properties.

The Integrate method does the actual work of numerically integrating a target function. It has three overloads. With no parameters, the method uses the values supplied through the Integrand, LowerBounds and UpperBounds properties.

The remaining two overloads take four or five parameters. The first parameter, if present, is a FuncT, TResult delegate that specifies the target function. The remaining four parameters are Double values that specify the lower and upper bounds of the integration region in the X direction followed by the bounds in the Y direction.

The AdaptiveIntegratorND class uses an adaptive method. It computes an approximation of the integral over a region, as well as an estimate for the error. The region is divided into two new regions. This process is repeated recursively until the estimated error in each region is small enough. The adaptive integrator can only integrate over rectangular regions. The following code shows how to integrate a function over the unit cube:

Func<Vector<double>, double> fN = x => Elementary.Sinc(x.Norm());
NumericalIntegratorND integrator = new AdaptiveIntegratorND();
Vector<double> lowerBounds = Vector.CreateConstant(3, 0.0);
Vector<double> upperBounds = Vector.CreateConstant(3, 1.0);
integrator.Integrate(fN, lowerBounds, upperBounds);
Console.WriteLine("  Value: {0}", integrator.Result);
Console.WriteLine("  Status: {0}", integrator.Status);
Console.WriteLine("  Estimated error: {0}", integrator.EstimatedError);
Console.WriteLine("  Iterations: {0}", integrator.IterationsNeeded);
Console.WriteLine("  Function evaluations: {0}", integrator.EvaluationsNeeded);

Dim fN As Func(Of Vector(Of Double), Double) =
    Function(x) Elementary.Sinc(x.Norm())
Dim integrator As NumericalIntegratorND = New AdaptiveIntegratorND()
Dim lowerBounds As Vector(Of Double) = Vector.CreateConstant(3, 0.0)
Dim upperBounds As Vector(Of Double) = Vector.CreateConstant(3, 1.0)
integrator.Integrate(fN, lowerBounds, upperBounds)
Console.WriteLine("  Value: {0}", integrator.Result)
Console.WriteLine("  Status: {0}", integrator.Status)
Console.WriteLine("  Estimated error: {0}", integrator.EstimatedError)
Console.WriteLine("  Iterations: {0}", integrator.IterationsNeeded)
Console.WriteLine("  Function evaluations: {0}", integrator.EvaluationsNeeded)

let fN = fun (x : Vector<double>) -> Elementary.Sinc(x.Norm())
let integrator = new AdaptiveIntegratorND()
let lowerBounds = Vector.CreateConstant(3, 0.0)
let upperBounds = Vector.CreateConstant(3, 1.0)
let result = integrator.Integrate(fN, lowerBounds, upperBounds)
Console.WriteLine("  Value: {0}", integrator.Result)
Console.WriteLine("  Status: {0}", integrator.Status)
Console.WriteLine("  Estimated error: {0}", integrator.EstimatedError)
Console.WriteLine("  Iterations: {0}", integrator.IterationsNeeded)
Console.WriteLine("  Function evaluations: {0}", integrator.EvaluationsNeeded)

The Integrate method does method always returns the best estimate for the integral. Successive calls to the Result property will also return this value, until the next call to Integrate.

If the ThrowExceptionOnFailure property is set to true, an exception is thrown if the algorithm has failed to obtain the integral with the desired accuracy. If false, the Integrate method returns the best approximation to the integral, regardless of whether it is within the requested tolerance.

The Status property indicates how the algorithm terminated. Its possible values and their meaning are listed below.