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?
9
-
3$\begingroup$ I am not following. h = 0 would be a divide by zero? $\endgroup$– paparazzoCommented Feb 6, 2017 at 16:54
-
10$\begingroup$ I am not following. What's $h$? $\endgroup$– David RicherbyCommented Feb 6, 2017 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$– EvilCommented Feb 6, 2017 at 19:40
-
5$\begingroup$ There's no $h$ in this function at all. What are you talking about? Do you mean $a$? $\endgroup$– WildcardCommented Feb 6, 2017 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. ♦Commented Feb 7, 2017 at 2:13
|
Show 4 more comments
3 Answers
$\begingroup$
$\endgroup$
3
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.
-
$\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$ Commented Mar 3, 2017 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$ Commented Mar 3, 2017 at 18:48
-
$\begingroup$ Thanks -- appreciated. Do feel free to ask in chat or meta if you need help. $\endgroup$ Commented Mar 3, 2017 at 18:58
$\begingroup$
$\endgroup$
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).
$\begingroup$
$\endgroup$
2
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.
-
$\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$– YakkCommented Feb 7, 2017 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$ Commented Dec 4, 2017 at 3:15