15

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))
Dipole
  • 873
  • 2
  • 10
  • 23

1 Answers1

32

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

Kirill
  • 11,438
  • 2
  • 27
  • 51
  • Fantastic! Can you recommend one such book were these techniques are outlined? – Dipole Oct 08 '14 at 08:55
  • 2
    @Jack "Accuracy and Stability of Numerical Algorithms" is a good upper-level book. Any introductory textbook will discuss this as well. – Kirill Oct 08 '14 at 11:43
  • I'd like to know whether you used the Wolfram Mathematica to draw this graph.THX:) – xyz Oct 05 '15 at 04:04
  • 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). – a06e Jul 22 '19 at 14:40
  • @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. – Kirill Jul 23 '19 at 19:04