Misc thoughts ... Maybe 1/(1+x) would be a better test function, since 4/(1+x^2) has only even-exponent terms in the power series? Try just using one cycle for the whole integral: 11, 141, 1331, etc. Is the convergence interesting? Does it? What happens if you stick with Simpson's 141, but move around the cut points (the abscissae)? Rich ------------ Quoting rwg@sdf.lonestar.org:
Merriam-Webster: simpson's rule Function: noun Usage: usually capitalized S Etymology: after Thomas Simpson died 1761 English mathematician
: a method used especially by naval architects for computing the approximate area bounded by a curve by adding the areas of a series of figures formed from an odd number of equally spaced ordinates to the curve and parabolas drawn through the points where these ordinates cut the curve
(I always wondered what they did in Course XIII.) This is a complicated description of the 1,4,2,4,...,2,4,1 method. However, I have to think that it was not Thomas, but rather Homer, who authored the "3/8 rule", 1,3,3,2,3,3,2,...,3,3,1. E.g., let's compute pi = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, ...
by averaging 4/(1+x^2) over [0,1]. With nine samples, Thomas gives cf(''(makelist(4/(1+x^2),x,(0..8)/8).[1,4,2,4,2,4,2,4,1]/24))
[3, 7, 15, 1, 186, 2, 6, 1, 1, 4, 21, 1, 73, 11]
and with *ten*, Homer gives cf(''(makelist(4/(1+x^2),x,(0..9)/9).[1,3,3,2,3,3,2,3,3,1]/24))
[3, 7, 15, 1, 127, 4, 2, 18, 6, 1, 78, 1, 4],
and for comparably many samples, never beats Thomas. Is there any reason to use Homer's rule besides being stuck with 6 n + 4 samples? Interestingly, the period 6 pattern 41, 216, 27, 272, 27, 216, 41+41, 216, ..., 41, which exactly integrates heptics,
expand([41, 216, 27, 272, 27, 216, 41]/840 .makelist(a*x^7+b*x^6,x,(0..6)/6))
a b - + - 8 7
overtakes Thomas on the pi example at nineteen samples: cf(''([1,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,1]/6/9 .makelist(4/(1+x^2),x,(0..18)/18))) [3, 7, 15, 1, 291, 2, 1, 4, 5, 1, 29, 2, 8, 1, 1, 1, 1, 2, 1, 1, 1, 8, 2, 672, 1, 2, 2, 2, 4, 89, 1, 1, 1, 1, 1, 2, 4, 22, 1, 11, 2, 1, 1, 1, 1, 6]
cf(''([41,216,27,272,27,216,41+41,216,27,272,27,216,41+41,216,27,272,27,216,41]/840/3 .makelist(4/(1+x^2),x,(0..18)/18))) [3, 7, 15, 1, 292, 1, 1, 4, 1, 37, 5, 3, 112, 1, 10, 2, 7, 2, 23, 1, 1, 3, 1, 1, 2, 1, 2, 2, 2, 2, 10, 1, 3, 2, 11, 1, 6, 3, 5, 4, 1, 2, 6, 2, 3, 1, 2]
But something really weird is going on. The period 4 weighting pattern 7, 32, 12, 32, 7+7, 32, 12,..., 32, 7, which exactly integrates quintics, has "exactly" (modulo double precision) -2 times Thomas's error with 201 samples! Perhaps this will make sense in the morning. Meanwhile, two applications of a 0th order Euler-Maclaurin expansion give an interesting formula for Thomas's error:
'integrate(g(k),k,0,1)-('sum(4*g(k/n+1/(2*n))+2*g(k/n),k,0,n-1)+g(1)-g(0))/(6*n) = ('sum('integrate((1/6-x)*'diff(g(x/n+i/n)+g((1-x)/n+i/n),x,1),x,0,1/2),i,0,n-1))/n
(d296) n - 1 ==== \ k 1 k > (4 g(- + ---) + 2 g(-)) + g(1) - g(0) 1 / n 2 n n / ==== [ k = 0 I g(k) dk - ------------------------------------------- = ] 6 n / 0
1 - n - 1 2 ==== / \ [ 1 d x i 1 - x i > I (- - x) (-- (g(- + -) + g(----- + -))) dx / ] 6 dx n n n n ==== / i = 0 0 --------------------------------------------------, n
(directly verifiable by integration by parts). The individual integrals on the right are small. Can someone explain why? Testing with g(x):=4/(1+x^2), n = 1..4,
(c346) makelist(subst(lambda([[l]],funmake_no_simp("+",l)),"[",lhs(e) = dfloat(rhs(e))),e,block([listarith:False],makelist(expand(eval(''(subst(makelist,nounify(sum),d296)))),n,1,4)))
47 (d346) [%pi - -- = + 0.00825932025646d0, 15
8011 %pi - ---- = 3.35550919902339d-4 - 3.11524781089349d-4, 2550
829597 %pi - ------ = 3.48105078317297d-5 + 5.09316464365226d-5 264069
152916620159 - 8.48695005182994d-5, %pi - ------------ = 48674874300
7.09747270259224d-6 + 1.325010624272d-5 + 4.8192196463089d-6
- 2.5015667505753d-5]
(c347) resimplify(dfloat(%))
(d347) [0.00825932025646d0 = 0.00825932025646d0, 2.40261388126939d-5 = 2.402613881299d-5, 8.72653749706131d-7 = 8.72653749952849d-7, 1.51131086312262d-7 = 1.51131085868171d-7]
(The inexactitudes are floating point artifacts.) --rwg DAMN YOUR ROUND YAM PEDIATRIC PATRICIDE SCHEMATIC CATECHISM ALLIGATOR WEED WATER GLADIOLE TRANSPORTEE PATERNOSTER PENETRATORS RHINOPLASTIES RELATIONSHIPS SUPERIMPOSED PSEUDOPRIME COMMERCIALIST MICROCLIMATES GRAPHITE EARTH PIG (Hmm, the spellchecker mailsquirrel is having a nosebleed.)
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun