Higher Dimensional Numerical Integration in C# QuickStart Sample

Illustrates numerical integration of functions in higher dimensions using classes in the Numerics.NET.Calculus namespace in C#.

This sample is also available in: Visual Basic, F#, IronPython.

Overview

This QuickStart sample demonstrates how to perform numerical integration of functions in two dimensions using Numerics.NET’s integration capabilities.

The sample shows several key aspects of multi-dimensional integration:

• Using the AdaptiveIntegrator2D class for efficient integration over rectangular regions
• Working with the Repeated1DIntegrator2D class that uses repeated one-dimensional integration
• Setting up integrand functions using delegates
• Integrating over arbitrary non-rectangular regions by specifying boundary functions
• Monitoring integration progress through properties like Status, EstimatedError, and function evaluation counts
• Comparing numerical results with known analytical solutions

The example includes integration of several test functions including:

• A rational function: 4/(1 + 2x + 2y) over the unit square
• A trigonometric function: (π²/4)sin(πx)sin(πy) over the unit square
• A polynomial function: x²y² over the unit disk

For each integration, the sample demonstrates how to set up the problem, execute the integration, and analyze the results including error estimates and computational effort required.

The code

using System;
// The numerical integration classes reside in the
// Numerics.NET.Calculus namespace.
using Numerics.NET.Calculus;
// Function delegates reside in the Numerics.NET
// namespace.
using Numerics.NET;

// Illustrates numerical integration in higher dimensions using
// classes in the Numerics.NET.Calculus namespace of Numerics.NET.

// The license is verified at runtime. We're using
// a 30 day trial key here. For more information, see
//     https://numerics.net/trial-key
Numerics.NET.License.Verify("your-trial-key-here");

//
// Two-dimensional integration
//

// The function we are integrating must be
// provided as a Func<double, double, double> delegate.

// Variable to hold the result:
double result;

// The AdaptiveIntegrator2D class is the most efficient
// 2D integrator in most cases. It uses an adaptive algorithm.

// Construct an instance of the integrator class:
AdaptiveIntegrator2D integrator1 = new AdaptiveIntegrator2D();

// An example of setting the integrand and bounds through properties
// is given below. Here, we put the integrand and the bounds
// of the integration region directly in the call to Integrate,
// which performs the calculation:
Func<double, double, double> f1 = (x, y) => 4 / (1 + 2 * x + 2 * y);
integrator1.Integrate(f1, 0, 1, 0, 1);
Console.WriteLine("4 / (1 + 2x + 2y) on [0,1] * [0,1]");
Console.WriteLine($"  Value:       {integrator1.Result:F15}");
Console.WriteLine($"  Exact value: {Math.Log(3125.0 / 729.0):F15} = Ln(3125 / 729)");
// To see whether the algorithm ended normally,
// inspect the Status property:
Console.WriteLine("  Status: {0}",
    integrator1.Status);
Console.WriteLine("  Estimated error: {0}",
    integrator1.EstimatedError);
Console.WriteLine("  Iterations: {0}",
    integrator1.IterationsNeeded);
Console.WriteLine("  Function evaluations: {0}",
    integrator1.EvaluationsNeeded);

// Another integrator uses repeated 1-dimensional
// integration:
Repeated1DIntegrator2D integrator2 = new Repeated1DIntegrator2D();

// You can set the order of integration, as well as
// the integration rules for the X and the Y direction:
integrator2.InitialDirection = Repeated1DIntegratorDirection.X;

// You can set the integrand and the bounds of the integration region
// by setting properties of the integrator object:
integrator2.Integrand = (x, y) => Math.PI * Math.PI / 4 * Math.Sin(Math.PI * x) * Math.Sin(Math.PI * y);
integrator2.XLowerBound = 0;
integrator2.XUpperBound = 1;
integrator2.YLowerBound = 0;
integrator2.YUpperBound = 1;

result = integrator2.Integrate();
Console.WriteLine("Pi^2 / 4 Sin(Pi x) Sin(Pi y)   on [0,1] * [0,1]");
Console.WriteLine($"  Value:       {integrator2.Result:F15}");
Console.WriteLine($"  Exact value: {1.0:F15}");
// To see whether the algorithm ended normally,
// inspect the Status property:
Console.WriteLine("  Status: {0}",
    integrator2.Status);
Console.WriteLine("  Estimated error: {0}",
    integrator2.EstimatedError);
Console.WriteLine("  Iterations: {0}",
    integrator2.IterationsNeeded);
Console.WriteLine("  Function evaluations: {0}",
    integrator2.EvaluationsNeeded);

//
// Integration over arbitrary regions
//

// The repeated 1D integrator can also be used to compute
// integrals over arbitrary regions. To do this, you need to
// supply function that return the lower bound and upper bound
// of the region as a function of x.

// Here, we integrate x^2 * y^2 over the unit disk.
integrator2.LowerBoundFunction = x => Math.Abs(x) >= 1.0 ? 0.0 : -Math.Sqrt(1.0 - x*x);
integrator2.UpperBoundFunction = x => Math.Abs(x) >= 1.0 ? 0.0 : Math.Sqrt(1.0 - x*x);
integrator2.XLowerBound = -1;
integrator2.XUpperBound = 1;

integrator2.Integrand = (x, y) => x * x * y * y;

result = integrator2.Integrate();
Console.WriteLine("x^2 * y^2 on the unit disk");
Console.WriteLine($"  Value:       {integrator2.Result:F15}");
Console.WriteLine($"  Exact value: {Math.PI / 24:F15} = Pi / 24");
// To see whether the algorithm ended normally,
// inspect the Status property:
Console.WriteLine("  Status: {0}",
    integrator2.Status);
Console.WriteLine("  Estimated error: {0}",
    integrator2.EstimatedError);
Console.WriteLine("  Iterations: {0}",
    integrator2.IterationsNeeded);
Console.WriteLine("  Function evaluations: {0}",
    integrator2.EvaluationsNeeded);

Console.Write("Press Enter key to exit...");
Console.ReadLine();