Newton-Raphson requires the calculation of the derivative. I’ve chosen to do so by using the central-difference formula.
The comments should make the files self explanatory. The output at the bottom of the page suggests that my function solves equations in about the same time as MATLAB’s built in fzero. Remember that Newton-Raphson only returns complex results if the initial guess is complex.

Remember that the indentation is not as it should be, this is because of the bad support WordPress has for MATLAB code.  If you wish to play around with the code, DOWNLOAD rar file here.


% The Newton-Raphson method for solving equations of the form f(x)=0.
% We use the central difference approximation of the derivative by calling
% the function central_diff.m.
% ------------------------------------------------------------------------
% INPUT: Anonymous function f, and the first guess at a solution, x0.
% e_diff is the number of significant digits wanted when calculating the
% derivative. For example 1e-5 gives 5 significant digits. e_newton is the
% accuracy of the result.
% OUTPUT: Approxomation of the root, r. Number of iterations, n.

function [r,fval,n] = newton(f,x0,e_diff,e_newton)

n = 0; % Variable used to count the number of iterations.
x1 = x0 - f(x0)/central_diff(f,x0,e_diff); % First approximation.

while 1==1
x2 = x1 - f(x1)/central_diff(f,x1,e_diff); % 'second' approximation.
% We wish to stop if the difference |x_{i+1}-x_i| is less than some
% predefined tolerance e_newton. We also want the value of f at x2 to
% be less than this tolerance.

if abs(x2-x1) < e_newton && abs(f(x2)) < e_newton
r = x2;
fval = f(r);
return;
else
x1 = x2;
n = n + 1;
end
end

</pre>
% Approximates the derivative of f(x0) by using the central difference
% formula.
% --------------------------------------------------------------------
% INPUT: Anonymous function f and the point of interest x0. e_diff is the
% amount of significant digits wanted.
% OUTPUT: The derivative of f at x0.
function [df] = central_diff(f,x0,e_diff)

h = 0.1; % Initial step size h.
df1 = (f(x0+h)-f(x0-h)) / (2*h); % First approximation of derivative.
while 1==1
 df2 = (f(x0+h)-f(x0-h)) / (2*h);
 if abs(df1-df2) < e_diff
 df = df2;
 return;
 else
 df1 = df2;
 h = h/2;
 end
end
<pre>

% Script to test the function newton.m. Compares the solution to Matlab's
% function fzero.
format long

x0 = 1;
f = @(x) x*cos(x)+x^3+4.5*x^2+exp(x);
e_diff = 1e-5;
e_newton = 1e-12; % Gives ten significant digits.

disp('Time used by fzero.m:');
tStart1=tic;
[x,fval_matlab] = fzero(f,x0);
toc(tStart1);

disp('Time used by newton.m:');
tStart2=tic;
[r,fval_dafeda,n] = newton(f,x0,e_diff,e_newton);
toc(tStart2);

fprintf('\nThe solution to f(x)=0 by the Newton-Raphson method is %5.15f: \n',r);
fprintf('\nThe value f(r) %5.15f \n',fval_dafeda);
fprintf('\nNumber of iterations needed: %i \n', n);

OUTPUT:

EDU>> newton_test
Time used by fzero.m:
Elapsed time is 0.002823 seconds.
Time used by newton.m:
Elapsed time is 0.002426 seconds.

The solution to f(x)=0 by the Newton-Raphson method is -4.538614468492257:

The value f(r) -0.000000000000004

Number of iterations needed: 85

About these ads