The method starts with some initial guess, called $x_0$. It then
applies a formula to produce an improved guess. This is repeated until
the improved guess is accurate enough or it is clear the algorithm
fails to work.
For the Babylonian method, the next guess, $x_{i+1}$, is derived from the current guess, $x_i$. In mathematical notation, this is the updating step:
```math
x_{i+1} = \frac{1}{2}(x_i + \frac{k}{x_i})
```
We use this algorithm to approximate the square root of $2$, a value known to
the Babylonians.
Start with $x$, then form $x/2 + 1/x$, from this again form $x/2 + 1/x$, repeat.
We represent this step using a function
```julia
babylon(x) = x/2 + 1/x
```
Let's look starting with $x = 2$ as a rational number:
```julia; hold=true
x₁ = babylon(2//1)
x₁, x₁^2.0
```
Our estimate improved from something which squared to $4$ down to something which squares to $2.25.$ A big improvement, but there is still more to come. Had we done one more step:
```julia;
x₂ = (babylon ∘ babylon)(2//1)
x₂, x₂^2.0
```
We now see accuracy until the third decimal point.
```julia;
x₃ = (babylon ∘ babylon ∘ babylon)(2//1)
x₃, x₃^2.0
```
This is now accurate to the sixth decimal point. That is about as far
as we, or the Bablyonians, would want to go by hand. Using rational
numbers quickly grows out of hand. The next step shows the explosion.
```julia;
reduce((x,step) -> babylon(x), 1:4, init=2//1)
```
(In the above, we used `reduce` to repeat a function call ``4`` times, as an alternative to the composition operation. In this section we show a few styles to do this repetition before introducing a packaged function.)
However, with the advent of floating point numbers, the method stays quite manageable:
The tangent line and the function nearly agree near $2$. So much so,
that the intersection point of the tangent line with the $x$ axis
nearly hides the actual zero of $f(x)$ that is near $2.1$.
That is, it seems that the intersection of the tangent line and the
$x$ axis should be an improved approximation for the zero of the
function.
Let $x_0$ be $2$, and $x_1$ be the intersection point of the tangent line
at $(x_0, f(x_0))$ with the $x$ axis. Then by the definition of the
tangent line:
```math
f'(x_0) = \frac{\Delta y }{\Delta x} = \frac{f(x_0)}{x_0 - x_1}.
```
This can be solved for $x_1$ to give $x_1 = x_0 - f(x_0)/f'(x_0)$. In general, if we had $x_i$ and used the intersection point of the tangent line to produce $x_{i+1}$ we would have Newton's method:
```math
x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}.
```
Using automatic derivatives, as brought in with the `CalculusWithJulia` package, we can implement this algorithm.
The algorithm above starts at $2$ and then becomes:
```julia;
f(x) = x^3 - 2x - 5
x0 = 2.0
x1 = x0 - f(x0) / f'(x0)
```
We can see we are closer to a zero:
```julia;
f(x0), f(x1)
```
Trying again, we have
```julia;
x2 = x1 - f(x1)/ f'(x1)
x2, f(x2), f(x1)
```
And again:
```julia;
x3 = x2 - f(x2)/ f'(x2)
x3, f(x3), f(x2)
```
```julia;
x4 = x3 - f(x3)/ f'(x3)
x4, f(x4), f(x3)
```
We see now that $f(x_4)$ is within machine tolerance of $0$, so we
call $x_4$ an *approximate zero* of $f(x)$.
> **Newton's method:** Let $x_0$ be an initial guess for a zero of
> $f(x)$. Iteratively define $x_{i+1}$ in terms of the just
> generated $x_i$ by:
> ```math
> x_{i+1} = x_i - f(x_i) / f'(x_i).
> ```
> Then for
> reasonable functions and reasonable initial guesses, the sequence of
> points converges to a zero of $f$.
On the computer, we know that actual convergence will likely never
occur, but accuracy to a certain tolerance can often be achieved.
In the example above, we kept track of the previous values. This is
unnecessary if only the answer is sought. In that case, the update
step could use the same variable. Here we use `reduce`:
```julia;hold=true;
xₙ = reduce((x, step) -> x - f(x)/f'(x), 1:4, init=2)
xₙ, f(xₙ)
```
In practice, the algorithm is implemented not by repeating the update step a fixed number of times, rather by repeating the step until either we
converge or it is clear we won't converge. For good guesses and most
This interactive graphic (built using [JSXGraph](https://jsxgraph.uni-bayreuth.de/wp/index.html)) allows the adjustment of the point `x0`, initially at ``0.85``. Five iterations of Newton's method are illustrated. Different positions of `x0` clearly converge, others will not.
For the function $f(x) = \cos(x) - x$, we see that SymPy can not solve symbolically for a zero:
```julia;
@syms x::real
solve(cos(x) - x, x)
```
We can find a numeric solution, even though there is no closed-form answer. Here we try Newton's method:
```julia; hold=true
f(x) = cos(x) - x
x = .5
x = x - f(x)/f'(x) # 0.7552224171056364
x = x - f(x)/f'(x) # 0.7391416661498792
x = x - f(x)/f'(x) # 0.7390851339208068
x = x - f(x)/f'(x) # 0.7390851332151607
x = x - f(x)/f'(x)
x, f(x)
```
To machine tolerance the answer is a zero, even though the exact answer is irrational and all finite floating point values can be represented as rational numbers.
##### Example
Use Newton's method to find the *largest* real solution to ``e^x = x^6``.
A plot shows us roughly where the value lies:
```julia; hold=true
f(x) = exp(x)
g(x) = x^6
plot(f, 0, 25, label="f")
plot!(g, label="g")
```
Clearly by ``20`` the two paths diverge. We know exponentials eventually grow faster than powers, and this is seen in the graph.
To use Newton's method to find the intersection point. Stop when the increment ``f(x)/f'(x)`` is smaller than `1e-4`.
We need to turn the solution to an equation into a value where a function is ``0``. Just moving the terms to one side of the equals sign gives ``e^x - x^6 = 0``, or the ``x`` we seek is a solution to ``h(x)=0`` with ``h(x) = e^x - x^6``.
```julia; hold=true; term=true
h(x) = exp(x) - x^6
x = 20
for step in 1:10
delta = h(x)/h'(x)
x = x - delta
@show step, x, delta
end
```
So it takes ``8`` steps to get an increment that small and about `10` steps to get to full convergence.
##### Example division as multiplication
[Newton-Raphson Division](http://tinyurl.com/kjj9w92) is a means to divide by multiplying.
Why would you want to do that? Well, even for computers division is
harder (read slower) than multiplying. The trick is that $p/q$ is
simply $p \cdot (1/q)$, so finding a means to compute a reciprocal by
multiplying will reduce division to multiplication.
Well suppose we have $q$, we could try to use Newton's method to find
$1/q$, as it is a solution to $f(x) = x - 1/q$. The Newton update step
simplifies to:
```math
x - f(x) / f'(x) \quad\text{or}\quad x - (x - 1/q)/ 1 = 1/q
```
That doesn't really help, as Newton's method is just $x_{i+1} = 1/q$.
That is, it just jumps to the answer, the one we want to compute by some other means!
Trying again, we simplify the update step for a related function:
$f(x) = 1/x - q$ with $f'(x) = -1/x^2$ and then one step of the process is:
It can be shown that we have for any $q$ in $[1/2, 1]$ with initial guess $x_0 =
48/17 - 32/17\cdot q$ that Newton's method will converge to ``16`` digits in no more
than this many steps:
```math
\log_2(\frac{53 + 1}{\log_2(17)}).
```
```julia;
a = log2((53 + 1)/log2(17))
ceil(Integer, a)
```
That is ``4`` steps suffices.
For $q = 0.80$, to find $1/q$ using the above we have
```julia; hold=true
q = 0.80
x = (48/17) - (32/17)*q
x = -q*x*x + 2*x
x = -q*x*x + 2*x
x = -q*x*x + 2*x
x = -q*x*x + 2*x
```
This method has basically $18$ multiplication and addition operations
for one division, so it naively would seem slower, but timing this
shows the method is competitive with a regular division.
## Wrapping in a function
In the previous examples, we saw fast convergence, guaranteed converge in ``4`` steps, and an example where ``8`` steps were needed to get the requested level of approximation. Newton's method usually converges quickly, but may converge slowly, and may not converge at all. Automating the task to avoid repeatedly running the update step is
a task best done by the computer.
The `while` loop is a
good way to repeat commands until some condition is met. With this, we
present a simple function implementing Newton's method, we iterate
until the update step gets really small (the `atol`) or the
convergence takes more than ``50`` steps. (There are other, better choices that could be used to determine when the algorithm should stop, these are just easy to understand.)
```julia;
function nm(f, fp, x0)
atol = 1e-14
ctr = 0
delta = Inf
while (abs(delta) > atol) && (ctr < 50)
delta = f(x0)/fp(x0)
x0 = x0 - delta
ctr = ctr + 1
end
ctr < 50 ? x0 : NaN
end
```
##### Examples
- Find a zero of $\sin(x)$ starting at $x_0=3$:
```julia;
nm(sin, cos, 3)
```
This is an approximation for $\pi$, that historically found use, as the convergence is fast.
- Find a solution to $x^5 = 5^x$ near $2$:
Writing a function to handle this, we have:
```julia;
k(x) = x^5 - 5^x
```
We could find the derivative by hand, but use the automatic one instead:
```julia;
alpha = nm(k, k', 2)
alpha, f(alpha)
```
### Functions in the Roots package
Typing in the `nm` function might be okay once, but would be tedious
if it was needed each time. Besides, it isn't as robust to different inputs as possible. The `Roots` package provides a `Newton`
method for `find_zero`.
To use a different method with `find_zero`, the calling pattern is `find_zero(f, x, M)` where `f` represent the function(s), `x` the initial point(s), and `M` the method. Here we have:
```julia
find_zero((sin, cos), 3, Roots.Newton())
```
Or, if a derivative is not specified, one can be computed using automatic differentiation:
```julia; hold=true
f(x) = sin(x)
find_zero((f, f'), 2, Roots.Newton())
```
The argument `verbose=true` will force a print out of a message summarizing the convergence and showing each step.
Find the intersection point between $f(x) = \cos(x)$ and $g(x) = 5x$ near $0$.
We have Newton's method to solve for zeros of $f(x)$, i.e. when $f(x) =
0$. Here we want to solve for $x$ with $f(x) = g(x)$. To do so, we
make a new function $h(x) = f(x) - g(x)$, that is $0$ when $f(x)$
equals $g(x)$:
```julia; hold=true
f(x) = cos(x)
g(x) = 5x
h(x) = f(x) - g(x)
x0 = find_zero((h,h'), 0, Roots.Newton())
x0, h(x0), f(x0), g(x0)
```
----
We redo the above using a *parameter* for the ``5``, as there are some options on how it would be done. We let `f(x,p) = cos(x) - p*x`. Then we can use `Roots.Newton` by also defining a derivative:
```julia; hold=true
f(x,p) = cos(x) - p*x
fp(x,p) = -sin(x) - p
xn = find_zero((f,fp), pi/4, Roots.Newton(); p=5)
xn, f(xn, 5)
```
To use automatic differentiation is not straightforward, as we must hold the `p` fixed. For this, we introduce a closure that fixes `p` and differentiates in the `x` variable (called `u` below):
```julia; hold=true
f(x,p) = cos(x) - p*x
fp(x,p) = (u -> f(u,p))'(x)
xn = find_zero((f,fp), pi/4, Roots.Newton(); p=5)
```
##### Example: Finding $c$ in Rolle's Theorem
The function $r(x) = \sqrt{1 - \cos(x^2)^2}$ has a zero at $0$ and one at ``a`` near ``1.77``.
```julia;
r(x) = sqrt(1 - cos(x^2)^2)
plot(r, 0, 1.77)
```
As $f(x)$ is differentiable between $0$ and $a$, Rolle's theorem says
there will be value where the derivative is $0$. Find that value.
This value will be a zero of the derivative. A graph shows it should be near $1.2$, so we use that as a starting value to get the answer:
```julia;
find_zero((r',r''), 1.2, Roots.Newton())
```
## Convergence rates
Newton's method is famously known to have "quadratic convergence." What
does this mean? Let the error in the $i$th step be called $e_i = x_i -
\alpha$. Then Newton's method satisfies a bound of the type:
```math
\lvert e_{i+1} \rvert \leq M_i \cdot e_i^2.
```
If $M$ were just a constant and we suppose $e_0 = 10^{-1}$ then $e_1$
would be less than $M 10^{-2}$ and $e_2$ less than $M^2 10^{-4}$,
$e_3$ less than $M^3 10^{-8}$ and $e_4$ less than $M^4 10^{-16}$ which
for $M=1$ is basically the machine precision when values are near
``1``. That is for some problems, with a good initial guess it will
take around ``4`` or so steps to converge.
To identify ``M``, let ``\alpha`` be the zero of ``f`` to be approximated. Assume
* The function ``f`` has at continuous second derivative in a neighborhood of ``\alpha``.
* The value ``f'(\alpha)`` is *non-zero* in the neighborhood of ``\alpha``.
Then this linearization holds at each $x_i$ in the above neighborhood:
Suppose $f(x)$ is increasing and concave up. From the tangent line representation: $f(x) = f(c) + f'(c)\cdot(x-c) + f''(\xi)/2 \cdot(x-c)^2$, explain why it must be that the graph of $f(x)$ lies on or *above* the tangent line.
```julia; hold=true; echo=false
choices = [
L"As $f''(\xi)/2 \cdot(x-c)^2$ is non-negative, we must have $f(x) - (f(c) + f'(c)\cdot(x-c)) \geq 0$.",
L"As $f''(\xi) < 0$ it must be that $f(x) - (f(c) + f'(c)\cdot(x-c)) \geq 0$.",
L"This isn't true. The function $f(x) = x^3$ at $x=0$ provides a counterexample"
This question can be used to give a proof for the previous two questions, which can be answered by considering the graphs alone. Combined, they say that if a function is increasing and concave up and ``\alpha`` is a zero, then if ``x_0 < \alpha`` it will be ``x_1 > \alpha``, and for any ``x_i > \alpha``, ``\alpha <= x_{i+1} <= x_\alpha``, so the sequence in Newton's method is decreasing and bounded below; conditions for which it is guaranteed mathematically there will be convergence.
###### Question
Let $f(x) = x^2 - 3^x$. This has derivative $2x - 3^x \cdot
\log(3)$. Starting with $x_0=0$, what does Newton's method converge on?
```julia; hold=true; echo=false
f(x) = x^2 - 3^x;
fp(x) = 2x - 3^x*log(3);
val = Roots.newton(f, fp, 0);
numericq(val, 1e-14)
```
###### Question
Let $f(x) = \exp(x) - x^4$. There are 3 zeros for this function. Which one does Newton's method converge to when $x_0=2$?
```julia; hold=true; echo=false
f(x) = exp(x) - x^4;
fp(x) = exp(x) - 4x^3;
xstar= Roots.newton(f, fp, 2);
numericq(xstar, 1e-1)
```
###### Question
Let $f(x) = \exp(x) - x^4$. As mentioned, there are 3 zeros for this function. Which one does Newton's method converge to when $x_0=8$?
```julia; hold=true; echo=false
f(x) = exp(x) - x^4;
fp(x) = exp(x) - 4x^3;
xstar = Roots.newton(f, fp, 8);
numericq(xstar, 1e-1)
```
###### Question
Let $f(x) = \sin(x) - \cos(4\cdot x)$.
Starting at $\pi/8$, solve for the root returned by Newton's method
```julia; hold=true; echo=false
k1=4
f(x) = sin(x) - cos(k1*x);
fp(x) = cos(x) + k1*sin(k1*x);
val = Roots.newton(f, fp, pi/(2k1));
numericq(val)
```
###### Question
Using Newton's method find a root to $f(x) = \cos(x) - x^3$ starting at $x_0 = 1/2$.
```julia; hold=true; echo=false
f(x) = cos(x) - x^3
val = Roots.newton(f,f', 1/2)
numericq(val)
```
###### Question
Use Newton's method to find a root of $f(x) = x^5 + x -1$. Make a quick graph to find a reasonable starting point.
```julia; hold=true; echo=false
f(x) = x^5 + x - 1
val = Roots.newton(f,f', -1)
numericq(val)
```
###### Question
```julia; hold=true;echo=false
##Consider the following illustration of Newton's method:
caption = """
Illustration of Newton's method. Moving the point ``x_0`` shows different behaviours of the algorithm.
This function does not have a small first derivative; or a large second derivative; and the bump up can be made as close to the origin as desired, so the starting point can be very close to the zero. However, even though the conditions of the error term are satisfied, the error term does not apply, as ``f`` is not continuously differentiable.
###### Question
Let $f(x) = \sin(x) - x/4$. Starting at $x_0 = 2\pi$ Newton's method will converge to a value, but it will take many steps. Using the argument `verbose=true` for `find_zero`, how many steps does it take:
Can we find which point on its graph has the largest $y$ value?
This would be straightforward *if* we could write $y(x) = \dots$, for then we would simply find the critical points and investiate. But we can't so easily solve for $y$ interms of $x$. However, we can use Newton's method to do so:
(Using automatic derivatives works for values identified with `find_zero` *as long as* the initial point has its type the same as that of `x`.)
###### Question
In the last problem we used an *approximate* derivative (forward difference) in place of the derivative.
This can introduce an error due to the approximation. Would Newton's method still converge if the derivative in the algorithm were replaced with an approximate derivative? In general, this can often be done *but* the convergence can be *slower* and the sensitivity to a poor initial guess even greater.
Three common approximations are given by the
difference quotient for a fixed $h$: $f'(x_i) \approx (f(x_i+h)-f(x_i))/h$;
the secant line approximation: $f'(x_i) \approx (f(x_i) - f(x_{i-1})) / (x_i - x_{i-1})$; and the
