Newton's method

Solving equations with calculus

It's good to understand that there are many equations (even very simple equations) that cannot be easily solved in closed form. As two seemingly simple examples, let's take a look at

  • $\cos(x)=x$, and
  • $x^3-3x-1=0$.

It's easy to see that these equations have solutions by simply looking at some graphs. We'll use a couple of Python packages to do this. Note that you can try the code out yourself at http://try.jupyter.org/!

In [1]:
%matplotlib inline
from sympy import *
from sympy.plotting import plot
x = symbols('x')
import warnings
warnings.filterwarnings('ignore')

Here's the graph of the first equation.

In [2]:
plot(cos(x), x, (x,-0.2,1.2))
Out[2]:
<sympy.plotting.plot.Plot at 0x106584748>

That point of intersection shows that there is a solution. Let's ask the computer for a symbolic representation of the solution.

In [3]:
solve(cos(x)-x)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-3-7a8e53eb2234> in <module>()
----> 1 solve(cos(x)-x)

/Users/mcmcclur/anaconda3/lib/python3.5/site-packages/sympy/solvers/solvers.py in solve(f, *symbols, **flags)
    907     ###########################################################################
    908     if bare_f:
--> 909         solution = _solve(f[0], *symbols, **flags)
    910     else:
    911         solution = _solve_system(f, symbols, **flags)

/Users/mcmcclur/anaconda3/lib/python3.5/site-packages/sympy/solvers/solvers.py in _solve(f, *symbols, **flags)
   1412     if result is False:
   1413         raise NotImplementedError(msg +
-> 1414         "\nNo algorithms are implemented to solve equation %s" % f)
   1415 
   1416     if flags.get('simplify', True):

NotImplementedError: multiple generators [x, cos(x)]
No algorithms are implemented to solve equation -x + cos(x)

Well, that's discouraging! How about a simple cubic polynomial?

In [4]:
plot(x**3-3*x-1, (x,-2,2))
Out[4]:
<sympy.plotting.plot.Plot at 0x1095004e0>

Looks like there are three roots. Let's find them!

In [5]:
init_printing(pretty_print=True)
solve(x**3-3*x-1)
Out[5]:
$$\left [ \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}} + \frac{1}{\left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}}}, \quad \frac{1}{\left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}}} + \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}}, \quad \frac{1}{\sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}}} + \sqrt[3]{\frac{1}{2} + \frac{\sqrt{3} i}{2}}\right ]$$

Crazy!

Newton's method

Newton'd method tells us that we can approximate the root of a function by iterating the related Newton's method iteration function: $$n(x) = x - \frac{f(x)}{f'(x)}.$$ We should start the iteration near the root we want to find. For example, the previous graph shows that $f(x)=x^3-3x-1$ has a root slightly less than $2$. If we set $x_1=2$, $x_2=n(x_1)$, and (generally) $x_i=f(x_{i-1})$ we should generate a sequence that converges to the root. Let's try it!

In [6]:
init_printing(pretty_print=False)
def f(x): return x**3 - 3*x - 1
n = lambdify(x, x-f(x)/diff(f(x),x))
xi = 2.0
cnt = 0
while abs(f(xi)) > 10.0**(-10):
    cnt = cnt+1
    xi = n(xi)
(xi,cnt)
Out[6]:
(1.8793852415718169, 4)

This tells us that in 4 iterations, Newton's method found an estimate to the solution $x_i$ so that $|f(x_i)|<1/10^10$ - i.e., a residual accurate to within 10 digits. Not bad!

We can apply this to $\cos(x)=x$ by first writing the equation as $\cos(x)-x=0$ and finding a root of $f(x)=\cos(x)-x$.

In [8]:
def f(x): return cos(x)-x
n = lambdify(x, x-f(x)/diff(f(x),x))
xi = 2.0
cnt = 0
while abs(f(xi)) > 10.0**(-10):
    cnt = cnt+1
    xi = n(xi)
(xi,cnt)
Out[8]:
(0.7390851332198146, 3)