Back to M331: Matlab Codes, Notes and Links

Newton's Method in Matlab

Suppose we want to find the first positive root of the function

g(x)=sin(x)+x cos(x).

Since

g'(x)=2cos(x)-xsin(x),

Newton's iteration scheme,

xn+1=xn-g(xn)/g'(xn)

takes the form

xn+1=xn-[sin(xn)+x cos(xn)]/[2cos(xn)-xsin(xn)].

To check out in which range the root is, we first plot g(x) in the range 0£x£2.5 using the command

>> ezplot('sin(x)+x*cos(x)',[0 2.5])

The graph (not shown) indicates that there is a zero in the range 1.5£x£2.5. We use Newton's iteration with a starting value in that range to approximate the root. We are happy if the difference between successive iterates is smaller than 10-5, |xn+1-xn|£10-5. We also do not want that the number of iterations gets too large, so we restrict n to n£25. Moreover, we want to know all the iterates in the sequence until the iteration is stopped.

If we use x0=1.5 as starting value, this iteration algorithm can be implemented in Matlab through the following commands (executed in a script):

x=1.5;                                         %set starting value
nmax=25;                                       %set max number of iterations
eps=1;                                         %initialize error bound eps
xvals=x;                                       %initialize array of iterates
n=0;                                           %initialize n (counts iterations)

while eps>=1e-5&n<=nmax                        %set while-conditions
    y=x-(sin(x)+x*cos(x))/(2*cos(x)-x*sin(x)); %compute next iterate
    xvals=[xvals;y];                           %write next iterate in array
    eps=abs(y-x);                              %compute error
    x=y;n=n+1;                                 %update x and n
end                                            %end of while-loop

The iterates are stored in the array xvals. The last value in this array is the desired approximation of the root and coincides with last update of x,

>> xvals

xvals =

   1.50000000000000
   2.31460495577897
   2.04271147451151
   2.02882064418941
   2.02875783941878
   2.02875783811043

>> x

x =

   2.02875783811043

If you are not interested in the full sequence of iterates, then the array xvals is not needed and the code can be shortened as follows:

x=1.5;
nmax=25;
eps=1;
n=0;

while eps>=1e-5&n<=nmax
    y=x-(sin(x)+x*cos(x))/(2*cos(x)-x*sin(x));
    eps=abs(y-x);
    x=y;n=n+1;
end

Here the last value of the iteration is stored in the variable x, but the previous iterates are no longer accessible.