Numerical Integration

The following double precision numerical integration or quadrature rules are supported in Math.NET Numerics under the MathNet.Numerics.Integration namespace. Unless stated otherwise, the examples below evaluate the integral $$\int_0^{10} x^2 \, dx = \frac{1000}{3} \approx 333.\overline{3}$$.

Simpson's Rule

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:  // Composite approximation with 4 partitions double composite = SimpsonRule.IntegrateComposite(x => x * x, 0.0, 10.0, 4); // Approximate value using IntegrateComposite with 4 partitions is: 333.33333333333337 Console.WriteLine("Approximate value using IntegrateComposite with 4 partitions is: " + composite); // Three point approximation double threePoint = SimpsonRule.IntegrateThreePoint(x => x * x, 0.0, 10.0); // Approximate value using IntegrateThreePoint is: 333.333333333333 Console.WriteLine("Approximate value using IntegrateThreePoint is: " + threePoint); 

Newton Cotes Trapezium Rule

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17:  // Adaptive approximation with a relative error of 1e-5 double adaptive = NewtonCotesTrapeziumRule.IntegrateAdaptive(x => x * x, 0.0, 10.0, 1e-5); // Approximate value of the integral using IntegrateAdaptive with a relative error of 1e-5 is: 333.333969116211 Console.WriteLine("Approximate value using IntegrateAdaptive with a relative error of 1e-5: " + adaptive); // Composite approximation with 15 partitions double composite = NewtonCotesTrapeziumRule.IntegrateComposite(x => x * x, 0.0, 10.0, 15); //Approximate value of the integral using IntegrateComposite with 15 partitions is: 334.074074074074 Console.WriteLine("Approximate value using IntegrateComposite with 15 partitions is: " + composite); // Two point approximation double twoPoint = NewtonCotesTrapeziumRule.IntegrateTwoPoint(x => x * x, 0.0, 10.0); //Approximate value using IntegrateTwoPoint is: 500 Console.WriteLine("Approximate value using IntegrateTwoPoint is: " + twoPoint); 

Double-Exponential Transformation

The Double-Exponential Transformation is suited for integration of smooth functions with no discontinuities, derivative discontinuities, and poles inside the interval.

 1: 2: 3: 4: 5:  // Approximate using a relative error of 1e-5. double integrate = DoubleExponentialTransformation.Integrate(x => x * x, 0.0, 10.0, 1e-5); // Approximate value using a relative error of 1e-5 is: 333.333333333332 Console.WriteLine("Approximate value using a relative error of 1e-5 is: " + integrate); 

Gauss-Legendre Rule

A fixed-order Gauss-Legendre integration routine is provided for fast integration of smooth functions with known polynomial order. The N-point Gauss-Legendre rule is exact for polynomials of order $$2N-1$$ or less. For example, these rules are useful when integrating basis functions to form mass matrices for the Galerkin method [GSL].

The basic idea of Gauss-Legendre integration is to approximate the integral of a function $$f(x)$$ using $$N$$ Weights $$w_i$$ and abscissas (or nodes) $$x_i$$.

$\int_a^b f(x) \, dx \approx \sum_{i = 0}^{N - 1} w_i f(x_i)$

This algorithm calculates the abscissas and weights for a given order and integration interval. For efficiency, pre-computed abscissas and weights for the orders $$N = 2 - 20, \, 32, \, 64, \, 96, 100, \, 128, \, 256, \, 512, \, 1024$$ are used. Otherwise, they are calculated on the fly using Newton's method. For more information on the algorithm see [Holoborodko, Pavel] .

Abscissas and Weights

We'll first use the abscissas and weights to approximate an integral using a 5-point Gauss-Legendre rule

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12:  // Create a 5-point Gauss-Legendre rule over the integration interval [0, 10] GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5); double sum = 0; // Will hold the approximate value of the integral for (int i = 0; i < rule.Order; i++) // rule.Order = 5 { // Access the ith abscissa and weight sum += rule.GetWeight(i) * rule.GetAbscissa(i) * rule.GetAbscissa(i); } // Approximate value is: 333.333333333333 Console.WriteLine("Approximate value is: " + sum); 

If you prefer direct access to the abscissas and weights, as opposed to using the methods

• double GetAbscissa(int i)
• double GetWeight(int i)

then use the properties Abscissas and Weights

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:  // Create a 5-point Gauss-Legendre rule over the integration interval [0, 10] GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5); double[] x = rule.Abscissas; // Creates a clone and returns array of abscissas double[] w = rule.Weights; // Creates a clone and returns array of weights double sum = 0; // Will hold the approximate value of the integral for (int i = 0; i < rule.Order; i++) // rule.Order = 5 { // Access the ith abscissa and weight sum += w[i] * x[i] * x[i]; } // Approximate value is: 333.333333333333 Console.WriteLine("Approximate value is: " + sum);; 

In addition to obtaining the abscissas and weights, the order and integration interval can be obtained

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:  // Create a 5-point Gauss-Legendre rule over the integration interval [0, 10] GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5); // The order of the rule is: 5 Console.WriteLine("The order of the rule is: " + rule.Order); // The lower integral bound is 0 Console.WriteLine("The lower integral bound is: " + rule.IntervalBegin); // The upper integral bound is 10 Console.WriteLine("The upper integral bound is: " + rule.IntervalEnd); 

Integrate Method

For convenience, we provide an overloaded static method double Integrate(...) which preforms 1D and 2D integration of a function. The first parameter to the method is a delegate of type Func<double, double> or Func<double, double, double> for 1D and 2D integration respectively. So for example

  1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:  // 1D integration using a 5-point Gauss-Legendre rule over the integration interval [0, 10] double integrate1D = GaussLegendreRule.Integrate(x => x * x, 0.0, 10.0, 5); // Approximate value of the 1D integral is: 333.333333333333 Console.WriteLine("Approximate value of the 1D integral is: " + integrate1D); // 2D integration using a 5-point Gauss-Legendre rule over the integration interval [0, 10] X [1, 2] double integrate2D = GaussLegendreRule.Integrate((x, y) => (x * x) * (y * y), 0.0, 10.0, 1.0, 2.0, 5); // Approximate value of the 2D integral is: 777.777777777778 Console.WriteLine("Approximate value of the 2D integral is: " + integrate2D); 

where we used $$\int_0^{10}\int_1^2 x^2 y^2 \,dydx = \frac{7000}{9} \approx 777.\overline{7}$$ for the 2D integral example.