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#.

View this sample in: Visual Basic F# IronPython

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;

namespace Numerics.NET.QuickStart.CSharp
{
    /// <summary>
    /// Illustrates numerical integration in higher dimensions using
    /// classes in the Numerics.NET.Calculus namespace of Numerics.NET.
    /// </summary>
    class NDIntegration
    {
        static void Main(string[] args)
        {
            // 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("64542-18980-57619-62268");

            //
            // 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();
        }

    }
}