Basic Integration QuickStart Sample (C#)

Illustrates the basic numerical integration (quadrature) classes in C#.

VB.NET code F# 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 the basic use of the numerical integration
    /// classes in the Extreme.Mathematics.Calculus namespace.
    /// </summary>
    class BasicIntegration
    {
        /// <summary>
        /// The main entry point for the application.
        /// </summary>
        [STAThread]
        static void Main(string[] args)
        {
            // Numerical integration algorithms fall into two
            // main categories: adaptive and non-adaptive.
            // This QuickStart Sample illustrates the use of
            // the non-adaptive numerical integrators.
            //
            // All numerical integration classes derive from
            // NumericalIntegrator. This abstract base class
            // defines properties and methods that are shared
            // by all numerical integration classes.

            //
            // The integrand
            //
            // The function we are integrating must be
            // provided as a RealFunction. For more
            // information about this delegate, see the
            // FunctionDelegates QuickStart sample.
            RealFunction f = 
                new RealFunction(Math.Sin);
            // Variable to hold the result:
            double result;
            
            //
            // SimpsonIntegrator
            // 

            // The simplest numerical integration algorithm
            // is Simpson's rule. 
            SimpsonIntegrator simpson = new SimpsonIntegrator();
            // You can set the relative or absolute tolerance
            // to which to evaluate the integral.
            simpson.RelativeTolerance = 1e-5;
            // You can select the type of tolerance using the
            // ConvergenceCriterion property:
            simpson.ConvergenceCriterion =
                ConvergenceCriterion.WithinRelativeTolerance;
            // The Integrate method performs the actual 
            // integration:
            result = simpson.Integrate(f, 0, 2);
            Console.WriteLine("sin(x) on [0,2]");
            Console.WriteLine("Simpson integrator:");
            // The result is also available in the Result 
            // property:
            Console.WriteLine("  Value: {0}", simpson.Result);
            // To see whether the algorithm ended normally,
            // inspect the Status property:
            Console.WriteLine("  Status: {0}",
                simpson.Status);
            // You can find out the estimated error of the result
            // through the EstimatedError property:
            Console.WriteLine("  Estimated error: {0}", 
                simpson.EstimatedError);
            // The number of iterations to achieve the result
            // is available through the IterationsNeeded property.
            Console.WriteLine("  Iterations: {0}",
                simpson.IterationsNeeded);
            // The number of function evaluations is available 
            // through the FunctionEvaluationsNeeded property.
            Console.WriteLine("  Function evaluations: {0}",
                simpson.FunctionEvaluationsNeeded);

            //
            // Gauss-Kronrod Integration
            //

            // Gauss-Kronrod integrators also use a fixed point 
            // scheme, but with certain optimizations in the 
            // choice of points where the integrand is evaluated.
            //
            // The Extreme.Mathematics.Calculus namespace contains a series
            // of Gauss-Kronrod integrators that are mainly used
            // as an integration rule for the adaptive 
            // integrator. No iteration is used.
            // Here's the 21-point rule:
            GaussKronrodIntegrator21 gk21 = new GaussKronrodIntegrator21();
            gk21.Integrate(f, 0, 2);
            Console.WriteLine("21 point Gauss-Kronrod rule:");
            Console.WriteLine("  Value: {0}", gk21.Result);
            Console.WriteLine("  Status: {0}",
                gk21.Status);
            Console.WriteLine("  Estimated error: {0}", 
                gk21.EstimatedError);
            Console.WriteLine("  Iterations: {0}",
                gk21.IterationsNeeded);
            Console.WriteLine("  Function evaluations: {0}",
                gk21.FunctionEvaluationsNeeded);

            // The NonAdaptiveGaussKronrodIntegrator uses a
            // succession of 10, 21, 43, and 87 point rules
            // to approximate the integral.
            NonAdaptiveGaussKronrodIntegrator nagk =
                new NonAdaptiveGaussKronrodIntegrator();
            nagk.Integrate(f, 0, 2);
            Console.WriteLine("Non-adaptive Gauss-Kronrod rule:");
            Console.WriteLine("  Value: {0}", nagk.Result);
            Console.WriteLine("  Status: {0}",
                nagk.Status);
            Console.WriteLine("  Estimated error: {0}", 
                nagk.EstimatedError);
            Console.WriteLine("  Iterations: {0}",
                nagk.IterationsNeeded);
            Console.WriteLine("  Function evaluations: {0}",
                nagk.FunctionEvaluationsNeeded);
            
            //
            // Romberg Integration
            //
            // Romberg integration combines Simpson's Rule
            // with a scheme to accelerate convergence.
            // This algorithm is useful for smooth integrands.
            RombergIntegrator romberg = new RombergIntegrator();
            result = romberg.Integrate(f, 0, 2);
            Console.WriteLine("Romberg integration:");
            Console.WriteLine("  Value: {0}", romberg.Result);
            Console.WriteLine("  Status: {0}",
                romberg.Status);
            Console.WriteLine("  Estimated error: {0}", 
                romberg.EstimatedError);
            Console.WriteLine("  Iterations: {0}",
                romberg.IterationsNeeded);
            Console.WriteLine("  Function evaluations: {0}",
                romberg.FunctionEvaluationsNeeded);
            
            // However, it breaks down if the integration
            // algorithm contains singularities or 
            // discontinuities.
            f = new RealFunction(HardIntegrand);
            result = romberg.Integrate(f, 0, 1);
            Console.WriteLine("Romberg on hard integrand:");
            Console.WriteLine("  Value: {0}", romberg.Result);
            Console.WriteLine("  Actual value: 100");
            Console.WriteLine("  Status: {0}",
                romberg.Status);
            Console.WriteLine("  Estimated error: {0}", 
                romberg.EstimatedError);
            Console.WriteLine("  Iterations: {0}",
                romberg.IterationsNeeded);
            Console.WriteLine("  Function evaluations: {0}",
                romberg.FunctionEvaluationsNeeded);

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

        /// <summary>
        /// Function that will cause difficulties to the
        /// simplistic integration algorithms.
        /// </summary>
        private static double HardIntegrand(double x)
        {
            // This is put in because some integration rules
            // evaluate the function at x=0.
            if (x <= 0)
                return 0;
            return Math.Pow(x,-0.9) * Math.Log(1/x);
        }
    }
}