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 |

x_{1} = -1 | y_{1} = f(-1) = 0 |

x_{2} = -0.5 | y_{2} = f(-0.5) = 0.86603 |

x_{3} = 0 | y_{3} = f(0) = 1 |

x_{4} = –x_{2} = 0.5 | y_{4} = f(0.5) = 0.86603 |

x_{5} = –x_{1} = 1 | y_{5} = 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 |

x_{1} = -0.90618 | y_{1} = f(-0.90618) = 0.42289 | w_{1} = 0.23693 |

x_{2} = -0.53847 | y_{2} = f(-0.53847) = 0.84265 | w_{2} = 0.47863 |

x_{3} = 0 | y_{3} = f(0) = 1 | w_{3} = 0.56889 |

x_{4} = –x_{2} = 0.53847 | y_{4} = f(0.53847) = 0.84265 | w_{4} = w_{2} = 0.47863 |

x_{5} = –x_{1} = 0.90618 | y_{5} = f(0.90618) = 0.42289 | w_{5} = w_{1} = 0.23693 |

Abscissas and weights are calculated by the following equations:

**Integral**:

`I`_{G} = 2·`w`_{1}·`y`_{1} + 2·`w`_{2}·`y`_{2} + `w`_{3}·`y`_{3} = 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 |

x_{1} = x(1) = -0.95137 | y_{1} = f(-0.95137) = 0.30806 | w_{1} = 0.11501 |

x_{2} = x(2) = -0.67427 | y_{2} = f(-0.67427) = 0.73848 | w_{2} = 0.48299 |

x_{3} = x(3) = 0 | y_{3} = f(0) = 1 | w_{3} = 0.7854 |

x_{4} = –x_{2} = 0.67427 | y_{4} = f(0.67427) = 0.73848 | w_{4} = w_{2} = 0.48299 |

x_{5} = –x_{1} = 0.95137 | y_{5} = f(0.95137) = 0.30806 | w_{5} = w_{1} = 0.11501 |

**Integral**:

**Value**:

`I`_{DE} = 2·`w`_{1}·`y`_{1} + 2·`w`_{2}·`y`_{2} + `w`_{3}·`y`_{3} = 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