Why is computation of this function numerically unstable?

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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.