Numerical method for finding the square-root.

Essentially if you are interesting in evaluating $\sqrt{a}$, the idea is to first find the greatest perfect square less than or equal to $a$. Say this is $b^2$ i.e. $b = \lfloor \sqrt{a} \rfloor \implies b^2 \leq a < (b+1)^2$. Then consider the function $$f(x) = b + \dfrac{a-b^2}{x+b}$$ $$f(b) = b + \underbrace{\dfrac{a-b^2}{2b}}_{\in [0,1]} \in [b,b+1]$$ $$f(f(b)) = b + \underbrace{\dfrac{a-b^2}{f(b) + b}}_{\in [0,1]} \in [b,b+1]$$ In general $$f^{(n)}(b) = \underbrace{f \circ f \circ f \circ \cdots f}_{n \text{times}}(b) = b + \dfrac{a-b^2}{f^{(n-1)}(b)+b}$$ Hence, $f^{(n)}(b) \in [b,b+1]$ always.

If $\lim\limits_{n \to \infty}f^{(n)}(b) = \tilde{f}$ exists, then $$\tilde{f} = b + \dfrac{a-b^2}{\tilde{f}+b}$$ Hence, $$\tilde{f}^2 + b \tilde{f} = b \tilde{f} + b^2 + a - b^2 \implies \tilde{f}^2 = a$$

To prove the existence of the limit look at $$(f^{(n)}(b))^2 - a = \left(b + \dfrac{a-b^2}{f^{(n-1)}(b)+b} \right)^2 - a = \dfrac{(a-b^2)(a-(f^{(n-1)}(b))^2)}{(b+f^{(n-1)}(b))^2} = k_{n-1}(a,b)((f^{(n-1)}(b))^2-a) $$ where $\vert k_{n-1}(a,b) \vert \lt1$. Hence, convergence is also guaranteed.

EDIT

Note that $k_{n-1}(a,b) = \dfrac{(a-b^2)}{(b+f^{(n-1)}(b))^2} \leq \dfrac{(b+1)^2 - 1 - b^2}{(b+b)^2} = \dfrac{2b}{(2b)^2} = \dfrac1{2b}$. This can be interpreted as larger the number, faster the convergence.

Comment: This method works only when you want to find the square of a number $\geq 1$.

EDIT

To complete the answer, I am adding @Hurkyl's comment. Functions of the form $$g(z) = \dfrac{c_1z+c_2}{c_3z+c_4}$$are termed Möbius transformations. With each of these Möbius transformations, we can associate a matrix $$M = \begin{bmatrix} c_1 & c_2\\ c_3 & c_4\end{bmatrix}$$ Note that the function, $$f(x) = b + \dfrac{a-b^2}{x+b} = \dfrac{bx + a}{x+b}$$ is a Möbius transformation.

Of the many advantages of the associated matrix, one major advantage is that the associate matrix for the Möbius transformation $$g^{(n)}(z) = \underbrace{g \circ g \circ \cdots \circ g}_{n \text{ times}} = \dfrac{c_1^{(n)} z + c_2^{(n)}}{c_3^{(n)} z + c_4^{(n)}}$$ is nothing but the matrix $$M^n = \begin{bmatrix}c_1 & c_2\\ c_3 & c_4 \end{bmatrix}^n = \begin{bmatrix}c_1^{(n)} & c_2^{(n)}\\ c_3^{(n)} & c_4^{(n)} \end{bmatrix}$$ (Note that $c_k^{(n)}$ is to denote the coefficient $c_k$ at the $n^{th}$ level and is not the $n^{th}$ power of $c_k$.)

Hence, the function composition is nothing but raising the matrix $M$ to the appropriate power. This can be done in a fast way since $M^n$ can be computed in $\mathcal{O}(\log_2(n))$ operations. Thereby we can compute $g^{(2^n)}(b)$ in $\mathcal{O}(n)$ operations.

The iteration to find $\sqrt k$ is $f(x) = \frac{d x+k}{x+d}$ where $d = \lfloor \sqrt k \rfloor$. The iterations start with $x = d$.

If $x$ is a fixed point of this, $x = \frac{d x+k}{x+d}$, or $x(x+d) = dx + k$ or $x^2 = k$, so any fixed point must be the square root.

Now wee see if the iteration increases or decreases. If $y = \frac{d x+k}{x+d}$, $$y - x = \frac{d x+k}{x+d} - x = \frac{d x+k - x(x+d)}{x+d} = \frac{k - x^2}{x+d} $$ so if $x^2 < k$, $y > x$ and if $x^2 > k$, $y < x$.

Also, proceeding like analyses of Newton's method, $y^2-k = \frac{(d x+k)^2}{(x+d)^2} - k = \frac{d^2 x^2 +2 d x k + k^2 - k(x+d)^2}{(x+d)^2} = \frac{d^2 x^2 +2 d x k + k^2 - k(x^2 + 2dx +d^2)}{(x+d)^2} = \frac{d^2 x^2 +2 d x k + k^2 - kx^2 - 2dkx -kd^2)}{(x+d)^2} = \frac{d^2 x^2 + k^2 - kx^2 -kd^2)}{(x+d)^2} = \frac{d^2 (x^2-k) + k^2 - kx^2)}{(x+d)^2} = \frac{d^2 (x^2-k) - k(x^2-k))}{(x+d)^2} = \frac{(d^2-k) (x^2-k)}{(x+d)^2} = (x^2-k)\frac{d^2-k}{(x+d)^2} $.

Since $d = \lfloor \sqrt k \rfloor$, $d < \sqrt k < d+1$ or $d^2 < k < d^2 + 2d +1$ or $-2d - 1 < d^2 - k < 0$, so $|d^2-k| < 2d+1$. Using this, $|y^2-k| < |x^2-k|\frac{2d+1}{(x+d)^2}| = |x^2-k|\frac{2d+1}{x^2+2dx+d^2} $, so $|y^2-k|< |x^2-k|$, and each iteration gets closer to the square root.

Since the starting iterate is $d$, all following iterates exceed $d$ so $|y^2-k| < |x^2-k|\frac{2d+1}{(d+d)^2}| < |x^2-k|\frac{2d+1}{4d^2}| < |x^2-k|\frac{1+1/(2d)}{2d}| \le 3|x^2-k|/4$ since $d \ge 1$.

This show that the iteration converges. However, this does not show that it converges quadratically like Newton's, only that it converges linearly.