13
$\begingroup$

I am having trouble implementing a function numerically. It suffers from the fact that at large input values the result is a very large number times a very small number. I am not sure if catastrophic cancellation is the correct term so please correct me if it is. Evidence of something going awry:

enter image description here

How can I avoid the oscillations and the assigment of 0.0 for larger inputs of 6?

Here is my function:

import numpy as np

def func(x):
    t = np.exp(-np.pi*x)
    return 1/t*(1-np.sqrt(1-t**2))
$\endgroup$
31
$\begingroup$

This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$ Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that formula in its standard form is numerically unstable for some parameter values). The way to get rid of $1-\sqrt{1-t^2}$ is to multiply and divide by $1+\sqrt{1-t^2}$. Here is a comparison:

enter image description here

$\endgroup$
  • $\begingroup$ Fantastic! Can you recommend one such book were these techniques are outlined? $\endgroup$ – Dipole Oct 8 '14 at 8:55
  • 2
    $\begingroup$ @Jack "Accuracy and Stability of Numerical Algorithms" is a good upper-level book. Any introductory textbook will discuss this as well. $\endgroup$ – Kirill Oct 8 '14 at 11:43
  • $\begingroup$ I'd like to know whether you used the Wolfram Mathematica to draw this graph.THX:) $\endgroup$ – xyz Oct 5 '15 at 4:04
  • $\begingroup$ Do you know of any references collecting and/or discussing similar tricks to rewrite mathematical expressions in mathematically equivalent ways which reduce loss of significance? I read the book by Higham, but the discussion is general and all the later chapters are about linear algebra (which is not my topic at the moment). $\endgroup$ – becko Jul 22 at 14:40
  • $\begingroup$ @becko It's pretty ad hoc in my experience. It's much easier to do if you have a way of testing your formula with correct answers (even if you just generate them with extra-precision arithmetic), so that you don't go looking for numerical instability without first having failing test cases. And if it works for all known inputs, there's no real problem whether numerical instability is present anywhere or not. $\endgroup$ – Kirill Jul 23 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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