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

[robve]

(04-10-2021 12:05 AM)Albert Chan Wrote:  

(04-08-2021 01:58 PM)Albert Chan Wrote:  I believe there is nothing to worry about.
Integrand suffered from subtraction cancellations.

Amazingly, qthsh (I picked sqrt version, for speed), does not get garbage from this.

lua> f = function(x) n=n+1; return 1-sqrt(2)*cosh(x)/sqrt(cosh(2*x)) end
lua> f01 = function(t) local y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y) end

With default setting of n=6, eps=1e-9, we have:

lua> Q = require 'quad'
lua> n=0; print(Q.qthsh(f01, 0, 1, 1), n) -- Previous Point
-0.6931471805996943       80
lua> n=0; print(Q.qthsh(f01, 0, 1, 0), n) -- Default Zero
-0.6931471805582439       77

Amazing. Confirmed with the C version, reusing previous point:

qthsh a NAN!!!
qthsh b NAN!!!
qthsh a NAN!!!
qthsh b NAN!!!
qthsh a NAN!!!
k=5 err=1.86609e-09 pts=80
qthsh = -0.693147180599693
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
k=6 err=5.67726e-07 pts=139
wp34s = -0.693147748281705

Using zero as default when NaN occurs:

qthsh a NAN!!!
qthsh b NAN!!!
k=5 err=1.80629e-09 pts=77
qthsh = -0.693147180558242
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
k=6 err=4.70357e-12 pts=133
wp34s = -0.693147180559957

Zero as default improves the result. Note that for the wp34s I had to set eps=1e-15 to prevent early termination. With qthsh with eps=1e-12:

qthsh a NAN!!!
qthsh b NAN!!!
k=6 err=1.10518e-14 pts=155
qthsh = -0.693147180559942

- Rob