0

What's the general procedure for finding a numerically stable form for a function?

Such as:

$$\frac{1}{1+2x}-\frac{1-x}{1+x}$$ when $x≈0$
or $$ln(x)-1$$when $x≈e$

I've only seen examples, but cannot grasp the "procedure" for what to look for in order to find a numerically stable alternative.

mavavilj
  • 7,270
  • Your examples are not good, the numerically stable form is mostly about removable singularities where the non-cancelled form exhibits catastrophic cancellation, like in $\frac{1-\cos x}{x^2}$, $\frac{x\cos x-\sin x}{x^3}$, $\frac{x^3+2x^2-3}{x^2-4x+3}$,... – Lutz Lehmann Mar 30 '16 at 13:53
  • @LutzL But I found these as examples of numerical instability in my notes. – mavavilj Mar 30 '16 at 13:55
  • 1
    Do your notes explain how these functions exhibit instability, or what the numerically stable forms are? It's not clear to me why these functions would be given as examples. – David K Mar 30 '16 at 13:57
  • This is strange, since the critical points of your formulas are poles where there is no result. $\ln(x)^{-1}$ could be interesting, as well as $x\ln(x)$, both at $x=0$ and close to it. – Lutz Lehmann Mar 30 '16 at 13:59
  • @DavidK I'm trying to find the numerically stable forms. I only know these are numerically unstable. – mavavilj Mar 30 '16 at 14:04
  • One can see that these exhibit numerical instability by plotting them. – mavavilj Mar 30 '16 at 14:17
  • http://www.wolframalpha.com/input/?i=1%2F%281%2B2%2Ax%29+-+%281-x%29%2F%281%2Bx%29+for+-0.1+%3C+x+%3C+0.1 looks OK to me. – David K Mar 30 '16 at 14:25
  • http://www.wolframalpha.com/input/?i=ln%28x%29+-+1+for+2.5%3Cx%3C3 also looks OK. – David K Mar 30 '16 at 14:26
  • @DavidK Is it not a problem if $ln(e)-1=0$, i.e. one could suspect that that WA plot doesn't show the problem? – mavavilj Mar 30 '16 at 15:08
  • To see the problem in WA, use "plot (ln(x) - 1)/(x-exp(1)) for exp(1)-1e-8<x<exp(1)+1e-8" (better use gnuplot or another function plotter locally). The rational expression probably gets transformed internally into one fraction, resolving the numerical instability without showing it. There you really need a dumber resp. more literal function plotter. – Lutz Lehmann Mar 30 '16 at 15:52
  • @LutzL Can you explain why plot that way? – mavavilj Mar 31 '16 at 09:04
  • Because "numerically stable" means that the relative error in the result should be a small multiple of the relative error in the argument. Now $f(x+h)=f(x)+f'(x)h+O(h^2)$, thus for $h^2<μ$ the relative error is $f(x)/(f(x_0)+f'(x_0)(x-x_0))-1$ which for $|f(x_0)|\ggμ$ should be a small multiple of $μ$. However, the relative error for $f(x_0)=0$ reduces to $f(x)/(f'(x_0)(x-x_0))-1$ for $|x-x_0|<\sqrtμ$ and can thus have magnitudes $μ/(x-x_0)$. You get, for instance, in the example relative errors of $\sim 10^{-3}$ for $x-e\sim 10^{-13}$, which is consistent with $μ\approx 10^{-16}$. – Lutz Lehmann Mar 31 '16 at 11:48

0 Answers0