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