Numerical Methods with C++ Part 1: Newton-Cotes Integration

Welcome to a new post on thoughts-on-cpp.com. This time I would like to start a new series, Numerical Methods. But don’t be disappointed if you’re expecting a new post on the n-body-problem, I’m still planning to continue the “My God, It’s Full of Stars” series. With this new series, I would like to talk about numerical methods which are important in my daily business, starting from Trapezoidal and Simpson integration methods of type Newton-Cotes.

The Newton-Cotes formula is a quadratic numerical approximation for integral calculations. The idea is to interpolate the function, which shall be integrated, by a polynomial with equidistant nodes. The polynomial can be of, for example, a form of aligned trapezoids or aligned parables. Common Newton-Cotes integration polynomial rules are

With this post, we will have a closer look at the first two integration rules, Trapezoidal and Simpson.

Trapezoidal Integral Approximation

The integral approximation by trapezoids is very simple to explain. All we need to do is to subdivide the example function

$I=\int_0^{\pi/2} \frac{5}{e^\pi-2}\exp(2x)\cos(x)dx=1.0$

we want to integrate into equidistant areas whose exact integrations we can sum up. Clearly, the accuracy of the approximation is depending on the number of subdivisions N. The width of each area is therefore defined by

$\Delta x = \frac{b-a}{N}$

with $a=0$ and $b=\pi/2$.

The area of each trapezoid can be calculated by

$A=\frac{a+b}{2}\cdot h$

and therefore the approximated integral $\widetilde{ I }$ of our function f(x) can be defined as

$\widetilde{ I } = \sum_{i=1}^{N+1} \frac{x_i-x_{i+1}}{2}(f(x_i)+f(x_{i+1}))$

The implementation of the Trapezoidal integration is taking 4 parameters, the range [a,b] of integration of the function f(x), the number of subdivisions n, and the function f(x).

Simpson Integral Approximation

The Simpson rule is approximating the integral of the function f(x) by the exact integration of a parable p(x) with nodes at ab, and $m=\frac{a+b}{2}$. In order to increase the approximation accuracy, the function can be subdivided by N, similar to the Trapezoidal integral approximation.

The exact integration can be done by summing up the area of 6 rectangles with equidistant width. The height of the first rectangle is defined by $f(x_a)$, the height of the next 4 rectangles is defined by $f(\frac{a+b}{2})$, and the height of the last rectangle is defined by $f(x_b)$. As a result, the formula of the approximated integral according to Simpson rule can be defined as

$\widetilde{ I } = \sum_{i=1}^{N+1} \frac{x_i-x_{i+1}}{6}(f(x_i)+4f(\frac{a+b}{2})+f(x_{i+1}))$

The implementation of the Simpson integration is, similar to the Trapezoidal based solution, taking 4 parameters, the range [a,b] of integration of the function f(x), the number of subdivisions n, and the function f(x)

Romberg Integral Approximation

If we integrate the function f(x) with the Trapezoidal approach but bisecting the step size based on the step size of the previous step, we get the approximation sequence

$T(1)=0.18755, \> T(\frac{1}{2})=0.724727, \> T(\frac{1}{4})=0.925565$

The idea of the Romberg integration is now to introduce a y-axis symmetric parable which is crossing the points $T(1), \> T(\frac{1}{2})$ and extrapolate to $x \rightarrow 0$.

Therefore every term of the first column R(n,0) of the Romberg integration is equivalent to the Trapezoidal integration, were every solution of the second column R(n,1) is equivalent to the Simpson rule, and every solution of the third column R(n,2) is equivalent to Boole’s rule. As a result, the Formulas for the Romberg integration are

$R(0,0)=h_{1}(f(a)+f(b))$

$R(n, 0)=\frac{1}{2} R(n-1,0)+h_{n} \sum_{k=1}^{2^{\pi-1}} f\left(a+(2 k-1) h_{n}\right)$

$R(n, m)=R(n, m-1)+\frac{R(n, m-1)-R(n-1, m-1)}{4^{m}-1}$

The repository numericalIntegration can be found at GitHub. It will also contain the other numerical integration methods, as the name suggests, later on.

Did you like the post?

Processing…
Success! You're on the list.

3 thoughts on “Numerical Methods with C++ Part 1: Newton-Cotes Integration”

1. stilgar says:

I know that multiplication is commutative and it technically is correct, but I think your trapezoidal equation is wrong. Well, maybe not wrong, but just hard to read. You wrote it first as (a+b)/2 * h and then as (x1-x0)/2 * (f(x1)+f(x0)). Someone might wonder why h = f(x1)+f(x0), when h is in fact x1 – x0.

It would change the second equation so it matches the first one. Just a minor nitpick 🙂

Like

1. Thanks for your suggestions, I will have a look

Like

This site uses Akismet to reduce spam. Learn how your comment data is processed.