Numerical Methods in 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  .

Approximation of integral of f(x) between 0 and pi/2 by trapezoids

Approximation of integral of f(x) between 0 and pi/2 by trapezoids

 

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.

Approximation of integral of f(x) between 0 and pi/2 by parabola

Approximation of integral of f(x) between 0 and pi/2 by a parabola

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  .

Romberg parable with nodes at T(1/2) and T(1/4)

Romberg parable with nodes at T(1/2) and T(1/4)

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}

Screen capture of program execution

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?

What are your thoughts?

Feel free to comment and share this post.

2 thoughts on “Numerical Methods in C++ Part 1: Newton-Cotes Integration

  1. 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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

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