Skip to article frontmatterSkip to article content

3.4Newton method

#%config InlineBackend.figure_format = 'svg'
from pylab import *
import sympy

3.4.1The method

The Newton-Raphson method is an iterative method to find a root of a function. It starts with an initial guess for the root and tries to improve it. Given an initial guess x0x_0 for the root, we have

f(x0)≠0f(x_0) \ne 0
(1)

Let us try to find a correction Δx0\Delta x_0 which brings it closer to the root

f(x0+Δx0)=0f(x_0 + \Delta x_0) = 0
(2)

Using Taylor formula, and truncating it at the linear term

f(x0)+Δx0f′(x0)≈0⟹Δx0≈−f(x0)f′(x0)f(x_0) + \Delta x_0 f'(x_0) \approx 0 \qquad \Longrightarrow \qquad \Delta x_0 \approx - \frac{f(x_0)}{f'(x_0)}
(3)

We obtain a new estimate for the root

x1=x0+Δx0=x0−f(x0)f′(x0)x_1 = x_0 + \Delta x_0 = x_0 - \frac{f(x_0)}{f'(x_0)}
(4)

Geometrically, we are approximating the function by the tangent line at x0x_0 which is given by

y=y(x)=f(x0)+(x−x0)f′(x0)y = y(x) = f(x_0) + (x-x_0) f'(x_0)
(5)

and finding the zero of the tangent line,

y(x1)=0⟹x1=x0−f(x0)f′(x0)y(x_1) = 0 \limplies x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}
(6)

see figure. We now repeat the same process for x1,x2,…x_1, x_2, \ldots until some convergence is achieved.

This algorithm requires that ff is differentiable and f′f' is not zero in a neighbourhood of the root. The iterations may not converge in a finite number of steps, and we need to choose some stopping criteria. We can use a check on the function value and another on the root; we stop when

∣f(xn+1)∣<δand/or∣xn+1−xn∣<∣xn+1∣ϵ|f(x_{n+1})| < \delta \qquad \textrm{and/or} \qquad |x_{n+1} - x_n| < |x_{n+1}| \epsilon
(8)

returns true. The motivation for the second criterion is explained below and it is a test for the error in the location of the root.

The algorithm for Newton-Raphson method is given below.

Let us apply this method on some functions.

3.4.2Convergence analysis

Let rr be a root of ff. From Taylor formula with remainder term

0=f(r)=f(xn)+(r−xn)f′(xn)+12(r−xn)2f′′(ξn),ξn between r and xn0 = f(r) = f(x_n) + (r - x_n) f'(x_n) + \half (r - x_n)^2 f''(\xi_n), \qquad \xi_n \mbox{ between $r$ and $x_n$}
(12)

Solve for rr from the linear term

