@def title="Root-finding without a calculator" @def tags = ["maths", "numerical methods"] @def mintoclevel=1
In this article, I will introduce techniques for approximating square and cube roots. The methods used to derive these techniques are very general and can be used to find the roots of any continuous real-valued function.
\tableofcontents
A popular technique for approximating the square root of a number is to use the tangent line approximation. Say we wish to find
If we have a function
[y = mx + c]
as it is a straight line. At the point
\begin{eqnarray} y &=& f(x) \\ y' &=& f'(x) \\ \end{eqnarray}
or in other words:
\begin{eqnarray} mx_0 + c &=& f(x_0) \\ m &=& f'(x_0) \end{eqnarray}
substituting the second of these equations into the first yields:
\begin{eqnarray} f'(x_0)x_0 + c &=& f(x_0) \\ c &=& f(x_0) - f'(x_0) x_0 \\ \therefore y &=& f'(x_0) x + f(x_0) - f'(x_0) x_0 \\ &=& f'(x_0)(x-x_0) + f(x_0). \end{eqnarray}
Which is the equation of the tangent line. If we let
\begin{eqnarray} f'(x) &=& \dfrac{1}{2\sqrt{x}} \\ y &=& \dfrac{1}{2\sqrt{x_0}} (x-x_0) + \sqrt{x_0} \\ &=& \sqrt{x_0} + \dfrac{x-x_0}{2\sqrt{x_0}}. \end{eqnarray}
If we let
\begin{eqnarray} y &=& \sqrt{b} + \dfrac{a-b}{2\sqrt{b}}. \end{eqnarray}
Which is our approximation for
While its approximation is often satisfactory, it can be substantially off if
In this case, you can use Newton's method to refine the approximation. Newton's method is a technique in which we essentially use the tangent line to approximate the zeros of a function. If one of that function's roots happens to be the square root you're searching for, applying Newton's method will have the result of giving you ever improving approximations to the square root. One function that is easy to exactly compute that's root is
\begin{eqnarray} x_{n+1} &=& x_n - \dfrac{f(x_n)}{f'(x_n)} \\ &=& x_n - \dfrac{x_n^2 - a}{2x_n} \\ &=& x_n - \dfrac{x_n}{2} + \dfrac{a}{2x_n} \\ &=& \dfrac{x_n}{2} + \dfrac{a}{2x_n}. \end{eqnarray}
Where
A classic example where the tangent line approximation is pretty poor is when we are trying to find
\begin{eqnarray} \sqrt{2} &\approx& \sqrt{1} + \dfrac{2-1}{2\sqrt{1}} \\ &=& 1 + \dfrac{1}{2} \\ &=& \dfrac{3}{2} \\ &=& 1.5. \end{eqnarray}
When of course
\begin{eqnarray} x_1 &=& \dfrac{x_0}{2} + \dfrac{2}{2x_0} \\ &=& \dfrac{\dfrac{3}{2}}{2} + \dfrac{2}{2\times \dfrac{3}{2}} \\ &=& \dfrac{3}{4} + \dfrac{2}{3} \\ &=& \dfrac{3\times 3}{4\times 3} + \dfrac{2\times 4}{3\times 4} \\ &=& \dfrac{9}{12} + \dfrac{8}{12} \\ &=& \dfrac{17}{12} \\ &\approx& 1.417. \end{eqnarray}
Now this is clearly not accurate to five decimal places, so another iteration of Newton's method is required.
\begin{eqnarray} x_2 &=& \dfrac{x_1}{2} + \dfrac{2}{2x_1} \\ &=& \dfrac{\dfrac{17}{12}}{2} + \dfrac{1}{\dfrac{17}{12}} \\ &=& \dfrac{17}{24} + \dfrac{12}{17} \\ &=& \dfrac{17 \times 17}{24 \times 17} + \dfrac{12 \times 24}{17 \times 24} \\ &=& \dfrac{289}{408} + \dfrac{288}{408} \\ &=& \dfrac{577}{408} \ &\approx& 1.4142157. \end{eqnarray}
This is accurate to four decimal places, the fifth is not as we would round it up to two, not one. Running Newton's method again yields:
\begin{eqnarray} x_3 &=& \dfrac{x_2}{2} + \dfrac{2}{2x_2} \\ &=& \dfrac{\dfrac{577}{408}}{2} + \dfrac{1}{\dfrac{577}{408}} \\ &=& \dfrac{577}{816} + \dfrac{408}{577} \\ &=& \dfrac{577 \times 577}{816 \times 577} + \dfrac{408 \times 816}{577 \times 816} \\ &=& \dfrac{332929}{470832} + \dfrac{332928}{470832} \\ &=& \dfrac{665857}{470832} \ &\approx& 1.41421356237469. \end{eqnarray}
Which is an accurate approximation of
The tangent line approximation yields:
\begin{eqnarray} \sqrt{76} &\approx& \sqrt{81} - \dfrac{5}{2\sqrt{81}} \\ &=& 9 - \dfrac{5}{2\times 9} \\ &=& 9 - \dfrac{5}{18} \\ &=& \dfrac{9\times 18}{18} - \dfrac{5}{18} \\ &=& \dfrac{162}{18} - \dfrac{5}{18} \\ &=& \dfrac{157}{18} \\ &\approx& 8.722. \end{eqnarray}
The actual value of
\begin{eqnarray} x_1 &=& \dfrac{x_0}{2} + \dfrac{76}{2x_0} \\ &=& \dfrac{157}{18 \times 2} + \dfrac{76}{2 \times \dfrac{157}{18}} \\ &=& \dfrac{157}{36} + \dfrac{38}{\dfrac{157}{18}} \\ &=& \dfrac{157}{36} + \dfrac{38 \times 18}{157} \\ &=& \dfrac{157 \times 157}{36 \times 157} + \dfrac{684 \times 36}{157 \times 36} \\ &=& \dfrac{24649}{5652} + \dfrac{24624}{5652} \\ &=& \dfrac{49273}{5652} \\ &\approx& 8.717799. \end{eqnarray}
Which is accurate to five decimal places.
The tangent line approximation yields:
\begin{eqnarray} \sqrt{5} &\approx& \sqrt{4} + \dfrac{1}{2\sqrt{4}} \\ &=& 2 + \dfrac{1}{4} \\ &=& \dfrac{9}{4} \\ &=& 2.25. \end{eqnarray}
The actual value of
\begin{eqnarray} x_1 &=& \dfrac{x_0}{2} + \dfrac{5}{2x_0} \\ &=& \dfrac{\dfrac{9}{4}}{2} + \dfrac{5}{2 \times \dfrac{9}{4}} \\ &=& \dfrac{9}{8} + \dfrac{10}{9} \\ &=& \dfrac{9 \times 9 + 8 \times 10}{9 \times 8} \\ &=& \dfrac{161}{72} \\ &=& 2.236\overline{11}. \end{eqnarray}
So this is only accurate to three decimal places. The next iteration of Newton's method yields:
\begin{eqnarray} x_2 &=& \dfrac{x_1}{2} + \dfrac{5}{2x_1} \\ &=& \dfrac{\dfrac{161}{72}}{2} + \dfrac{5}{2 \times \dfrac{161}{72}} \\ &=& \dfrac{161}{144} + \dfrac{5 \times 36}{161} \\ &=& \dfrac{161 \times 161 + 5 \times 36 \times 144}{161 \times 144} \\ &=& \dfrac{25921 + 25920}{23184} \\ &=& \dfrac{51841}{23184} \\ &\approx& 2.23606797791. \end{eqnarray}
Which is accurate to nine decimal places, or eight if we account for the effect of rounding.
Similarly the cube root of a number can be approximated by a tangent line approximation and refined using Newton's method. If we let
\begin{eqnarray} y &=& f'(x_0)(x-x_0) + f(x_0) \\ &=& \dfrac{1}{3(\sqrt[3]{x_0})^2} (x-x_0) + \sqrt[3]{x_0} \\ &=& \sqrt[3]{x_0} + \dfrac{x-x_0}{3(\sqrt[3]{x_0})^2}. \end{eqnarray}
Letting
\begin{eqnarray} a_c &\approx& \sqrt[3]{b} + \dfrac{a-b}{3(\sqrt[3]{b})^2}. \end{eqnarray}
And to refine this approximation we can use Newton's method with
\begin{eqnarray} x_{n+1} &=& x_n - \dfrac{f(x_n)}{f'(x_n)} \\ &=& x_n - \dfrac{x_n^3-a}{3x_n^2} \\ &=& x_n - \dfrac{x_n}{3} + \dfrac{a}{3x_n^2} \\ &=& \dfrac{2x_n}{3} + \dfrac{a}{3x_n^2}. \end{eqnarray}
Here we use
\begin{eqnarray} \sqrt[3]{63} &\approx& \sqrt[3]{64} - \dfrac{1}{3(\sqrt[3]{64})^2} \\ &=& 4 - \dfrac{1}{3\times 16} \\ &=& 4 - \dfrac{1}{48} \\ &=& \dfrac{191}{48} \\ &\approx& 3.9791\overline{6}. \end{eqnarray}
The correct value of
\begin{eqnarray} x_1 &=& \dfrac{2x_0}{3} + \dfrac{63}{3x_0^2} \\ &=& \dfrac{2x_0}{3} + \dfrac{21}{x_0^2} \\ &=& \dfrac{2 \times \dfrac{191}{48}}{3} + \dfrac{21}{\left(\dfrac{191}{48}\right)^2} \\ &=& \dfrac{191}{72} + \dfrac{21 \times 48^2}{191^2} \\ &=& \dfrac{191}{72} + \dfrac{21 \times 2304}{36481} \\ &=& \dfrac{191 \times 36481 + 21 \times 2304 \times 72}{72 \times 36481} \\ &=& \dfrac{6967871 + 3483648}{2626632} \\ &=& \dfrac{10451519}{2626632} \\ &\approx& 3.97905721. \end{eqnarray}
Which is accurate to seven decimal places, or eight if we round off the remaining digits from there.
Newton's method can be used for finding the roots of any continuous real-valued function, or even system of continuous real-valued functions, although it has some shortcomings. One is that it requires an initial guess as to the root, and that its results can depend heavily on this initial guess. Additionally, it can sometimes fail to converge, such as when the root occurs at a stationary point for the function (maxima or minima), as there the derivative is zero. For example, it fails for the polynomial
Newton's method belongs to a class of methods known as the Householder's methods. Newton's method is the first order Householder's method, Halley's method is the second order Householder's method and it is possible to derive the $d$th order Householder's method using the formula:
\begin{eqnarray} x_{n+1} = x_n + d\dfrac{(1/f)^{(d-1)}(x_n)}{(1/f)^{(d)}(x_n)} \end{eqnarray}
where the bracketed superscripts denotes derivatives of the superscript's order. For example
Getting the initial guesses for the roots, especially of polynomials that have multiple real roots, can be a challenge. One technique is to use the bisection method. In this technique one takes an interval over which the sign of the function changes from positive to negative or vice versa (which implies that at least one root must lie within this interval) and continuously subdivides this interval and re-evaluates the function at its endpoints until we find the root. This method is very slow, so usually one would only apply it a few times, and then once our interval is acceptably small we would apply Newton's method to get the root. One can also evaluate the function at a long list of points within an interval and find where in the interval the sign changes, as that is where a root will be. If the interval is large enough and the function has multiple real roots there may be multiple sign change points and hence multiple roots we can converge to using Newton's method. Like Newton's method, however, it does not converge when the root is a minima or maxima.
These methods can be done by hand, but for many problems they can be the quite tedious; as such I wrote this [Julia](https://en.wikipedia.org/wiki/Julia_(programming language)) script that implements the bisection method and Newton's method to find the roots of a given function:
#!/usr/bin/env julia
"""
bisection(f::Function, N::Integer, a::Number, b::Number)
Uses the bisection method with N subdivisions of the specified domain [a, b]
to find the roots of the function f within said domain. Used to find the
initial guess that Newton's method then uses to get a more precise estimate of
the roots.
"""
function bisection(f, N, a, b)
x = LinRange(a, b, N+1)
fval = f.(x)
xv = x[2:end]
initGuess = xv[diff(sign.(fval)).!= 0]
return initGuess
end
"""
newtons(f::Function, h::Float, tol::Float, itMax::Integer,
initGuess::Vector{Float})
Uses Newton's method (and if needed, up to third order Householder's methods)
to refine initGuess, our initial guess of the root(s) of f, until either itMax
iterations has been performed or the relative magnitude of the update Newton's
provides to the estimate is below tol. h is the step size used to approximate
the derivative.
"""
function newtons(f, h, tol, itMax, initGuess)
# Derivatives
function fd(x, h)
return (f(x+h)-f(x-h))/(2*h)
end
function f2d(x, h)
return (f(x+h)-2*f(x)+f(x-h))/(h^2)
end
function f3d(x, h)
return (f(x+2*h)-2*f(x+h)+f(x-h)-f(x-2*h))/(2*h^3)
end
# Check to see if (relative) difference in update is below tol
function tolCheck(diff, sol, tol)
if ! (sol ≈ 0)
check = abs(diff/sol) > tol
else
check = abs(diff) > tol
end
return check
end
# Initialize sol
sol = initGuess
count = Int64.(zeros(size(initGuess)))
for j=1:length(initGuess)
diff = 1
while (tolCheck(diff, sol[j], tol) && count[j] < itMax)
fn = f(sol[j])
der = fd(sol[j], h)
der2 = f2d(sol[j], h)
der3 = f3d(sol[j], h)
if (der ≈ 0) && (der2 ≈ 0)
# Third order Householder's method
diff = (6*fn*der^2 - 3*fn^2*der2)/(6*der^3-6*fn*der*der2 + fn^2*der3)
elseif (der ≈ 0)
# Second order, Halley's method
diff = 2*fn*der/(-fn*der2 + 2*der^2)
else
# First order, Newton's method
diff = fn/der
end
sol[j] -= diff
count[j] += 1
end
if (count[j] == itMax)
print("Maximum iterations exceeded and the amount by which ")
println("Newton's last updated the solution was: ", diff)
end
end
# Delete NaN entries from sol and corresponding count entries
# While loop is used because sol's length will change over this loop
j = 1;
while j < length(sol)
if isnan(sol[j])
deleteat!(sol, j)
deleteat!(count, j)
end
j += 1;
end
return sol, count
end
"""
findRoot(f::Function, h::Float, tol::Float, itMax::Integer, a::Number, b::Number, N::Integer)
Uses the bisection method to get an initial guess of the root(s) of f on the
domain [a, b] with N subdivisions, then uses Newton's method with a maximum of
itMax iterations and a relative error tolerance of tol. h is the step size
used to approximate the derivative.
"""
function findRoot(f, h, tol, itMax, a, b, N)
initGuess = bisection(f, N, a, b)
sol, count = newtons(f, h, tol, itMax, initGuess);
return sol, count
end
"""
printSoln(sol::Vector{Float}, count::Vector{Integer})
Print the solution and the number of iterations used.
"""
function printSoln(sol, count)
if (length(sol) == 1)
println("Root = ", sol[1], ", count = ", count[1])
else
for i=1:length(sol)
println("The $(i)th root = ", sol[i], ", count = ", count[i])
end
end
end
# This is where you specify the function you want to
# find the root of
function f(x)
return x^4 + x^3 - 10x^2 - 4x + 16
end
# Initialize required variables
h = 1e-10
tol = 1e-15
itMax = 1000
a = -100
b = 100
N = 100000
sol, count = findRoot(f, h, tol, itMax, a, b, N)
printSoln(sol, count)
.