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

[robve]

(04-12-2021 12:20 PM)Albert Chan Wrote:  

robve Wrote:I now see what the difference is between qthsh and quad, it's quite simple: qthsh uses f(a+x) while quad uses f(c-x). When a=0, the check a+x>a always succeeds for nonzero positive x, but c-x>a fails due to cancellation with c=(a+b)/2.

This means that qthsh approaches the zero point aggressively, something we already saw before. It also means that qthsh is not "symmetric" in the bounds, something we also saw before. This could be desirable actually, when functions exhibit interesting behavior close to x=0.

Yes, it is more accurate interpolate from the closest edge.
https://www.hpmuseum.org/forum/thread-16...#pid140084

Let me add to my PM some context: a notable difference are the function evaluation guards at both endpoints in the different implementations which are a+d*r>a in qthsh versus c-d*r and b-d*r<b in qthsh versus c+d*r<b. If the finite endpoints a or b are zero or close to zero, then the corresponding endpoint may be more closely approached by qthsh with additional abscissas. This is important. More accurate results can be achieved with qthsh when the integrand has a significant area close to \( x=0 \), for example \( \int_0^1 1/\sqrt{x}\, dx \) gives:
qthsh points=63 relative error=1e-13
quad (old version using x=c+d*r similar to WP-34S) points=115 relative errror=6e-9

Likewise, WP-34S "pure" Tanh-Sinh (C version with IEEE 754) with eps=1e-15 has a high relative error of 5e-7 after 25 points "convergence" at k=2. With eps=1e-20 WP-34S has a 2e-10 relative error after 217 points and k=5 convergence.

qthsh converged quickly at level k=3 whereas quad using x=c+d*r converged at level k=4 with a larger error, requiring one more level in an attempt to produce a result within eps=1e-9. This is illustrated by the following chart (click to enlarge):

   

This shows the point distributions of qthsh and quad for this integrand, using a log scale to emphasize the points approaching the zero endpoint a=0.

(04-12-2021 12:20 PM)Albert Chan Wrote:  Even if interpolate closed to center !
(c=(a+b)/2 might have rounding errors; a, b are user inputs, thus can considered exact)

And, accuracy affect both edges, not just the "zero" side.
Example, (code snippet from previous post)

Yes, c can cause rounding errors and loss of precision in a+b. Whenever we add or subtract there can be a (significant) loss. In the end, IEEE 754 double precision does not help. When c is used we should also compare to c (at least in this case). To do so in the code, the Exp-Sinh and Sinh-Sinh quadratures should be split from the Tanh-Sinh computation, which I've already done in yesterday's updated document (that will continue to evolve.)

(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.

BTW. one trick to get a and b numerically corrected to prevent noise with c:
c=(a+b)/2
d=(b-a)/2
a=c-d
b=c+d
This helps, but it is "iffy" to use for this algorithm, because it (slightly) changes the endpoints! In other contexts these tricks can be important, for example in numerical differentiation e.g. h=1e-3, x=a-h then h is NOT the difference between a and x exactly, but the adjusted h=a-x is.

Note that the Exp-Sinh and Sinh-Sinh quadratures described in the document use a Michalski & Mosig simplification to drop \( \pi/2 \). These rules are transformations by substitution anyway. Having a few more points in the area between the endpoints looks advantageous, judging from the 818 integrals tested in the spreadsheet compared to a "pure" Tanh-Sinh rule. Other parameterizations are possible to control the point distributions, e.g. Michalski & Mosig in their paper start with h=1.5 instead of 1.

An update of the document and code is forthcoming.

- Rob