14
$\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))
Improve this question
$\endgroup$

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
5
  • $\begingroup$ Fantastic! Can you recommend one such book were these techniques are outlined? $\endgroup$
    – Dipole
    Commented Oct 8, 2014 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
    Commented Oct 8, 2014 at 11:43
  • $\begingroup$ I'd like to know whether you used the Wolfram Mathematica to draw this graph.THX:) $\endgroup$
    – xyz
    Commented Oct 5, 2015 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$
    – a06e
    Commented Jul 22, 2019 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
    Commented Jul 23, 2019 at 19:04

Your Answer

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

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