Stand-alone integration rules
Integration rules are used to approximate the integral of a function over a given interval, most often using a weighted sum of function values at certain specific points. The points are called nodes, and the weights are called weights.
The order of a rule is the degree of the polynomial that the rule can integrate exactly. For example, a rule of order 3 can integrate exactly any polynomial of degree 3 or less.
A rule is nested in another rule if the nodes of the smaller rule are a subset of the nodes of the larger rule. Nested rules are useful because they provide two approximations of the integral with different accuracies, which provides a way to estimate the error of the approximation.
In the previous section we discussed some very basic rules that use only one or two function values. In order to compute a result with sufficient accuracy, the integration interval is divided in many subintervals to which the rule is applied.
In this section we will discuss more advanced rules that compute the integral over the entire interval in one step. Integration rules can be created by calling the static methods of the IntegrationRules class.
We will start by discussing the rules from the previous section in terms of a more general class to which they belong: Newton-Cotes rules
Newton-Cotes rules
Newton-Cotes rules use the values of the function at equally spaced points in the interval to approximate the integral. The simplest rules use only one function value. The left and right point rules use either of the end points, while the midpoint rule uses the point at the middle of the interval.
More complex rules can be built by combining the simple rules. The trapezoidal rule uses the average of the left and right point rules, while Simpson's rule uses the average of the left and right point rules and the midpoint rule.
Using only a few points to approximate an integral is not very accurate. It is possible to compute the integral with higher accuracy by using more points. Unfortunately, the error of the Newton-Cotes rules generally does not decrease when we increase with the number of points.
The root cause is that fixing the points at which the function is evaluated does not allow for the best possible approximation.
Clenshaw-Curtis Rules
Clenshaw-Curtis quadrature approximates a function by expanding it in terms of Chebyshev polynomials. The integral is then approximated by the weighted sum of the integrals of the Chebyshev polynomials. The equations can be re-written in terms of the function values at the nodes.
There are two variants of Clenshaw-Curtis quadrature: the open and closed rules. The open rule uses the function values at the zeros of the Chebyshev polynomials, while the closed rule uses the function values at the extrema, which includes the end points.
The advantage of the closed rule is that it can be nested, and so they can provide an estimate of the error.
To create a Clenshaw-Curtis rule, call the ClenshawCurtis method. The argument is the total number of nodes to use, including the end points. This creates a closed rule that can be nested.
To create an open Clenshaw-Curtis rule, call the ClenshawCurtisOpen method. The argument is again the total number of nodes to use.
Gauss Quadrature
In Gauss quadrature, the nodes and weights are chosen to be optimal in the sense that the rule integrates exactly all polynomials of the highest possible degree, which is \(2n - 1\) where \(n\) is the number of points. This basic rule is called the Gauss-Legendre rule.
Gauss rules can be created for specific weight functions. Instead of integrating the function \(f(x)\), these Gauss rules integrate \(w(x) f(x)\), where \(w(x)\) is the weight function. The nodes and weights are still chosen to be optimal, but the rules are exact for the weight function times a polynomial of degree \(2n - 1\). This is especially useful when the weight function has singularities or is otherwise difficult to integrate.
The construction of Gauss integration rules uses the polynomials that are orthogonal with respect to the weight function on the interval of integration. The rules get their name from the class of orthogonal polynomials.
The standard Gauss rules are for the weight function \(w(x) = 1\), and the corresponding polynomials are the Legendre polynomials. And so the basic Gauss rule is the Gauss-Legendre rule.
All Gauss rules can be created by calling the static methods of the IntegrationRules class. The methods are listed below:
Rule | Interval | Weight function |
---|---|---|
$$[-1,1]$$ | $$1$$ | |
$$[-1,1]$$ | $$\frac{1}{\sqrt{1-x^2}}$$ | |
$$[-1,1]$$ | $$\sqrt{1-x^2}$$ | |
$$[-1,1]$$ | $$(1-x^2)^\alpha$$ | |
$$[-1,1]$$ | $$(1-x)^\alpha(1+x)^\beta$$ | |
$$[0,+\infty]$$ | $$x^\alpha e^{-x}$$ | |
$$[-\infty,+\infty]$$ | $$e^{-x^2}$$ |
Gauss rules are open: they do not include the end points of the interval. A variant of the Gauss rules that includes the end points is called Gauss-Lobatto rules. Fixing two of the nodes at the end points of the interval means that the order is also reduced by two: an \(n\) point rule has order \(2n-3\).
Gauss-Lobatto rules can be created by calling the GaussLobatto method for the standard rule, or GaussJacobiLobatto for the Jacobi rule.
Computing Integrals
To compute an integral using an integration rule, call the rule's Evaluate method. The first argument is the function to integrate. The second and third arguments are optional. They are the lower and upper bounds of the integration interval.
The return value is of type IntegrationRuleResult. The Value property contains the result of the integration, while the Error property contains an estimate for the error, if it is available.
As an example, let's integrate the function \(f(x) = e^x / (x + 1)^{1/2}\) over the interval \([-1,1]\). We'll use the following setup:
double exactValue = 2.4602620138961554780;
Console.WriteLine($"Integrate Exp(x) / Sqrt(1 + x) over [-1,+1]");
Console.WriteLine($"Exact result: {exactValue}");
We can use the Gauss-Legendre rule with 20 points:
Func<double, double> f = x => Math.Exp(x) / Math.Sqrt(1 + x);
var gl = IntegrationRules.GaussLegendre(20);
var result = gl.Evaluate(f);
var error = result.Value - exactValue;
Console.WriteLine($"20 pt. Gauss-Legendre: {result.Value} ({error})");
This returns a result of 2.438... Unfortunately, even though we used a decent number of points, the error is quite large and far from acceptable.
A closer look at the integrand reveals that we can improve the accuracy by using a Gauss rule with a suitable weight function. If we choose a Gauss-Jacobi rule with \(\alpha = 0\) and \(\beta = -1/2\), then the integrand just becomes \(e^x\). This is a much simpler function and should be easier to integrate. Let's try with just 6 points:
var gj = IntegrationRules.GaussJacobi(6, 0.0, -0.5);
result = gj.Evaluate(Math.Exp);
error = result.Value - exactValue;
Console.WriteLine($"6 pt. Gauss-Jacobi: {result.Value} ({error})");
The error is only about 2e-12, an excellent result given that we only used 6 function evaluations.
When a rule is evaluated over a non-standard interval, then the nodes are transformed and the weights are scaled accordingly.
It is important to note that any weight function is transformed as well. If the width of the new interval does not equal the width of the standard interval, then the weight function will have a different shape. This is not an issue for integration rules like Gauss-Legendre and Clenshaw-Curtis that don't have a weight function (or whose the weight function is equal to 1).
References
C. W. Clenshaw and A. R. Curtis, A method for numerical integration on an automatic computer, Numerische Mathematik 2, 197, (1960).
C. F. Gauss, Methodus nova integralium valores per approximationem inveniendi., Comm. Soc. Sci. Göttingen Math. Vol. 3. S. 29–76, 1814.
C. G. J. Jacobi, Ueber Gauß' neue Methode, die Werthe der Integrale näherungsweise zu finden, Journal für die Reine und Angewandte Mathematik. 1. S. 301–308, 1826.
R. Piessens, E. de Doncker, C. Uberhuber and D. Kahaner, QUADPACK, A Subroutine Package for Automatic Integration , Springer-Verlag, New York, 1983.