Why is computing the function at certain ranges of h $$f(x) = \dfrac{\sqrt{-x+a} - \sqrt{2x+a}}{ 4a}$$ unstable? How would I rewrite this function so that it is stable?

  • 3
    $\begingroup$ I am not following. h = 0 would be a divide by zero? $\endgroup$ – paparazzo Feb 6 '17 at 16:54
  • 10
    $\begingroup$ I am not following. What's $h$? $\endgroup$ – David Richerby Feb 6 '17 at 19:05
  • 6
    $\begingroup$ I get that $h$ is the constant, but now your update introduced another variable and invalidated already existing answers - that is bad. $\endgroup$ – Evil Feb 6 '17 at 19:40
  • 5
    $\begingroup$ There's no $h$ in this function at all. What are you talking about? Do you mean $a$? $\endgroup$ – Wildcard Feb 6 '17 at 22:40
  • 2
    $\begingroup$ Cross-posted: math.stackexchange.com/q/2131244/14578, cs.stackexchange.com/q/69931/755. Please do not post the same question on multiple sites. Each community should have an honest shot at answering without anybody's time being wasted. $\endgroup$ – D.W. Feb 7 '17 at 2:13

The question has been changed since this answer.

Another alternative is $$\begin{align}f(x)&=\frac{\sqrt{x+h}-\sqrt{x-h}}{2h}\\ &=\frac{\sqrt{x+h}-\sqrt{x-h}}{2h}\cdot\frac{\sqrt{x+h}+\sqrt{x-h}}{\sqrt{x+h}+\sqrt{x-h}}\\ &=\frac{1}{\sqrt{x+h}+\sqrt{x-h}}\end{align}$$

Now, even if $h$ becomes very small, the value approaches $$\frac{1}{2\sqrt{x}}$$ as it should.

| cite | improve this answer | |
  • $\begingroup$ Hey, there. I see you've just started reviewing posts to the site. Please actually read the post before clicking one of the action buttons. You reviewed four posts in 36 seconds and that just isn't possible. In my opinion, three of those questions had serious problems with them, yet you clicked "no action needed" without even reading them. Please do not do that: it is harmful to the site. $\endgroup$ – David Richerby Mar 3 '17 at 18:41
  • $\begingroup$ @DavidRicherby Thank you for the feedback. I did read them, however, I genuinely thought that no action was needed. Perhaps I am not clear enough of the site's policies yet. I will refrain from reviewing till I get more familiar with the do's and don't's. I do not wish to cause any harm. Thank you again. $\endgroup$ – GoodDeeds Mar 3 '17 at 18:48
  • $\begingroup$ Thanks -- appreciated. Do feel free to ask in chat or meta if you need help. $\endgroup$ – David Richerby Mar 3 '17 at 18:58

What you are experiencing is called 'loss of significant digits'. As h gets smaller, the two square-root terms become closer and closer together, in the sense that their values start with more and more matching digits. Since floating-point cannot represent more than a finite prefix of the infinite-digit actual value for each term, the result of subtracting the terms eventually loses all meaningful information.

One way to rewrite this is to expand the two terms as Taylor series and subtract them symbolically; some terms will cancel, the rest will yield a useable series for evaluation (probably more stable, not necessarily more accurate).

| cite | improve this answer | |

We can write it as $$f(x) = \dfrac{\sqrt{-x+a} - \sqrt{2x+a}}{ 4a}= \dfrac{\sqrt{-\frac xa+1} - \sqrt{\frac{2x}a+1}}{ 4\sqrt a}$$ If $x \ll a$ both of the square roots will be just about $1$. When you add the $\frac xa$ terms to $1$ you will lose many bits of $\frac xa$, then when you subtract the square roots the $1$s will cancel. If your floats have $53$ bits in the mantissa and $\frac xa$ is about $2^{-30}$ you will only keep about $23$ bits of $\frac xa$. Only the top $23$ bits of the numerator will be accurate. You need to analytically cancel the $1$s to get the accuracy you want. You can do that with the Taylor expansion $$\sqrt {-\frac xa+1} \approx 1-\frac x{2a}, \sqrt {\frac {2a}x+1} \approx 1+\frac xa\\\sqrt{-\frac xa+1} - \sqrt{\frac{2x}a+1}\approx -\frac {3x}{2a}$$ Now you don't lose the precision.

| cite | improve this answer | |
  • $\begingroup$ Stopping at the first term may also lead to some issues; you should at least generate a numerically stable error term, so you'd be aware if you have tiny error or not. $\endgroup$ – Yakk Feb 7 '17 at 14:39
  • $\begingroup$ @Yakk: I think this answers the question of why there is a problem. The answer is stable, if perhaps not accurate. Yes, to make a good algorithm one would need to check the magnitude of $\frac xa$ and decide if more terms were needed. $\endgroup$ – Ross Millikan Dec 4 '17 at 3:15

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