[wp34s] New integration program

[Albert Chan]

(03-24-2017 05:12 PM)emece67 Wrote:  \[\int_{-1}^\infty{\ln(x+1)\over x^2+1}dx={3\pi\ln(2)\over8}\approx0.8165947838638507989377583368391052\ldots\]

With this integral transformed, I noted some weekness of tanh-sinh.
(Or, it might be just a fluke with this example ...)

Rewrite above integral, then let x=e^y:

\(\displaystyle \int _0 ^∞ {\ln(x) \over 1+(x-1)^2} dx = \int _{-∞} ^∞ {y\;e^y \over 1+(e^y-1)^2} dy\)

>>> from mpmath import *
>>> import double_exponential # emece67 source from here
>>> quadde = lambda f,(a,b): double_exponential.double_exponential(f,a,b)

>>> f = lambda y: y*exp(y) / (1 + expm1(y)**2)
>>> quadde(f, [-inf,inf]) # method sinh-sinh
(mpf('0.81659478386382767'), mpf('7.7731406533665393e-8'), 73, 4, 2)

To use tanh-sinh, just convert back intervals: [-∞,+∞] -> [0,1]

>>> def f01(t): y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y)
...
>>> quadde(f01, [0,1]) # method tanh-sinh
(mpf('0.81659478390520546'), mpf('1.0424940875997102e-5'), 77, 5, 0)

It seems sinh-sinh does better than tanh-sinh, for interval [-∞,+∞]

---

>>> quadde(f, [-40,40])
(0, mpf('0.83732598514658618'), 87, 5, 0)

A bit unexpected, getting area of 0.0, with huge error estimate.

plot(f, [-40,40]), showed it has a shape of a heartbeat (inside ±5)
We should instead move dominant part to the edge.

>>> fold = lambda x: f(x) + f(-x)
>>> quadde(fold, [0,40])
(mpf('0.81659478386385076'), mpf('2.0486468077507425e-10'), 161, 5, 0)

Result is similar to exp-sinh, for [0,+∞]

>>> quadde(fold, [0,inf])
(mpf('0.81659478386385087'), mpf('4.1397899552819695e-9'), 131, 4, 1)

Edit: I misinterpreted last argument, should be 0=tanh-sinh, 1=exp-sinh, 2=sinh-sinh