How does this numerical integral approximation stack up against the fastest converging Numerical Integration techniques

I came up with a numerical integration approximation sometime 2005 or so.

It employs Simpson's rule in computing the error incurred in computing the numerical integral of a function, say f(x) using a polynomial of best fit.

(DISCLAIMER: I do not know if anyone has come up with this formula before. I would like to know, please)

My technique computes the polynomial using a simple and I believe well known method...breaking up the curve into small units and finding the y values along those points.

Then solving for the coefficients of the resulting polynomial using a matrix method.

So if we wish to approximate y=sin(x), for example: say with a polynomial of degree n = 5,

between x = 1 and x = 3.

Then we would have:

dx = (3 - 1)/5; say dx = 0.4;


So x1 = 1, x2 = 1.4, x3 = 1.8, x4 = 2.2, x5 = 2.6, x6 = 3.

So the approximating polynomial gives:

 h(1.0) = sin(1)   = a0 + 1.0.a1 + 1.02.a2 + 1.03.a3 + 1.04.a4 + 1.05.a5
h(1.4) = sin(1.4) = a0 + 1.4.a1 + 1.42.a2 + 1.43.a3 + 1.44.a4 + 1.45.a5
h(1.8) = sin(1.8) = a0 + 1.8.a1 + 1.82.a2 + 1.83.a3 + 1.84.a4 + 1.85.a5
h(2.2) = sin(2.2) = a0 + 2.2.a1 + 2.22.a2 + 2.23.a3 + 2.24.a4 + 2.25.a5
h(2.6) = sin(2.6) = a0 + 2.6.a1 + 2.62.a2 + 2.63.a3 + 2.64.a4 + 2.65.a5
h(3.0) = sin(3.0) = a0 + 3.0.a1 + 3.02.a2 + 3.03.a3 + 3.04.a4 + 3.05.a5


This gives the polynomial:

h(x) = 0.05621527890686895 + 0.8070522394224585.x + 0.2659311952212848.x2 - 0.35406795117551315.x3 + 0.0696942723181288.x4 - 0.003354049885331352.x5

Then I compute the integral using polynomial integration of h(x), which is simple enough.

I then use Simpson's rule to estimate the error incurred in computing the integral, using the approximation:

Ae = S(f(x)) - S(h(x)).

Where S(f(x)) is Simpson's approximation for f(x) and S(h(x)) is Simpson's approximation for h(x).

By adding this error area to h(x) and summing all along the path, I get this formula

Please check out a detailed coverage of the technique here, if more detail is needed, as I dont want the question to spiral out of control.

I believe that this technique will always converge faster than Simpson's rule or the Polynomial approximation that I used. (even though it will always be more complex to program or use.) Is this true? If not, under what circumstances will it not be true?

Also, how does this technique match up against Gaussian quadrature and other macho-techniques accepted as best practices today?

Finally, I believe that any kind of polynomial approximation can be used to replace the type I used.

Is this true? what better polynomial approximations are there that can effectively replace the one I used in the formula?

Thanks a lot.

• How exactly does your method differ from a Newton-Cotes method? High order Newton-Cotes methods have long been recognized as being bad for a variety of reasons. The main one is simply that the evenly-spaced polynomial interpolant of a continuous function $f$ doesn't even converge to $f$ in general. Nor does it even if you assume $f$ is smooth. – Ian Aug 4 '16 at 15:13
• As for your "correction", this is very similar in spirit to Richardson extrapolation, which in the case of integration rules is called Romberg integration. – Ian Aug 4 '16 at 15:19
• I believe any polynomial approximation can be used. But for this particular one that I am using. I want to know if the error estimate that I used is efficient enough.I used Simpson's rule to estimate the error in computation using the polynomial method. I have tested it on many integrals and I have seen it increase the accuracy a lot. Between polynomial degrees 5 and 10, for many examples the integral, converges fast it seems. For y=sin(x) between x =1 and x=3, using a polynomial of order 5, this technique gives:1.5302948297254457. Real value is -1.5302948024685852 of -1.5302948024685852 – gbenroscience Aug 4 '16 at 15:25
• Okay..."Similar in spirit"? Are they absolutely the same? I mean the Romberg integration part. as in will they give exactly same results if applied to the same problems? – gbenroscience Aug 4 '16 at 15:29
• No, they're not quite the same, in fact your technique is not even guaranteed to improve the order of the original method. But the idea is the same: do some further calculation to get the "leading order term" of the error and then add that back in to hopefully reduce the error. – Ian Aug 4 '16 at 15:37

It looks like you are integrating a Lagrange polynomial on some evenly-spaced samples $(x_i,f(x_i))$ of the function. With $n$ points, your technique integrates polynomials of up to $(n-1)$th order exactly.
This is not optimal: by cleverly (and non-uniformly) choosing the locations of the sampling points $x_i$, an $(2n+1)$th order accurate integration can be done using only $n$ points. See Gaussian quadrature.