Estimate 1 / tanh (x) - 1 / x for very small x
I need to calculate the value
1/tanh(x) - 1/x
for x > 0
, where it x
can be both very small and very large.
Asymptotically for small x
we have
1/tanh(x) - 1/x -> x / 3
and for large x
1/tanh(x) - 1/x -> 1
Anyway, when evaluating an expression, already with 10^-7
and less rounding errors, output an expression that will evaluate to exactly 0:
import numpy
import matplotlib.pyplot as plt
x = numpy.array([2**k for k in range(-30, 30)])
y = 1.0 / numpy.tanh(x) - 1.0 / x
plt.loglog(x, y)
plt.show()
source to share
For very small x
one can use the Taylor extension 1/tanh(x) - 1/x
around0
,
y = x/3.0 - x**3 / 45.0 + 2.0/945.0 * x**5
The error has an order of magnitude O(x**7)
, so if the breakpoint is selected 10^-5
, the relative and absolute error will be well below the machine precision.
import numpy import matplotlib.pyplot as plt x = numpy.array([2**k for k in range(-50, 30)]) y0 = 1.0 / numpy.tanh(x) - 1.0 / x y1 = x/3.0 - x**3 / 45.0 + 2.0/945.0 * x**5 y = numpy.where(x > 1.0e-5, y0, y1) plt.loglog(x, y) plt.show()
source to share
Use python package mpmath
for arbitrary decimal precision. For example:
import mpmath
from mpmath import mpf
mpmath.mp.dps = 100 # set decimal precision
x = mpf('1e-20')
print (mpf('1') / mpmath.tanh(x)) - (mpf('1') / x)
>>> 0.000000000000000000003333333333333333333333333333333333333333311111111111111111111946629156220629025294373160489201095913
It becomes extremely accurate.
mpmath
Have a look at plotting . mpmath
works well with matplotlib
which one you are using, so this should fix your problem.
Here's an example on how to integrate mpmath into the code you wrote above:
import numpy
import matplotlib.pyplot as plt
import mpmath
from mpmath import mpf
mpmath.mp.dps = 100 # set decimal precision
x = numpy.array([mpf('2')**k for k in range(-30, 30)])
y = mpf('1.0') / numpy.array([mpmath.tanh(e) for e in x]) - mpf('1.0') / x
plt.loglog(x, y)
plt.show()
source to share
Probably a simpler solution to overcome this is to change the data type in which numpy is running:
import numpy as np import matplotlib.pyplot as plt x = np.arange(-30, 30, dtype=np.longdouble) x = 2**x y = 1.0 / np.tanh(x) - 1.0 / x plt.loglog(x, y) plt.show()
Usage longdouble
as the datatype gives the correct solution without rounding errors.
I modified your example, in your case the only thing you need to change is:
x = numpy.array([2**k for k in range(-30, 30)])
in
x = numpy.array([2**k for k in range(-30, 30)], dtype=numpy.longdouble)
source to share