# Newton's Method

Kuan Xu, 24th October 2012

(Chebfun example roots/NewtonRaphson.m)
[Tags: #root, #Newton]

Newton's method, as the most fundamental root-finding algorithm, usually appears no later than Chapter 2 in most numerical analysis textbooks. It uses the first two terms of the Taylor series of a function $f(x)$ in the vicinity of a suspected root to find successively better approximations to the root, using the formula $$x^{(k+1)} = x^{(k)} - {f(x^{(k)})\over f'(x^{(k)})}.$$

Let's consider $f(x) = x^3-3x^2+2$, which has several roots, as we see in the next plot.

LW = 'linewidth'; lw = 2;
MS = 'MarkerSize'; ms = 18;
dom = [-3 3];
f = chebfun('x.^3-3*x.^2+2', dom);
plot(f, LW, lw), hold on
plot(dom, [0 0], 'k'), hold off

ans =
v: 0
h: 3


Here are the roots.

roots(f)

ans =
-0.732050807568881
1.000000000000003
2.732050807568879


If we try to locate the leftmost root by Newton's method, we need to pick an initial guess, for example $-3$.

fprime = diff(f);
d = norm(f,inf);
tol = 1e-8;
xold = -2;
x = [];
i = 0;

plot(f, LW, lw), xlim([-2.5 0]), hold on
plot(dom, [0 0], 'k')
while(d > tol)
x = [x xold];
xnew = xold - f(xold)/fprime(xold);
d = abs(xnew - xold);

plot(xold, f(xold), 'ok')
strx = ['x_{' num2str(i) '}'];
text(xold - 0.05, 1.2, strx,'fontsize',12)
plot(xold, 0, '.k', MS, ms)
plot([xold xold], [0 f(xold)], '--k', LW, lw)
plot([xold xnew], [f(xold) 0], '-.k', LW, lw)

xold = xnew;
i = i+1;
end
hold off
root1 = xnew

root1 =
-0.732050807568878


In the above plot, the solid black dots are the successive approximations of the root, while the circles are their projections on the curve, from which the Newton's method locates the next approximation along the tangent (black dash-dot lines). The following table tells us with no surprise that Newton's method is quadratically convergent.

n = size(x,2);
res = abs(x - xnew);
LogRes = log(res);

disp('iterations     Logarithm of the step size')
for i = 1:n
fprintf('%5d  %25.8f\n', i, LogRes(i))
end

iterations     Logarithm of the step size
1                 0.23740079
2                -0.65787813
3                -1.98646163
4                -4.28606060
5                -8.73432460
6               -17.61270707
7               -35.12736266


You may expect that the same order of convergence to appear when we approximate the middle root in this example.

d = norm(f,inf);
xold = 0.5;
x = [];
while(d > tol)
x = [x xold];
xnew = xold - f(xold)/fprime(xold);
d = abs(xnew - xold);
xold = xnew;
end
hold off
root2 = xnew

n = size(x,2);
res = abs(x - xnew);
LogRes = log(res);

disp('iterations     Logarithm of the step size')
for i = 1:n
fprintf('%5d  %25.8f\n', i, LogRes(i))
end

root2 =
0.999999999999999
iterations     Logarithm of the step size
1                -0.69314718
2                -2.19722458
3                -6.98471632
4               -21.35961332


But now, this time we are evidently achieving cubic convergence! Is there something wrong? No, this makes perfect sense, if you notice that $f^{''}(1) = 0$ in this example.