Numerical Integration in Two Dimensions

The NumericalIntegrator2D 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 FuncT1, T2, TResult taking two numbers as arguments and returning a real number. It can be set through the Integrand property. The boundaries of the integration region are determined by the XLowerBound, XUpperBound, YLowerBound, and YUpperBound 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, XLowerBound, XUpperBound, YLowerBound, and YUpperBound 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.

Two algorithms for 2D integration are available. The AdaptiveIntegrator2D 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 square:

private static double Integrand1(double x, double y) 
{ return 4 / (1 + 2*x + 2*y); }

BivariateRealFunction f1 = new BivariateRealFunction(Integrand1);
AdaptiveIntegrator2D integrator1 = new AdaptiveIntegrator2D();
double result = integrator1.Integrate(f1, 0, 1, 0, 1);

Function Integrand1(ByVal x As Double, ByVal y As Double) As Double
    Return 4 / (1 + 2 * x + 2 * y)
End Function

Dim f1 As New BivariateRealFunction(AddressOf Integrand1)
Dim integrator1 As New AdaptiveIntegrator2D()
Dim result As Double = integrator1.Integrate(f1, 0, 1, 0, 1)

The Repeated1DIntegrator2D class computes the integral by repeated integration. The XIntegrationRule and YIntegrationRule properties specify the 1D integrator that is to be used to perform the integration in the X or the Y direction, respectively. They are of type NumericalIntegrator. The InitialDirection property specifies in which direction the integration should be performed first. It is of type Repeated1DIntegratorDirection. By default, integration is done first in the X direction. The following code illustrates the above.

private static double Integrand2(double x, double y) 
{ return Math.PI * Math.PI / 4 * Math.Sin(Math.PI * x) * Math.Sin(Math.PI * y) }

BivariateRealFunction f2 = new BivariateRealFunction(Integrand2);
Repeated1DIntegrator2D integrator2 = new Repeated1DIntegrator2D();

integrator2.InitialDirection = Repeated1DIntegratorDirection.X;
integrator2.Integrand = f2;
integrator2.XLowerBound = 0;
integrator2.XUpperBound = 1;
integrator2.YLowerBound = 0;
integrator2.YUpperBound = 1;

result = integrator2.Integrate();

Function Integrand2(ByVal x As Double, ByVal y As Double) As Double
    Return Math.PI * Math.PI / 4 * Math.Sin(Math.PI * x) * Math.Sin(Math.PI * y)
End Function

Dim f2 As New BivariateRealFunction(AddressOf Integrand2)
Dim integrator2 As New Repeated1DIntegrator2D()
integrator2.Integrand = f2
integrator2.XLowerBound = 0.0
integrator2.XUpperBound = 1.0
integrator2.YLowerBound = 0.0
integrator2.YUpperBound = 1.0
result = integrator2.Integrate()

The Repeated1DIntegrator2D class can also perform integration over non-rectangular regions. Specifically, it can integrate over regions whose lower and upper bound in one dimension can be expressed as a function of the coordinate in the other dimension. A large class of regions can be defined in this way, including disks and triangles. The lower and upper bounds are set using the LowerBoundFunction and UpperBoundFunction properties. These are FuncT, TResult delegates. These bounds always apply to the initial integration direction. The following code shows how to integrate a function over a disk shaped region.

private static double Integrand1(double x, double y) 
{ return x * x * y * y; }
private static double DiskLowerBound(double x)
{
    x = Math.Abs(x);
    return (x >= 1) ? 0 : -Math.Sqrt(1 - x * x);
}
private static double DiskUpperBound(double x)
{
    x = Math.Abs(x);
    return (x >= 1) ? 0 : +Math.Sqrt(1 - x * x);
}

integrator2.Integrand = new BivariateRealFunction(Integrand3);
integrator2.LowerBoundFunction = new RealFunction(DiskLowerBound);
integrator2.UpperBoundFunction = new RealFunction(DiskUpperBound);
integrator2.XLowerBound = -1;
integrator2.XUpperBound = +1;

result = integrator2.Integrate();

Function DiskLowerBound(ByVal x As Double) As Double
    x = Math.Abs(x)
    If (x >= 1) Then
        Return 0
    End If
    Return -Math.Sqrt(1 - x * x)
End Function
Function DiskUpperBound(ByVal x As Double) As Double
    x = Math.Abs(x)
    If (x >= 1) Then
        Return 0
    End If
    Return Math.Sqrt(1 - x * x)
End Function
Function Integrand3(ByVal x As Double, ByVal y As Double) As Double
    Return x * x * y * y
End Function

integrator2.Integrand = f3
integrator2.LowerBoundFunction = AddressOf DiskLowerBound
integrator2.UpperBoundFunction = AddressOf DiskUpperBound
integrator2.XLowerBound = -1.0
integrator2.XUpperBound = 1.0

result = integrator2.Integrate()

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.