# C ++ and Numerical Methods: Approximate Newton-Cotes Integration

Newton-Cotes methods are a combination of approximate integration techniques based on:

• splitting the integration interval into equal intervals;
• approximations of the integrand on selected intervals by polynomials;
• finding the total area of ​​the obtained curved trapezoid.

This article will cover several Newton-Cotes methods:

• trapezoid method;
• Simpson's method;
• Romberg method.

## Trapezoid method

The trapezoid method is the simplest of those considered. As an example, take the following integral:

$image$

The accuracy of the approximation depends on the number N of segments into which the integration interval is divided. Thus, the gap length:

$image$

The area of ​​the trapezoid can be calculated by the formula:

$image$

Summarizing all of the above, the approximate value of the integral is calculated by the formula:

$image$

The function that calculates the integral by the trapezoid method should take 4 parameters:

• boundaries of the integration segment;
• integrand function;
• the number N of intervals of the partition.

``double trapezoidalIntegral(double a, double b, int n, const std::function<double (double)> &f) { const double width = (ba)/n; double trapezoidal_integral = 0; for(int step = 0; step < n; step++) { const double x1 = a + step*width; const double x2 = a + (step+1)*width; trapezoidal_integral += 0.5*(x2-x1)*(f(x1) + f(x2)); } return trapezoidal_integral; }` `

## Simpson Method

The Simpson method consists in integrating the second-degree interpolation polynomial of the function f (x) with interpolation nodes a, b and m = (a + b) / 2 - parabolas p (x). To increase the accuracy, it makes sense to divide the integration segment into N equal intervals ( by analogy with the trapezoid method), on each of which apply the Simpson method.

The area of ​​the parabola can be found by summing the areas of 6 rectangles of equal width. The height of the first of them should be equal to f (a), from the third to the fifth - f (m), the sixth - f (m). Thus, the approximation by the Simpson method is found by the formula:

$image$

` `double simpsonIntegral(double a, double b, int n, const std::function<double (double)> &f) { const double width = (ba)/n; double simpson_integral = 0; for(int step = 0; step < n; step++) { const double x1 = a + step*width; const double x2 = a + (step+1)*width; simpson_integral += (x2-x1)/6.0*(f(x1) + 4.0*f(0.5*(x1+x2)) + f(x2)); } return simpson_integral; }` `

## Romberg Method

Let T (x) be the integral approximation obtained by the trapezoid method with step x. We get 3 such approximations, reducing the step size by 2 times with each calculation.

$image$

We now construct a parabola symmetrical with respect to the y axis, passing through the points T (1) and T (1/2) in order to extrapolate the obtained values ​​for x tending to 0.

Therefore, each member of the first column R (n, 0) of the Romberg approximations is equivalent to the solutions obtained by the trapezoidal method, and each solution of the second column of R (n, 1) is equivalent to the Simpson method. Thus, the formulas for approximate integration by the Romberg method:

$image$

$image$
$image$

C ++ implementation:

` `std::vector<std::vector<double>> rombergIntegral(double a, double b, size_t n, const std::function<double (double)> &f) { std::vector<std::vector<double>> romberg_integral(n, std::vector<double>(n)); romberg_integral.front().front() = trapezoidalIntegral(a, b, 1, f); double h = ba; for(size_t step = 1; step < n; step++) { h *= 0.5; double trapezoidal_integration = 0; size_t stepEnd = pow(2, step - 1); for(size_t tzStep = 1; tzStep <= stepEnd; tzStep++) { const double deltaX = (2*tzStep - 1)*h; trapezoidal_integration += f(a + deltaX); } romberg_integral[step].front() = 0.5*romberg_integral[step - 1].front() + trapezoidal_integration*h; for(size_t rbStep = 1; rbStep <= step; rbStep++) { const double k = pow(4, rbStep); romberg_integral[step][rbStep] = (k*romberg_integral[step][rbStep-1] - romberg_integral[step-1][rbStep-1])/(k-1); } } return romberg_integral; }` `

All Articles