Integral Approximation - Simpson’s Rule¶
Definition¶
Suppose \(f(x)\) is defined on the interval \([a, b]\). Then Simpson’s rule on the entire interval approximates the definite integral of \(f(x)\) on the interval by the formula
The idea is that if \(f(x) = 1, x, x^2\), this formula is an exact equality. So Simpson’s rule gives the correct integral of any quadratic function. In general, Simpson’s rule approximates \(f(x)\) by a parabola through the points on the graph of \(f(x)\) with \(x\)-coordinates \(a, \frac{a+b}{2}, b\).
Simpson’s rule is usually applied by breaking the interval into \(N\) equal-sized subintervals, where \(N\) is an even number, and approximating the integral over each pair of adjacent subintervals using the above estimate.
That is, let \(x_0=a, x_1=a+\frac{b-a}{N}, x_2=a+2\frac{b-x}{N}, \dots, x_{N-1}=a+(N-1)\frac{b-a}{N}, x_N=b\). Then
and so on. Adding these up gives
Usage¶
Imagine that we want to integrate the following expression:
Then the code will look like this:
// example_integral_simpson.cpp
#include <iostream>
#include "../src/numerary.hpp"
using namespace std;
using namespace numerary;
/* Function to be integrated */
double f(double x) {
return 5*pow(x, 3) + 2*cos(x);
}
/* The main function */
int main() {
const double from = 0; // Lower bound of integral
const double to = 1; // Upper bound of integral
const string method = "simpson"; // Numerical method we will use for integration ("simpson" by default)
const double eps = 1.e-9; // eps value for integration (1.e-9 by default)
double *I = Numerary::integrate(f, from, to, method, eps);
cout << "ans = " << I[0] << endl; // Value of calculated integral
cout << "err = " << I[1] << endl; // Error of calculated integral value
return 0;
}