Below is a simple MATLAB m-file to solve sin(x) = 0:
function [x,fail] = newton(f,df,x,eps,N)
%
% Applies Newton's method to the function f, with derivative in df,
% starting at the point x. It produces a sequences of points.
% Stops when the sequence reaches the point x where
% |f(x)| < eps or when two consecutive points in the sequence
% have distance less than eps apart from each other.
% Failure is reported (fail = 1) when the accuracy requirement
% is not satisfied in N steps; otherwise fail = 0 on return.
%
% Example :
% >> [x,fail] = newton('sin','cos',pi/4,1e-8,10)
%
fprintf('running the method of Newton...\n');
fx = feval(f,x);
for i = 1:N
dx = fx/feval(df,x);
x = x - dx;
fx = feval(f,x);
fprintf(' x = %e, dx = %e, f(x) = %e\n', x, dx, fx);
if (abs(fx) < eps) | (abs(dx) < eps)
fail = 0;
fprintf('succeeded after %d steps\n', i);
return;
end
end
fprintf('failed requirements after %d steps\n', N);
fail = 1;
The output of this algorithm on sin(x) = 0 starting at pi/4 is
>> [x,fail] = newton('sin','cos',pi/4,1e-8,10)
running the method of Newton...
x = -2.146018e-01, dx = 1.000000e+00, f(x) = -2.129584e-01
x = 3.356262e-03, dx = -2.179581e-01, f(x) = 3.356256e-03
x = -1.260225e-08, dx = 3.356274e-03, f(x) = -1.260225e-08
x = 1.654361e-24, dx = -1.260225e-08, f(x) = 1.654361e-24
succeeded after 4 steps
x =
1.6544e-24
fail =
0
Observe the fast convergence.
We can interpret Newton's method as a fixed-point iteration method. This leads to "cobweb pictures". The m-file below produces such a picture with MATLAB:
function [x,fail] = cobweb(g,x,eps,N,xmin,xmax,ymin,ymax)
%
% On entry :
% g function defining the x = g(x);
% x start point for the fixed-point iteration;
% eps accuracy requirement: stops when two consecutive
% approximations in the sequence are
% N maximal number of allowed iterations;
% xmin minimal value for x to view the graph;
% xmax maximal value for x to view the graph;
% ymin minimal value for g(x) to view the graph;
% ymax maximal value for g(x) to view the graph.
%
% On return :
% x last approximation in the sequence;
% fail 0 means success, 1 means failure.
%
% Example :
% >> cobweb('myfun',0.3,1.0e-6,10,0.1,10,0,10)
%
fprintf('illustration of fixed-point iteration...\n');
xr = xmin:0.1:xmax; % sample range to plot
yr = feval(g,xr); % evaluate at the samples
figure; % new figure window
axis([xmin xmax ymin ymax]); % adjust the scale
hold on; % overlay plot
plot(xr,yr); % plot the function y = g(x)
plot(xr,xr); % plot the diagonal y = x
for i = 1:N % perform N iterations
y = feval(g,x); % y = g(x)
plot([x x],[x y],'g--'); % vertical line
plot([x y],[y y],'r--'); % horizontal line
fprintf(' x = %e, y = %e, |x-y| = %e\n', x, y, abs(x-y));
if (abs(x-y) < eps)
fprintf('succeeded after %d steps\n', i);
x = y;
fail = 0;
return; % terminate with success
end;
x = y;
end;
fail = 1; % terminate with failure
where the content of myfun.m is
function y = myfun(x) % % Returns the function value y = (x^2 + 5)/(2*x). % The dots before the exponentiation and division operator below % enforce the operation to be carried out componentwise. % These dots are necessary when x is a vector. % y = (x.^2 + 5)./(2*x);
Note that with myfun is obtained from applying Newton's method on the equation x^2 - 5 = 0.
The nongraphical output is
>> cobweb('myfun',0.3,1.0e-6,10,0.1,10,0,10)
illustration of fixed-point iteration...
x = 3.000000e-01, y = 8.483333e+00, |x-y| = 8.183333e+00
x = 8.483333e+00, y = 4.536362e+00, |x-y| = 3.946971e+00
x = 4.536362e+00, y = 2.819283e+00, |x-y| = 1.717079e+00
x = 2.819283e+00, y = 2.296392e+00, |x-y| = 5.228916e-01
x = 2.296392e+00, y = 2.236860e+00, |x-y| = 5.953157e-02
x = 2.236860e+00, y = 2.236068e+00, |x-y| = 7.921835e-04
x = 2.236068e+00, y = 2.236068e+00, |x-y| = 1.403255e-07
succeeded after 7 steps
ans =
2.2361
>> ans^2
ans =
5.0000
>>
The picture we see is