# How do I globally change the precision of a piece of code in Python to debug it?

I am solving a system of non-linear equations using the Newton-Raphson method in Python. This involves using the solve(Ax,b) function (spsolve in my case, which is for sparse matrices) iteratively until the error or update reduces below a certain threshold. My specific problem involves calculating functions such as $$x/(e^x - 1)$$, which are badly calculated for small $$x$$ by Python, even using np.expm1().

Despite these difficulties, it seems like my solution converges, because the error becomes of the order of $$10^{-16}$$. However, the dependent quantities, do not behave physically, and I suspect this is due to the precision of these calculations. For example, I am trying to calculate the current due to a small potential difference. When this potential difference becomes really small, this current begins to oscillate, which is wrong, because currents must be conserved.

I would like to globally increase the precision of my code, but I'm not sure if that's a useful thing to do since I am not sure whether this increased precision would be reflected in functions such as spsolve. I feel the same about using the Decimal library, which would also be quite cumbersome. Can someone give me some general advice on how to go about this or point me towards a relevant post?

• Have you checked your solution when the error is larger? For example when it is in the range of $10^{-4}$ or $10^{-5}$ for example? You don't need to put the effort of calculating $\frac{x}{e^{x}-1}\Big|_{x=0}$ on Python when you know the answer as it must be 1. You can just create a function and exclude $x=0$ if that's the problem, which I doubt it. Also a quick comment about this argument "because currents must be conserved." The current is not conserved. The charge must be conserved. – Alone Programmer Dec 19 '19 at 3:51

In this answer, I recommended using mpmath Python library for arbitrary precision.