Roots of $x^2-2ax+b=0,~a^2>>b$, stable algorithm

Question

Find a stable algorithm for the computation of the roots of the equation $x^2-2ax+b=0$ if $a,~b>0$ and $a^2>>b$.

Attempt. The roots are $x_1=a-\sqrt{a^2-b}$ and $x_2=a+\sqrt{a^2-b}$, where we have a cancelation error for the computation of $x_1$, since $a^2-b\cong a^2$. So we compute $\displaystyle x_1=\frac{b}{a+\sqrt{a^2-b}}$ (stable algorithm).

My question: is it correct to say that $\displaystyle \frac{b}{2a}$ is a "good" approximation for $x_1$? Or should we remain on the algorithm $\displaystyle x_1=\frac{b}{a+\sqrt{a^2-b}}$?

Thanks in advance!

Answer 1

• About one part of your question, no purely mathematical answer can be given because $\displaystyle x_1=\frac{b}{a+\sqrt{a^2-b}}$ is strictly the same as $a-\sqrt{a^2-b}.$ It is a floating point issue.

• About the other part, yes, $\frac{b}{2a}$ is a good approximation.

This can be attested by applying Taylor expansion :

$$\sqrt{1-x}=1-\frac{1}{2}x-\frac{1}{8}x^2-\cdots$$

in the following way:

$$x_1=a-\sqrt{a^2-b}=a-a\sqrt{1-\frac{b}{a^2}}=a-a\left(1-\frac{1}{2}\frac{b}{a^2}-\frac{1}{8}\frac{b^2}{a^4}- \cdots\right)$$

Giving :

$$x_1 \ \approx \ \frac{b}{2a}+\frac{b^2}{8a^3}$$

Thus, the order of the error is approximately: $\frac{b^2}{8a^3}$,

which is very small because $b << a^2.$

Answer 2

$$\frac1{x_1}=\frac{a+\sqrt{a^2-b}}{b}=\frac ab+\frac ab\sqrt{1-\frac b{a^2}}\approx\frac ab{\left(2-\frac b{2a^2}\right)}$$ so that the relative truncation error is indeed small.

You can safely ignore the correction when $\dfrac b{a^2}<4\text{ ulp}$, (https://en.wikipedia.org/wiki/Unit_in_the_last_place). For higher ratios, it will depend on your own accuracy requirements.