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

[Albert Chan]

(04-06-2021 06:48 AM)emece67 Wrote:  I've just tested re-using the previous fpl or fmi when the exception raises, but the end results are the same.

I preferred original default of 0.0, instead of previous point.

My first post here mentioned problems if u-transformation is applied.
Instead of sharp point, when f(0) or f(1) = inf, it widened the problematic edge.

(03-28-2021 06:59 PM)Albert Chan Wrote:  >>> from mpmath import *
>>> f = lambda x: 1/sqrt(-ln(x))
>>> quadts(f, [0,1], error=True)
(mpf('1.7724538503750329'), mpf('1.0e-10'))
>>>
>>> uf = lambda u: 6*u*(1-u) * f(u*u*(3-2*u))
>>> quadts(uf, [0,1], error=True)
Traceback (most recent call last):
...
ZeroDivisionError

This is not just an issue of u-transformation.
Perhaps integrand already have fuzzy bad edge(s).

Previous point for default will have "upward pressure" for the fuzzy edge.
(this assumed previous point is next to current point; false if calculating points in parallel)

Default of 0.0 will have "downward pressure" for the fuzzy edge.
However, with double-exponential transform, I think more of the bad edge is closer to zero.

Example, using above f, uf, and we apply u-transform to uf, getting uuf.

Run double_exponential() with different behavior for "bad" edge, then square the result (should be pi)
First number is points evaluated, reported from double_exponential() third entry.

Code:

        Default Zero           Previous Point
f   25->3.14159258765585   25->3.14159258765585
uf  47->3.14159262045834   49->3.14159263570079
uuf 39->3.14159264565517  191->3.14160003440478

---

I think setting 0.0 for any exceptions may be overkill.
Example: if f(x) = x^p, and p is not yet defined, quadrature should raised NameError.

Trap ArithmeticError is good enough (= OverflowError, ZeroDivisionError, FloatingPointError)

On the other hand, Inf/NaN can occur without raising exceptions.

>>> ln(0), atanh(1), exp(inf), inf-inf
(mpf('-inf'), mpf('+inf'), mpf('+inf'), mpf('nan'))

We should also test for these, like this patch for double_exponential.py

Code:

      try: y = f(bpa2 + bma2*r)
      except ArithmeticError: y = inf
      p = 0.0 if y-y else y*w

      try: y = f(bpa2 + bma2/r) if expsinh else f(bpa2 - bma2*r)
      except ArithmeticError: y = inf
      if y-y == 0: p += y/w if expsinh else y*w

      p *= ch
      wsl += p
      tnfe += 2