Higher Dimensional Numerical Integration QuickStart Sample (C#)

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

Visual Basic code F# code IronPython code Back to QuickStart Samples

using System;

namespace Extreme.Mathematics.QuickStart.CSharp
{
	// The numerical integration classes reside in the
	// Extreme.Mathematics.Calculus namespace.
	using Extreme.Mathematics.Calculus;
	// Function delegates reside in the Extreme.Mathematics
	// namespace.
	using Extreme.Mathematics;

	/// <summary>
	/// Illustrates numerical integration in higher dimensions using
	/// classes in the Extreme.Mathematics.Calculus namespace of the Extreme
	/// Optimization Numerical Libraries for .NET.
	/// </summary>
	class NDIntegration
	{
		/// <summary>
		/// The main entry point for the application.
		/// </summary>
		[STAThread]
		static void Main(string[] args)
		{
			//
			// 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:       {0:F15}", integrator1.Result);
            Console.WriteLine("  Exact value: {0:F15} = Ln(3125 / 729)", Math.Log(3125.0 / 729.0));
            // 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:       {0:F15}", integrator2.Result);
            Console.WriteLine("  Exact value: {0:F15}", 1.0);
            // 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:       {0:F15}", integrator2.Result);
            Console.WriteLine("  Exact value: {0:F15} = Pi / 24", Math.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();
		}

    }
}