(PC-12xx~14xx) qthsh Tanh-Sinh quadrature

[Albert Chan]

(04-02-2021 02:46 PM)robve Wrote:  Functions that are typically integrated over non-[0,1] finite domains are normalized to the [0,1] domain in experiments.
It is simpler that way to conduct the tests and for comparisons.

Here is a link to convert integral limits to [0, 1]: http://fmnt.info/blog/20180818_infinite-integrals.html

Example, my previous post e-transform, without changing integral limits.

Code:

def e_transform2(f, exp=exp):  # interval [0,1] -> [0,1]
    def ef(z):
        try: return f(z)*z
        except: return 0
    return lambda t: ef(exp(1-1/t))/(t*t)

>>> f = lambda x: x**-0.99
>>> ef2 = e_transform2(f)

>>> quad(f, [0,1])
mpf('35.65470993432163')
>>> quad(ef2, [0,1])
mpf('99.999999999999915')

Comment: although convenient to have unit interval, integral may be harder to compute.
Conversion of [0,∞] -> [0,1], compress the curve, possibly making it spiky.

Cas> f(x) := x^-0.99
Cas> ef1(t) := f(e^-t)*e^-t;     // [0,1] -> [0,inf]
Cas> ef2(t) := ef1(1/t-1)/(t*t); // [0,inf] -> [0,1]

ef1(t) = exp(-0.01*t), no spike, thus easy to integrate.
ef2(t) has a spike at t = 0.005

Cas> ef2(0.007) → 4939.99121224
Cas> ef2(0.005) → 5467.81701782
Cas> ef2(0.003) → 4003.61366011
Cas> ef2(0.001) → undef
Cas> limit(e2(t), t=exact(0.001)) → 1000000*exp(98901/100)/exp(999)

Numerically, e2(0.001) will be hit with overflow, ∞/∞ → undef (*)

(*) mpmath does not overflow, nor underflow

>>> 1000000*exp(98901/100)/exp(999)
mpf('45.856206642206899')

Thus, mpmath can easily plot the shape of spike (Cas cannot)
>>> plot(ef2, [0, 0.01])