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))

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

  • $\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

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.