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?
Thank you!
You can try using mpmath, but YMMV. generally scipy uses double precision. For a vast majority of cases, analyzing the sources of numerical errors is more productive than just trying to reimplement everything with higher widths floats.
Can you change the precision globally of a piece of code in Python, as a way of debugging it?