5
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Mithridates the Great
    Dec 19, 2019 at 3:51

1 Answer 1

Reset to default
3
$\begingroup$

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

Notice, that since matrices in mpmath are implemented as dictionaries:

Only non-zero values are stored, so it is cheap to represent sparse matrices.

thus, this particular library seems like a good fit for your purpose of debugging.

I find this library (and approach in general) very useful for debugging purposes, mostly, to convince myself that my stubborn bug doesn't come from the finite precision of floating-point arithmetic. Because if it does, it is an order of magnitude more complicated problem to solve.

Note, that using an arbitrary-precision library will significantly slow down your code; therefore, if you find that your problem is resolved, you still have a huge incentive to come up with a solution using regular double precision.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.