[wp34s] New integration program
04-04-2021, 03:33 PM (This post was last modified: 04-10-2021 10:45 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 2,231 Joined: Jul 2018
RE: [wp34s] New integration program
(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 [-∞,+∞]

---

(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)
(mpf('0.81659478386385076'), mpf('2.0486468077507425e-10'), 161, 5, 0)

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