Adaptive Simpson and Romberg integration methods

[Albert Chan]

(03-25-2021 01:58 AM)robve Wrote:  Weird that HP-71B Math Pac uses trapezoidal quadratures when trying to exclude endpoints.

That may be because trapezoids gives better estimate than mid-point rules ?

Lets try \(\displaystyle \int_1^e {dx \over x} = \ln(e) - \ln(1) = 1 \)

To simplify, we skip u-transformation, and supplied first trapezoid (t1) instead.
(we got the same conclusion either way, it is just easier to visualize 1/x)

lua> I = require'integ'       -- code posted here
lua> e = exp(1)
lua> f = function(x) return 1/x end
lua> t1 = (f(1)+f(e))/2 * (e-1)
lua> g = I.tm(f,1,e,t1)      -- generate trapezoids and mid-points area (both using the *same* points)

lua> for i=1,6 do t,m = g(); print(2^i, (1-m)/(1-t)) end

2   -1.5239160915264984
4   -1.8129951785930594
8   -1.9427734883703107
16  -1.984685010021192
32  -1.996097578762125
64  -1.99901957558112

Trapezoids seems to be twice as good as mid-points (in the opposite direction).
When we extrapolate with Romberg, trapezoids are *much* better

lua> I.integ(I.simpson(f,1,e,t1))

n    T = extrapolated-t   M = extrapolated-m   error(M)/error(T)
4    1.0003409768323182   0.9937356332778096   --> -18
8    1.000007899493941    0.9996852308223878   --> -40
16   1.0000000801614912   0.9999923219175764   --> -96
32   1.0000000003184022   0.9999999206312571   --> -249
64   1.00000000000046     0.9999999996826735   --> -690
...

n is sub-intervals for trapezoids, which is halved, for M
Ignoring end-points, both T and M required n-1 function evaluations.

With same sub-intervals, both scheme have comparable errors. Example:

error(T32) = 1 - 1.0000000003184022 ≈ −3.18e-10
error(M64) = 1 - 0.9999999996826735 ≈ +3.17e-10

But, getting M64 required doubled the time than T32.