Gene's Gaussian quadrature observations largely trivialize Simpson, but we still have the mystery of how 1 4 1 (3rd order) beats 7 32 12 32 7 (5th order) on 4/(1+x^2).
rcs>Maybe 1/(1+x) would be a better test function, since 4/(1+x^2) has
only even-exponent terms in the power series? Before even trying it, I'll say no, since we're not 0-centered. E.g., we could pretend we did 1/(1 + i x) and took the realpart afterwards. OK, now let's compute 1 / [ 1 I ----- dx = log(2) = K(0, 1, 2, 3, 1, 6, 3, 1, 1, 2, 1, 1, 1, 1, 3, ] x + 1 / 0
10, 1, 1, 1, 2, 1, 1, 1, 1, 3, 2, 3, 1, 13, 7, 4, 1, 1, 1, 7, 2, 4, ...)
by averaging 1/(1+x) over [0,1]. With nine samples, Thomas gives cf(''(makelist(1/(1+x),x,(0..8)/8).[1,4,2,4,2,4,2,4,1]/24))
[0, 1, 2, 3, 1, 6, 4, 1, 1, 2, 2, 5, 3, 1, 1, 1, 1, 1, 2]
and with *ten*, Homer gives cf(''(makelist(1/(1+x),x,(0..9)/9).[1,3,3,2,3,3,2,3,3,1]/24))
[0, 1, 2, 3, 1, 6, 5, 8, 2, 1, 2, 1, 8, 2, 1, 1, 5]
Worse, as predicted. On the other hand, the bizarreness with 7 32 12 32 7 giving -2 times the error with 1 4 1 may well depend on which function we're integrating, so 1/(1+x) will be good to test.
Sho nuff. I'd've bet my Heine this was just experimenter malfeasance, but, by Gauss, it's real. First, just to prove 7 32 12 32 7 really integrates quintics, (c41) sum(a[k]*x^k,k,0,6) 6 5 4 3 2 (d41) a x + a x + a x + a x + a x + a x + a 6 5 4 3 2 1 0 (c42) block([fancy_display:false],print(expand(integrate(%,x,0,1) <makelist(''%,x,(0..4)/4).[7,32,12,32,7]/90)),0)$ a a a a a a 55 a a a a a a 6 5 4 3 2 1 6 5 4 3 2 1 -- + -- + -- + -- + -- + -- + a < ----- + -- + -- + -- + -- + -- + a 7 6 5 4 3 2 0 384 6 5 4 3 2 0 (Note that the sides differ by a_6/384/7.) Here is 7 32... giving more than triple the absolute error for 32 samples of 4/(1+x^2): (c22) dfloat([makelist(4/(1+x^2),x,(0..32)/32) . ([32,12,32,7+7,32,12,32,7+7,32,12,32,7+7,32,12,32],append([7],%%,[14],%%,[7])/90/4/2),makelist(4/(1+x^2),x,(0..32)/32) . ([4,2,4,1+1,4,2,4,1+1,4,2,4,1+1,4,2,4],append([1],%%,[2],%%,[1])/12/4/2)]-%pi) (d22) [1.18244525282307d-10, - 3.69566599545123d-11] And for 64 samples, the error ratio gets slightly worse! (c24) dfloat([makelist(4/(1+x^2),x,(0..64)/64) . ([32,12,32,7+7,32,12,32,7+7,32,12,32,7+7,32,12,32],append(%%,[14],%%),append([7],%%,[14],%%,[7])/90/4/2/2),makelist(4/(1+x^2),x,(0..64)/64) . ([4,2,4,1+1,4,2,4,1+1,4,2,4,1+1,4,2,4],append(%%,[2],%%),append([1],%%,[2],%%,[1])/12/4/2/2)]-%pi) (d24) [1.84785520218611d-12, - 5.76871883595232d-13] (c25) d22[1]/d22[2] # %[1]/%[2] (d25) - 3.19954577680578d0 # - 3.20323325635104d0 Now try 1/(1+x). (c26) dfloat([makelist(1/(1+x),x,(0..64)/64).([32,12,32,7+7,32,12,32,7+7,32,12,32,7+7,32,12,32],append(%%,[14],%%),append([7],%%,[14],%%,[7])/90/4/2/2),makelist(1/(1+x),x,(0..64)/64).([4,2,4,1+1,4,2,4,1+1,4,2,4,1+1,4,2,4],append(%%,[2],%%),append([1],%%,[2],%%,[1])/12/4/2/2)]-log(2)) (d26) [3.61832785955584d-12, 1.86150950209907d-9] 7 32 ... did about as well, but Simpson (1 4 1) did far worse! So the question is not why 7 32 ... is lousy, but rather why is 1 4 1 so "lucky" on 4/(1+x^2). Obviously, we look at (c47) ('integrate(1/(1+x^3),x,0,1),%% = expand(apply_nouns(%%))) 1 / [ 1 log(2) sqrt(3) %pi (d47) I ------ dx = ------ + ----------- ] 3 3 9 / x + 1 0 (c48) dfloat([makelist(1/(1+x^3),x,(0..64)/64).([32,12,32,7+7,32,12,32,7+7,32,12,32,7+7,32,12,32],append(%%,[14],%%),append([7],%%,[14],%%,[7])/90/4/2/2),makelist(1/(1+x^3),x,(0..64)/64).([4,2,4,1+1,4,2,4,1+1,4,2,4,1+1,4,2,4],append(%%,[2],%%),append([1],%%,[2],%%,[1])/12/4/2/2)]-log(2)/3-sqrt(3)*%pi/9) (d48) [1.22610255282041d-12, 2.60732344048442d-9] Was Rich right about even:odd after all? 1 / [ 1 sqrt(2) (log(2 sqrt(2) + 3) + %pi) (d55) I ------ dx = ---------------------------------- ] 4 8 / x + 1 0 (c56) dfloat([makelist(1/(1+x^4),x,(0..64)/64) . ([32,12,32,7+7,32,12,32,7+7,32,12,32,7+7,32,12,32],append(%%,[14],%%),append([7],%%,[14],%%,[7])/90/4/2/2),makelist(1/(1+x^4),x,(0..64)/64) . ([4,2,4,1+1,4,2,4,1+1,4,2,4,1+1,4,2,4],append(%%,[2],%%),append([1],%%,[2],%%,[1])/12/4/2/2)]-rhs(%)) (d56) [6.0729199446996d-14, 1.98681771035325d-9] No, there's something magic about 1/(1+x^2) and 1 4 1 ! --rwg PS, to its credit, Mma 7.0 gets 1 1 1 1 / psi (--- + -) - psi (---) [ 1 0 2 n 2 0 2 n I ------ dx = ------------------------- ] n 2 n / x + 1 0 1 1 psi (-) - psi (---) - log(2) 0 n 0 2 n = ---------------------------- n by specializing a 2F1.