Comparison of three quadrature rules

Sometimes, in science and engineering, we need to calculate some definite integrals of the type:

It is equivalent to evaluating the area enclosed between the function graph and x-axis, within the interval [a, b]. There are some integrals (like elliptical, Fresnel, etc.) that cannot be expressed in analytical form. In such cases, we have to use numerical methods, also called quadrature rules. They approximate the integral by a finite sum like follows:

where x_i are the abscissas of the nodes, where the function is calculated, and w_i are the respective weights. Most quadrature rules use the same formula, but with different values for abscissas and weights.

Probably, most engineers are already familiar with Newton’s or Gauss formulas for numerical integration. However, there are other contemporary methods that can be more efficient for this purpose.

In this post, we will compare three quadrature rules for numerical integration. They will be applied to the same function, and with the same number of nodes. The calculations will be performed with Calcpad, as follows:

Test function (half-circle):

Theoretical value of the integral:

Number of nodes: n = 5

Boole’s rule

This is one of the Newton-Cotes‘ family of quadrature formulas, with order of n = 4. Nodes are located at equal spacing h, and weights are calculated by approximating the function with Lagrange polynomials.

Spacing – h = 0.5

Abscissas	Ordinates
x1 = -1	y1 = f(-1) = 0
x2 = -0.5	y2 = f(-0.5) = 0.86603
x3 = 0	y3 = f(0) = 1
x4 = –x2 = 0.5	y4 = f(0.5) = 0.86603
x5 = –x1 = 1	y5 = f(1) = 0

Integral:

Error:

Scheme:

Gauss-Legendre’s formula

This is the most simple Gaussian quadrature formula. The function is approximated by Legendre polynomials in the interval of [-1; 1]. Abscissas are determined as the roots of the respective polynomials, so the nodes are not equally spaced.

Abscissas	Ordinates	Weights
x1 = -0.90618	y1 = f(-0.90618) = 0.42289	w1 = 0.23693
x2 = -0.53847	y2 = f(-0.53847) = 0.84265	w2 = 0.47863
x3 = 0	y3 = f(0) = 1	w3 = 0.56889
x4 = –x2 = 0.53847	y4 = f(0.53847) = 0.84265	w4 = w2 = 0.47863
x5 = –x1 = 0.90618	y5 = f(0.90618) = 0.42289	w5 = w1 = 0.23693

Abscissas and weights are calculated by the following equations:

Integral:

I_G = 2·w₁·y₁ + 2·w₂·y₂ + w₃·y₃ = 2·0.23693·0.42289 + 2·0.47863·0.84265 + 0.56889·1
I_G = 1.5759

Error:

Scheme:

Tanh-Sinh quadrature

This is one of the double exponential integration formulas, proposed by Takahasi and Mori (1974).

Boundary of the interval – t_a = 1

Step:

Parameter – t_k = –t_a + (k – 1) ·h

Formula for calculation of abscissas:

Formula for calculation of weights:

Abscissas	Ordinates	Weights
x1 = x(1) = -0.95137	y1 = f(-0.95137) = 0.30806	w1 = 0.11501
x2 = x(2) = -0.67427	y2 = f(-0.67427) = 0.73848	w2 = 0.48299
x3 = x(3) = 0	y3 = f(0) = 1	w3 = 0.7854
x4 = –x2 = 0.67427	y4 = f(0.67427) = 0.73848	w4 = w2 = 0.48299
x5 = –x1 = 0.95137	y5 = f(0.95137) = 0.30806	w5 = w1 = 0.11501

Integral:

Value:

I_DE = 2·w₁·y₁ + 2·w₂·y₂ + w₃·y₃ = 2·0.11501·0.30806 + 2·0.48299·0.73848 + 0.7854·1
I_DE = 1.5696

Error:

Scheme:

Summary of the results

The relative errors, from the above three quadrature rules are presented in the following chart:

We can see that the Gaussian quadrature rule is much more precise than the Newton-Cotes’ one of the same order, which is expected. However, Tanh-Sinh quadrature shows the best results from all three. And this is not by chance. According to Wikipedia, this is probably the most efficient numerical integration rule ever known.

That is why, it’s adaptive version is implemented in Calcpad as the main numerical integration method ($Integral command). Alternatively, you can use the adaptive Gauss-Lobatto quadrature ($Area command). It is much slower, but unlike the Tanh-Sinh one, it works for functions that are not smooth and continuous.

For the above example, both of them provides extremely high accuracy, that is comparable to the floating point precision of the double (float64) data type:

You can download the Calcpad source file from GitHub:

https://github.com/Proektsoftbg/Calcpad/blob/main/Examples/Mathematics/Numerical%20Integration.cpd