r=xn−f(xn)f′(xn)⏟xn+1−12(r−xn)2f′′(ξn)f′(xn)r = \underbrace{x_n - \frac{f(x_n)}{f'(x_n)}}_{x_{n+1}} - \half (r - x_n)^2 \frac{f''(\xi_n)}{f'(x_n)}
(13)

The error in the root is

r−xn+1=−(r−xn)2f′′(ξn)2f′(xn)r - x_{n+1} = - (r - x_n)^2 \frac{f''(\xi_n)}{2 f'(x_n)}
(14)

Atleast in a neighborhood of the root, we will assume that

∣f′′(ξn)2f′(xn)∣≤M\left| \frac{f''(\xi_n)}{2 f'(x_n)} \right| \le M
(15)

Then the error en=∣r−xn∣e_n = |r - x_n| satisfies

en+1≤Men2e_{n+1} \le M e_n^2
(16)

To show convergence xn→rx_n \to r, we have to show that en→0e_n \to 0.

3.4.3Error estimate

From mean value theorem

f(xn)=f(xn)−f(r)=f′(ξn)(xn−r)f(x_n) = f(x_n) - f(r) = f'(\xi_n) (x_n - r)
(29)

where ξn\xi_n is between xnx_n and rr, so that

r−xn=−f(xn)f′(ξn)r - x_n = - \frac{f(x_n)}{f'(\xi_n)}
(30)

If f′f' is not changing rapidly between xnx_n and rr then

f′(ξn)≈f′(xn)f'(\xi_n) \approx f'(x_n)
(31)

We get the error estimate

r−xn≈−f(xn)f′(xn)=xn+1−xnr - x_n \approx - \frac{f(x_n)}{f'(x_n)} = x_{n+1} - x_n
(32)

and the relative error estimate

r−xnr≈xn+1−xnxn+1\frac{r - x_n}{r} \approx \frac{x_{n+1} - x_n}{x_{n+1}}
(33)

Hence we can use the test

∣xn+1−xn∣<ϵ∣xn+1∣,0<ϵ≪1|x_{n+1} - x_n| < \epsilon |x_{n+1}|, \qquad 0 < \epsilon \ll 1
(34)

to decide on convergence. This is used in the code examples to stop the iterations.

3.4.4Implication of quadratic convergence

What does quadratic convergence imply in practical terms ? For Newton method

ek+1≤Mek2e_{k+1} \le M e_k^2
(39)

Suppose M≈1.0M \approx 1.0 and e0=10−1e_0 = 10^{-1}. Then

e1≈10−2e2≈10−4e3≈10−8e4≈10−16(limit of double precision)\begin{align*} & e_1 \approx 10^{-2} \\ & e_2 \approx 10^{-4} \\ & e_3 \approx 10^{-8} \\ & e_4 \approx 10^{-16} \qquad \textrm{(limit of double precision)} \end{align*}
(40)

If the root r=O(1)r = O(1), then

x1 accurate to 2 digitsx2 accurate to 4 digitsx3 accurate to 8 digitsx4 accurate to 16 digits\begin{align*} & x_1 \textrm{ accurate to 2 digits} \\ & x_2 \textrm{ accurate to 4 digits} \\ & x_3 \textrm{ accurate to 8 digits} \\ & x_4 \textrm{ accurate to 16 digits} \end{align*}
(41)

Newton’s method doubles the number of accurate figures in each iteration. Of course this happens only after the method enters the quadratic convergence regime. In the first few iterations, there may be no reduction of error or less than quadratic reduction.

Observe the number of correct decimal places in the square root example, which roughly doubles in each iteration.

3.4.5Inflection point

If we have a function with an inflection point, see figure, then the Newton iterates may cycle indefinitely without any convergence. The Newton method is

Δxk=−f(xk)f′(xk),xk+1=xk+Δxk\Delta x_k = - \frac{f(x_k)}{f'(x_k)}, \qquad x_{k+1} = x_k + \Delta x_k
(46)

We can use ∣f(x)∣|f(x)| as a metric to measure convergence towards the root: ∣f(xk)∣|f(x_k)| must decrease monotonically towards zero. We will accept the step Δxk\Delta x_k if ∣f(xk+1)∣<∣f(xk)∣| f(x_{k+1})| < |f(x_k)|, otherwise halve the step size Δxk←12Δxk\Delta x_k \leftarrow \half \Delta x_k. We can halve the step size repeatedly until reduction in ∣f(x)∣|f(x)| is obtained; such a step size always exists. This can be summarised in this algorithm.

  • Δx=−f(x)/f′(x)\Delta x = -f(x)/f'(x)
  • While ∣f(x+Δx)∣>∣f(x)∣|f(x+\Delta x)| > |f(x)|
    • Δx=Δx/2\Delta x = \Delta x/2

Moreover, to prevent wild oscillations of the iterates, the step size can also be restricted, e.g., dont allow step size to more than double in successive steps

∣Δxk∣≤2∣Δxk−1∣|\Delta x_k| \le 2 |\Delta x_{k-1}|
(47)

3.4.6Implicit functions

Suppose we have a function y=g(x):R→Ry=g(x) : \re \to \re, which is defined implicitly by

G:R×R→R,G(x,y)=0G: \re \times \re \to \re, \qquad G(x,y) = 0
(51)

We may not be able find an explicit expression g(x)g(x), instead we find the value of yy for some discrete set of xx values. This leads to a root finding problem for each xx and we can apply Newton method.

In this context, recall the implicit function theorem.

For a given value of xx and make an initial guess y0y_0. This may not satisfy G(x,y0)=0G(x,y_0) = 0 so we want to find a correction y0+Δy0y_0 + \Delta y_0 such that G(x,y0+Δy0)=0G(x,y_0+\Delta y_0) = 0. By Taylor expansion wrt second argument of GG

G(x,y0)+Gy(x,y0)Δy0=0⟹Δy0=−G(x,y0)Gy(x,y0)G(x,y_0) + G_y(x,y_0) \Delta y_0 = 0 \quad\Longrightarrow\quad \Delta y_0 = - \frac{G(x,y_0)}{G_y(x,y_0)}
(52)

and hence we get an improved estimate

y1=y0+Δy0=y0−G(x,y0)Gy(x,y0)y_1 = y_0 + \Delta y_0 = y_0 - \frac{G(x,y_0)}{G_y(x,y_0)}
(53)

This can be used to generate an iterative scheme

yk+1=yk−G(x,yk)Gy(x,yk),k=0,1,2,…y_{k+1} = y_k - \frac{G(x,y_k)}{G_y(x,y_k)}, \qquad k=0,1,2,\ldots
(54)

If we have to solve this problem for several values of xx, then we can use the previous solution as an initial guess for the root, which can speed up the convergence of Newton method.

3.4.7System of equations

Consider a function f:RN→RNf : \re^N \to \re^N such that for x∈RNx \in \re^N

f(x)=[f1(x),…,fN(x)]⊤∈RNf(x) = [f_1(x), \ldots, f_N(x)]^\top \in \re^N
(55)

We want to find an r∈RNr \in \re^N such that f(r)=0f(r) = 0. This is a system of NN non-linear equations for the NN unknowns. We start with an initial guess x0∈RNx^0 \in \re^N and in general f(x0)≠0f(x^0) \ne 0. We will try to find a correction Δx0\Delta x^0 to x0x^0 such that

f(x0+Δx0)=0f(x^0 + \Delta x^0) = 0
(56)

Approximating by the first two terms of the Taylor expansion

f(x0)+f′(x0)Δx0=0f(x^0) + f'(x^0) \Delta x^0 = 0
(57)

we can find the correction by solving the matrix problem

f′(x0)Δx0=−f(x0)f'(x^0) \Delta x^0 = - f(x^0)
(58)

where f′∈RN×Nf' \in \re^{N \times N} is given by

[f′]ij=∂fi∂xj[f']_{ij} = \df{f_i}{x_j}
(59)

Once we solve for Δx0=−[f′(x0)]−1f(x0)\Delta x^0 = -[f'(x^0)]^{-1} f(x^0), we get a new estimate of the root

x1=x0+Δx0=x0−[f′(x0)]−1f(x0)x^1 = x^0 + \Delta x^0 = x^0 - [f'(x^0)]^{-1} f(x^0)
(60)

The above steps are applied iteratively. Choose an initial guess x0x^0 and set k=0k=0.

  1. Solve the matrix problem

    f′(xk)Δxk=−f(xk)f'(x^k) \Delta x^k = - f(x^k)
    (61)

  2. Update the root

    xk+1=xk+Δxkx^{k+1} = x^k + \Delta x^k
    (62)

  3. k=k+1k=k+1

  4. If not converged, go to Step 1.

Numerical Analysis
Bisection method
Numerical Analysis
Secant Method