1
$\begingroup$

Naive evaluation of $\sqrt{a + x} - \sqrt{a}$ when $|a| >> |x|$ suffers from catastrophic cancellation and loss of significance.

WolframAlpha gives the Taylor series for $\sqrt{a+x}-\sqrt{a}$ as: $$\frac{x}{2 \sqrt{a}} - \frac{x^2}{8 a^{3/2}} + \frac{x^3}{16 a^{5/2}} - \frac{5 x^4}{128 a^{7/2}} + \frac{7 x^5}{256 a^{9/2}} + O(x^6)$$ which (I think) equals: $$\sqrt{a} \left( \frac{1}{2} \left(\frac{x}{a}\right) - \frac{1}{8} \left(\frac{x}{a}\right)^2 + \frac{1}{16} \left(\frac{x}{a}\right)^3 - \frac{5}{128} \left(\frac{x}{a}\right)^4 + \frac{7}{256} \left(\frac{x}{a}\right)^5 + O\left(\left(\frac{x}{a}\right)^6\right) \right)$$

How quickly do the coeffients decrease?

How many terms are needed to reach $53$ bits of accuracy (IEEE double precision) in the result given that $10^{-300} < \left|\frac{x}{a}\right| < 1$ is known?

Alternatively, what are the threshold values of $\left|\frac{x}{a}\right|$ where the number of terms changes?

What about rounding errors, assuming each value is stored in double precision?

$\endgroup$
  • $\begingroup$ @Winther shouldn't that be $x / \left(\sqrt{a + x} + \sqrt{a}\right)$? thanks, investigating... $\endgroup$ – Claude May 13 at 19:00
2
$\begingroup$

To avoid cancellation error the first thing to do is to write:

$$ \sqrt{a+x}-\sqrt{a}=\frac{x}{\sqrt{a+x}+\sqrt{a}}=\sqrt{a}\frac{x}{a} \frac{1}{1+\sqrt{1+\frac{x}{a}}} $$

then with $y=\frac{x}{a}$ you must approximate this $$ \sqrt{a}\frac{y}{1+\sqrt{1+y}} $$ fonction for $y\in[10^{-300},1]$. This function has nothing pathological and IMHO can be computed in a straightforward way.

If you really want to use Taylor series for $y\sim 0$ $$ \sqrt{a}\frac{y}{1+\sqrt{1+y}}=\sqrt{a}(\frac{y}{2}-\frac{y^2}{8}+\frac{y^3}{16}-\frac{5 y^4}{128}+\frac{7 y^5}{256}+O\left(y^6\right)) $$ I assume that the series is alternating, hence the error term $e$ is majored by $|e|<\sqrt{a}\frac{7y^5}{256}$. For instance if you want $|e|<10^{-q}$ you can use the Taylors series for $0\le y \le y_*$ where $y_*$ is such that $$\sqrt{a}\frac{7y_*^5}{256}<10^{-q}$$ which gives $$y_*<10^{-q/5}(\frac{256}{7\sqrt{a}})^{1/5}$$

Example: with $q=5$, $a=3$

We get $y_*<0.184042$.

That means that you can use the Taylors series for $ y_*=\frac{x_*}{a}<0.184042$, hence $x_*<3\times 0.184042 \approx 0.552125$.

Let's try with $x=0.55$.

With the initial formula we find: $$ \sqrt{a+x}-\sqrt{a}\approx 0.152094 $$

With the Taylor series, with $y=\frac{0.55}{3}$ we get $$ \sqrt{a}(\frac{y}{2}-\frac{y^2}{8}+\frac{y^3}{16}-\frac{5 y^4}{128}+\frac{7 y^5}{256})\approx 0.152095 $$

We see that the error $|e|=|0.152094-0.152095|\approx 1.17957\times 10^{-6}$ is less that $10^{-q}=10^{-5}$ as expected

$\endgroup$
1
$\begingroup$

The Taylor series is

$$ \sqrt{a+x} - \sqrt{a} = \sum_{k=1}^\infty (-1)^{k+1} \frac{(2k)!}{(k!)^2(2k-1)} 4^{-k} a^{1/2-k} x^k$$ If $|x/a| < 1$, the absolute values of the terms decrease, since if $c_k = (2k)!/((k!)^2 (2k-1) 4^k)$, $$ \frac{c_{k+1}}{c_k} = \frac{2k-1}{2k+2} < 1$$ Thus if $a > x > 0$ the absolute value of the error is always less than that of the next term. However, if $x/a$ is close to $1$ the convergence is rather slow: $$c_k \sim \frac{1}{2 \sqrt{\pi} k^{3/2}}$$ so that won't be less than $2^{-53}$ unless $k > 1.862 \times 10^{10}$ approximately.

$\endgroup$
0
$\begingroup$

It was pointed out by Robert Israel that the series does badly when $|x| \approx |a|$, but in that case the loss of significance of the naive evaluation is small.

It was also suggested by Winther (and a since-deleted answer) to rewrite as $$\frac{x}{\sqrt{a+x}+\sqrt{a}}$$ The series for the denominator is similar to the series in the question, only with a leading constant term. This means that when $\left|\frac{x}{a}\right|$ is small enough, the addition of terms eventually becomes insignificant in double arithmetic.

If $\left|\frac{x}{a}\right| < 2^{-52}$, $1$ term is sufficient. Otherwise

If $\left|\frac{x}{a}\right| < 2^{-25}$, $2$ terms are sufficient. Otherwise

If $\left|\frac{x}{a}\right| < 2^{-16}$, $3$ terms are sufficient. Otherwise

If $\left|\frac{x}{a}\right| < 2^{-11}$, $4$ terms are sufficient. Otherwise

If $\left|\frac{x}{a}\right| < 2^{-9}$, $5$ terms are sufficient. Otherwise

If $\left|\frac{x}{a}\right| > 2^{-9}$, the loss of precision in the addition $a + x$ is relatively small.

But in fact, perhaps $$\frac{x}{\sqrt{a + x} + \sqrt{a}}$$ evaluated in double precision is good enough for all $|x| << |a|$ and series are unnecessary?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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