The Simpson Method

The Simpson method uses quadratic approximations of the function on each subinterval (three function values, at , and are used to determine uniquely the quadratic polynomial). We add up all the integrals of the quadratic polynomials on all the subintervals, to obtain an approximation of the function integral over . This approximation is (the Simpson formula)

Note that the weights are

An efficient implementation of the algorithm that computes for different values of might use 3 variables:

The integral is

To implement the adaptive Simpson method, we need to compute efficiently, re-using the function values that were calculated for . The key is to observe that, when doubling the number of subintervals, the old node points are a subset of the new node points become are a subset of the node points .