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

[Albert Chan]

(04-12-2021 03:06 PM)robve Wrote:  

(04-12-2021 12:20 PM)Albert Chan Wrote:  To reduce noise, we should quit summing, if we actually hit the finite edge.

Good point! It makes a lot more sense to stop Exp-Sinh at this point when x hits c for x=c+d/r and avoid using a like qthsh. The weight gets very small when approaching the finite endpoint.

More importantly, infinite "edge" has weight that gets very very big.

For exp-sinh, with limits of [c, d*inf], q before scaled up by cosh(j*h):

q = f(c + d/r)/r + f(c + d*r)*r

When r is huge, second term might generate magnified noise, messing up the sum.

Example, ∫(x^-0.99, x=0..1)

\(\displaystyle \int_0^1 x^{-0.99}\, dx
= \int_0^∞ (e^{-y})^{-0.99} (e^{-y}\, dy)
= \int_0^∞ e^{-0.01y} \, dy
\)

If we apply transformation without last simplifed step, integrand generate much noise.

>>> Q = require 'quad'
>>> f = function(x) return x^-0.99 end
>>> f1 = function(y) y=exp(-y); return f(y)*y end

lua> Q.quad(Q.count(f), 0, 1), Q.n
99.92043566319589       654
lua> Q.quad(Q.count(f1), 0, huge), Q.n
99.93721615545337       459

Now, try simplified form, avoiding the noise

lua> f2 = function(y) return exp(-0.01*y) end
lua> Q.quad(Q.count(f2), 0, huge), Q.n
99.9999999878134        237

Nice, but we can do better ! (see previous post)

exp(0) = 1
exp(-0.01y) = 0.001       → -y/100 = ln(10^-3) = -3*ln(10) ≈ -3*2.3 ≈ -7

\(\int_0^∞ = \int_0^{700} + \int_{700}^∞ :\)

lua> Q.quad(Q.count(f2), 0, huge, 700), Q.n
100        71