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

[Albert Chan]

(04-18-2021 12:56 PM)robve Wrote:  There is one bad case exp(-3*(x*x)-4/(x*x)) with 228 more points. The predictor pushes too hard on this one where the diff is close to zero already at the begin of the search, meaning there is no need to search for d in this case.

Note: this quote was removed, but it is instructive to see why it failed ...

Close to zero at the beginning of search is expected, for any f(x).
When v=1, f(v)*v - f(1/v)/v = f(1) - f(1) = 0

This is just the trivial solution, and cannot be used to predict d = 1.

To simplify, let v = e^t, instead of v = 2^t:

\(\displaystyle \int_0^∞ f(x)\,dx
= \int_{-∞ }^∞ f(v)·v\,dt
= \int_0^∞ \left[f(v)·v + f({1\over v})·{1\over v}\right]\,dt
\)

Except for v=1, if last 2 term are same size, we assumed both will be tiny.
This should covered most of the (likely bell-shaped) area, thus can estimate optimal d.

---

Back to why f(x) = exp(-3*(x*x)-4/(x*x)) failed sign change search.

f(v)*v = f(1/v)/v, where v = e^t

Take log both side:

-3*exp(2t) - 4*exp(-2t) + t = -3*exp(-2t) - 4*exp(2t) - t
2t + exp(2t) - exp(-2t) = 0
t + sinh(t) = 0

On the left, both terms are increasing function, only solution is t=0.
Thus, except for the trivial v=1, there is no solution.

And, there are many more bad cases, say, g(x) = exp(-m*x), m ≥ 1

g(v)*v = g(1/v)/v

Again, substitute v=e^t, and take log both side:

-m*e^t + t = -m*e^-t - t
2t = m*(e^t - e^-t)
t = m * sinh(t)

sinh(t)    = t + t^3/3! + t^5/5! + ...
sinh(t)/t = 1 + t^2/3! + t^4/5! + ... ≥ 1

Removed trivial solution, (t = sinh(t) ⇒ t = 0), we have m < 1

---

Perhaps estimate of d based on peak of f(v)*v is safer ?
(there is always a peak, unless integral diverges)

f(x) = exp(-3*(x*x)-4/(x*x))
Searching for peak of f(v)*v (we can also do numerically):

diff(f(v)*v, v) = 0
f(v) + v*f'(v) = 0
f(v) + v*f(v)*(-6v + 8/v^3) = 0
1 - 6v^2 + 8/v^2 = 0
6v^4 - v^2 - 8 = 0                     // assumed v ≠ 0
v^2 = (1+√193)/12                   // we throw away the negative root, v^2 > 0
v = √((1+√193)/12) ≈ 1.114     // again, throw away negative root, v=e^t > 0

lua> Q = require 'quad'
lua> f = function(x) x=x*x; return exp(-3*x-4/x) end
lua> Q.quad(Q.count(f), 0, huge, 1.114)
0.0005013071292635826       1.4667498799376226e-010
lua> Q.n
41

Same treatment for g(x) = exp(-m*x):

diff(g(v)*v, v) = 0
g(v)*1 + (g(v)*-m)*v = 0
1 - m*v = 0
v = 1/m       -- when g(v)*v peaked

Peak for d might not be "optimal", but may be good enough.