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:


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



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


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


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:


 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.


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:



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; } 

Source: https://habr.com/ru/post/479202/

All Articles