Comments and discussion on Valentin's 4th "Then and Now" - Area

[Albert Chan]

(02-12-2023 05:00 PM)Gjermund Skailand Wrote:  I found that double exponential integration is also well suited to evaluate the integrals.
...
Pushing the limits to e.g 1E-30 does results in a reported accuracy of 4E-32, but the results does not really change.
2.07669834394 7601111175686517522 91 (16.7 seconds, 941 function evaluations)

Thanks for DE integrate code!

A weakness with DE algorithm is we cannot have fuzzy integrand.
Above calculated area only had about half-precision accuracy.

Let's compare the same calculations, fuzzy vs "sharp" integrand:

>>> from mpmath import *
>>> mp.pretty = 1
>>> mp.dps = 34
>>> fnc = lambda x: sign(x)*abs(x)**(1/3)
>>> fnd = lambda y,s: fnc(y+s) - fnc(y-s)
>>> fne = lambda y: y*y*100/3007 + expm1(-sin(y))
>>> fnf = lambda y: fnd(y, sqrt(-log1p(fne(y))))
>>> main = [0, findroot(lambda y: fne(y)-expm1(-y*y),[0.5,1]), findroot(fne,[2.5,3])]
>>> speck = [findroot(fne, [-4.09,-4.08]), findroot(fne, [-4.05,-4.04])]

>>> main
[0, 0.8319711499790767979993060165734035, 2.82740261412956009186848044583741]
>>> speck
[-4.085146747638309747433733595149126, -4.049212264396950476554495084637107]

Note: mpmath quad use DE algorithm by default. Also, main had 2 parts, treated as 2 integrals.

>>> quad(fnf, main) + quad(fnf, speck)
2.076698343947601079521669308070185

The fuzziness comes from catastrophic cancellation of fnd
We redefine fnd, using identity: a^3 - b^3 = (a-b) * (a^2 + a*b + b^2)

>>> fnd = lambda y,s: 2*s/(fnc((y+s)**2) + fnc((y+s)*(y-s)) + fnc((y-s)**2))

>>> quad(fnf, main) + quad(fnf, speck)
2.076698343947601116938636757001192