It's possible to compute $\sqrt{a}$, $a>0$, using only addition and
division and to compute $1/b$, $b>0$ through addition and
multiplication by solving the equationsi) $x^2-a=0$
ii) $\frac{1}{x}-b=0$
using Newton's method. Write down the iteration formulas. What is a
good starting value $x_0?$
The iteration formula is $$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, \quad k=1,2,\ldots$$
How do I a choose a good starting value in the case of $f(x)=x^2-a$?
Best Answer
The goal is to compute $\sqrt{a}$ with a small relative error using a small number of iterations.
I will now analyse Newton's method for the equation $$ f(x) = x^2 - a = 0$$ and explain how to construct an excellent initial guess from scratch. I conclude with some general observations.
Let $$e_n = \sqrt{a} - x_n$$ denote the error and let $$r_n = \frac{e_n}{\sqrt{a}}$$ denote the relative error.
Newton's method takes the form $$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{x_n^2 - a}{2x_n} = \frac{x_n^2 + a}{2x_n} = g(x_n) $$ where $g : (0, \infty) \rightarrow \mathbb{R}$ is given by $$ g(x) = \frac{1}{2} \left( x + \frac{a}{x} \right). $$ Any sensible choice of $x_0$ must satisfy $x_0 > 0$ and it is clear that $x_0 > 0$ implies that $x_n > 0$ for all $n$. We have $$ e_{n+1} = \sqrt{a} - x_{n+1} = - \frac{x_n^2 - 2 a x_n + a}{2x_n} = - \frac{(x_n - \sqrt{a})^2}{2x_n} = - \frac{e_n^2}{2 x_n} $$ which implies $$ r_{n+1} = - \frac{\sqrt{a}}{2 x_n} \frac{r_n^2}{2}. $$ An elementary analysis reveals that $$ g(x) \ge \sqrt{a}, \quad x > 0. $$ It follows that $$ x_n \ge \sqrt{a}, \quad n \ge 1, $$ regardless of the choice of $x_0 > 0$. We can therefore assume that $x_n \ge \sqrt{a}$. If necessary, we simply discard $x_0$ and renumber, i.e., $x_n' = x_{n+1}$ for $n \ge 0$. It follows that $$ |r_{n+1}| \leq \frac{1}{2} r_n^2. $$ Let $y_n = \frac{1}{2} |r_n|$, then $$ y_{n+1} \leq y_n^2 $$ By induction on $n$ we discover, that $$ y_n \leq \left(y_0\right)^{2^n} $$ or equivalently $$ |r_n| \leq 2 \left(\frac{|r_0|}{2} \right)^{2^n}. $$ We conclude that Newton's method converges quadratically, if the initial relative error is less than $2$.
We now turn to the question of constructing a good initial guess. This is not easy for a general $a > 0$. However, the general case can be reduced to the special case of $a \in [1,4)$. To this end, we consider the binary representation of $a > 0$. We have $$ a = f \cdot 2^m,$$ where $f \in [1,2)$ and $m$ is an integer. Let $m = 2k + r$, where $r \in \{0,1\}$. Then, $$ \sqrt{a} = (\sqrt{f \cdot 2^r}) 2^k.$$ We conclude that any square root can be computed provided that we have the ability to compute the square root of any number in the interval $[1,4]$.
An initial guess can be constructed form the best uniform approximation of the square root on this interval. It is easy to verify that the function $h$ given by $$ h(s) = \sqrt{s} - \left(\frac{1}{3} s + \frac{17}{24}\right) $$ satisfies $$ \forall \: s \in [1,4] \: : \quad -\frac{1}{24} \leq h(s) \leq \frac{1}{24}. $$ The initial guess $$ x_0(a) = \frac{1}{3} a + \frac{17}{24} $$ will ensure that the initial residual satisfies $$ |r_0| \leq \frac{1}{24}. $$ In this case, we find that $n=4$ is the smallest integer such that $$ |r_n| \leq u, $$ where $u = 2^{-53} \approx 10^{-16}$ is the unit roundoff in IEEE double precision floating point arithmetic.
Be wary of the difference between exact and computer arithmetic. In exact arithmetic, Newton's method applied to the equation $$ f(x) = x^2 - a = 0$$ will converge for any initial guess $x_0 > 0$. This is not true in floating point arithmetic. The product $x \cdot x$ needed for the computation of $f(x)$ will overflow if $$x > \sqrt{\Omega},$$ where $\Omega$ is the largest representable number. The initial guess must therefore satisfy $$x_0 \leq \sqrt{\Omega}$$ or your implementation will fail immediately. This significantly reduces the choice of $x_0$. In particular, $x_0 = a$ is not a good choice for large values of $a$.
It is not irrelevant how you implement Newton's method. The expressions $$ x_{n+1} = x_n - \frac{x_n^2 - a}{2x_n}$$ and $$ x_{n+1} = \frac{x_n^2 + a}{2x_n}$$ and $$ x_{n+1} = \frac{1}{2} \left( x_n + \frac{a}{x_n} \right)$$ are equivalent in exact arithmetic, but not in floating point arithmetic. Here the first expression is of the preferred form, i.e., a good approximation plus a small correction. It does not matter that $x_n^2 - a$ suffers from catastrophic cancellation when $x_n \approx \sqrt{a}$. If $x_n$ is a good approximation, then you do not need to compute the correction with a tiny relative error. The first expression also has the virtue of utilizing the residual $x_n^2 - a$. This number is often used to determine if the iteration has converged.
Even in exact arithmetic it is not true that avoiding points where $f'(x_0) \not = 0$ is sufficient to ensure convergence. An example is the equation $$f(x) = \sin(x) = 0, \quad |x| < \frac{1}{2} \pi.$$ On this interval $f$ has one zero namely $x = 0$ and $f'(x) = \cos(x)$ has no zeros. Newton's method takes the form $$x_{n+1} = x_n - \tan(x_n).$$ Now if $x_0 = \nu$, where $\nu$ solves the equation $$ \tan(x) = 2x $$ then $$ x_1 = - \nu, \quad x_2 = \nu$$ and we are back where we started. It is easy to verify that $\nu \approx 1.1656$ is one such choice.