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.