Talk:Numerical integration/Gauss-Legendre Quadrature

How can we tell if an implementation of this task is correct? (Ok, if we copy the lisp code and get the same answers that would probably be correct, but that does not address all the abstractions raised in the task description.) --Rdm 13:12, 30 May 2011 (UTC)

Analytically, this is true:

$\int _{-3}^{3}e^{x}\,dx=e^{3}-e^{-3}\approx 20.035749854819805$

That last value is what I get when I compute it directly (using my platform's libc's implementation of the exponential function for double-precision IEEE arithmetic). A correct implementation should tend toward that value in this case as the order of the Legendre polynomial being used is increased, up to a limit where it no longer helps. I observe that with the Tcl code, where I've tested up to 13 it indeed improves, but beyond that it ceases to help, becoming slow but getting no more accurate. I guess that this is running into the limits of the accuracy of the method for this type of equation (and note that it is delivering about 9 decimal digits of accuracy at that point; pretty good…) –Donal Fellows 09:31, 2 June 2011 (UTC)

How did you decide on the number of Newton-Raphson iterations to use when computing the node values? Wouldn't that also affect the accuracy? (I passed the required precision as a parameter and iterated until two consecutive estimates agreed up to that precision.) Exact floating point representations of the coefficients are also possible with sufficient precision. Increasing the precision always seems to give more accurate results, but having too many nodes can strangely make it worse. --Sluggo

This is a two-fold problem. In principal, the higher order you go, the better you can decompose your function into a sum of polymials, so with infinite order you can get exact result -- on paper. In a computer, every time you do math with a floating point number, you get an error due to precision (up to $\sim 10^{-16}$ of the value with IEEE 64 bit in general), and later you use the result to do more math, the error propagates and gets larger and larger. Higher orders require a lot more arithmetic operations (about $O(n^{2})$ I think), so at some point, the precision error dominates and more terms only make results worse. --Ledrug 07:39, 28 June 2011 (UTC)

Though, of course, you could use something other than ieee floating point to represent your numbers -- perhaps arbitrary precision rationals? --Rdm 11:06, 28 June 2011 (UTC)

Yeah but 1) Legendre polynomial roots are not always rational; 2) it's sloooow, if you go down that route, you might as well use some other way to do integrals. --Ledrug 14:32, 28 June 2011 (UTC)