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

[robve]

(04-03-2021 03:17 AM)emece67 Wrote:  These are the results for my Python script:

Thank for sharing. I rewrote the Python code today in C with exactly the same functionality and updated the document with the results and other additions. In this way I can compare apples to apples using IEEE floating point.

I noticed that I had to set eps=1e-15 in the C code I wrote based on your Python implementation, because the error threshold variable with eps=1e-9 is only thr=3e-4, returning results with large errors even for the easy integrals. I also used n=6 levels to make the comparison fair.

The values with the C implementation with eps=1e-15 are the same as your numbers, except when the number of points is larger due to my choice of n=6 levels:

Code:

# (points, est.err)
1 (47,8e-12)
2 (49,5e-12)
3 (49,1e-8+)
4 (49,2e-13)
5 (49,7e-12)
6 (49,5e-8)
7 (49,2e-12)
8 (25,5e-7)
9 (25,5e-7)
10 (395,2e-4)
11 (395,2e-4)
12 (49,8e-8)
13 (395,2e-4)
14 (25,4e-7)
15 (25,4e-7)
16 (335,1e-16)
17 (49,7e-10)
18 (49,2e-12)
19 (395,1e-3)
20 (395,6e-3)
21 (395,6e-3)

Integral 16 is the only one that differs for a different reason, because it converges with 137 points in your case, but with IEEE double precision and eps=1e-15 in the C version, it takes 335 points.

I noticed a few things in the Python code, including:

      try:
        fpl = f(bpa2 + bma2*r)
      except Exception:
        fpl = 0;
      p = fpl*w
      try:
        fmi = f(bpa2 + bma2/r) if expsinh else f(bpa2 - bma2*r)
      except Exception:
        fmi = 0

I may be mistaken, but this sets fpl and fmi to zero when an exception occurs. This may put some downward pressure on the quadrature towards zero and could lead to delayed convergence if the weight isn't yet insignificant, in turn requiring more points that in turn could cause more exceptions in the same area. I did some comparisons and it looks like interpolating/extrapolating the missing function values is best (rather than terminating the iterations). A simple way to interpolate is to reuse the previous fpl and fmi, although more elaborate schemes could be used. Think about functions that have a large nonzero limit at an endpoint.

On the other hand, the Python code does not get as close to the endpoints (due to tmax) as the other implementations. By comparison, the VB implementation uses a fixed tmax=6.56, which I am looking into.

- Rob