Finding Roots à la Newton Raphson: nr.py

A standard numerical problem asks for the solution to the equation

\begin{displaymath}
f(x) = 0
\end{displaymath}

for some (usually nonlinear) differentiable function $f(x)$. Among the many methods for solving this problem, the Newton-Raphson method is very efficient, but may require starting from a reasonably good guess for the solution.

Suppose that $x = p$ solves the problem. This means

\begin{displaymath}
f(p) = 0.
\end{displaymath}

We don't know $p$, but suppose we have guessed an $x$ that is close to $p$. We can expand $f(p)$ in a Taylor series about $x$ up to terms linear in $p$:

\begin{displaymath}
f(p) = f(x) + f^\prime(x)(p - x) + \ldots{}
\end{displaymath}

Using $f(p) = 0$, we can solve for $p$:

\begin{displaymath}
p = x - f(x)/f^\prime(x)
\end{displaymath}

The result is only approximate, since we dropped quadratic and higher terms in the Taylor series. If these neglected terms are small, $p$ is a better approximation to the root than $x$, but it isn't guaranteed to be the exact solution. So this formula can be used as the basis for an iterative algorithm: take the value of $p$ from the lhs, plug it in for $x$ on the rhs, and repeat. As we get closer to the root, the higher order terms in the Taylor expansion become even less important and the method converges very rapidly.

How do we know when the method has converged? That is, when do we tell our program to stop? One stopping criterion is to monitor the improvement $\vert p-x\vert$ and stop when it is less than a specified tolerance. Another is to monitor the value of $\vert f(p)\vert$ and stop when it is small enough. Which one is best depends on our application, i.e. for what purpose are we solving this equation in the first place? Do we require knowing $p$ with high accuracy or does it suffice to get a suitably small $f(p)$? Note that even if the change $\vert p-x\vert$ is smaller than some tolerance, there are no assurances that we have the root to that accuracy. However, for most situations, that method is sufficient.