HiveBrain v1.2.0
Get Started
← Back to all entries
patternMinor

Visualizing Newton-Raphson method for finding zeroes of a function

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
visualizingmethodraphsonfunctionnewtonforfindingzeroes

Problem

I have created a program to visualize the working of Newton-Raphson method to find the zeroes of a function:

newton.m:

``
% Dummy statement to avoid writing function in the first line and making it a 'function file' instead of a 'script file'
1;

% The function to find zeroes of.
% The function is specifically chosen to not have any zeroes
% so as to show the weakness of Newton's method.
function y = f(x)
y = (x - 5).^2 + 5;
endfunction

% The derivative of f(x)
function y = fd(x)
y = 2 * (x - 5);
endfunction

% Initial guess
x0 = 1.5;

% Max number of iterations
itermax = 20;

% Epsilon value initialized to a very large value
eps = 1;

% A vector for storing the history of the approximate roots
xvals = x0;

% Number of iterations done
itercount = 0;

% Required for plotting f(x) vs x
x = linspace(0, 10, 100);

% Create a figure whose output is not rendered on the screen
% Not working currently; supposedly a bug in Octave
% A workaround is to use gnuplot instead of qt -
graphics_toolkit gnuplot`
% but this is very slow.
% Uncomment the following to activate the feature once the bug is fixed
% figure('Visible','off');

% The main loop
while eps >= 1e-5 && itercount <= itermax
% x1 = New value of root
% x0 = Current value of root
x1 = x0 - f(x0) / fd(x0);

% Plot f(x)
% Plot a line passing through points [x0, f(x0)] and [x1, 0]
% Plot a line passing through points [x1, 0] and [x1, f(x1)]
% Plot a line passing through points [x0, 0] and [x0, f(x0)]
plot(x, f(x), ";f(x);", [x0 x1], [f(x0) 0], "-r;f'(x);", [x1 x1], [0 f(x1)], ":r", [x0 x0], [0 f(x0)], ":r");
title('f(x) = (x-5)^2 + 5');

% Set limits for the axes shown in the plots
xlim([0 10]);
ylim([0 30]);

% Label the two consecutive zeroes on the X-axis
text(x0, -2, sprintf('x%d', itercount), 'color', 'red');
text(x1, -2, sprintf('x%d', itercount+1), 'color', 'red');

% Print the plot to a file
filename = sprintf('output/%05d.jpg', ite

Solution

Conventions, coding patterns and readability

Overall, your code looks nice! I have never seen a code where you insert a line in the beginning to avoid the script being interpreted as a function. It works, but I wouldn't do it like that.

The two functions you have are very simple. There is 1: no recursion. 2: no loops 3: no handling of invalid input. For such functions you can simply create anonymous functions. This way you can avoid the first line of code in your script.

% The function to find zeroes of.
% The function is specifically chosen to not have any zeroes
% so as to show the weakness of Newton's method.
f = @(x) (x - 5).^2 + 5; 

% The derivative of f(x)
fd = @(x) 2 * (x - 5);


You have very descriptive variable names. However, eps is not one of them, as eps is a built-in function and variable.

You have a growing vector in your loop xvals = [xvals; x1];. This is not recommended, as it is very slow. It's faster to pre-allocate memory for a vector that will be "large enough", and trim it down afterwards:

xvals = zeros(20,1);
while ...
...
end
xvals = xvals(1:iter);  % or xvals(iter+1:end) = [];


Also, if you don't want to pre-allocate memory, you can increase performance by substituting xvals = [xvals; x1] with:

xvals(end+1) = x1;


You should use

  • Don't do a = b = c = 0 because it's easier.



  • Don't do inline assignments b = (a=3) * a just because you can.



  • And it might not be recommended by everyone, but I recommend you to use end instead of endfunction, endif` etc to keep compatibility with MATLAB.



This is all I can do for now... Good luck!

Code Snippets

% The function to find zeroes of.
% The function is specifically chosen to not have any zeroes
% so as to show the weakness of Newton's method.
f = @(x) (x - 5).^2 + 5; 

% The derivative of f(x)
fd = @(x) 2 * (x - 5);
xvals = zeros(20,1);
while ...
...
end
xvals = xvals(1:iter);  % or xvals(iter+1:end) = [];
xvals(end+1) = x1;

Context

StackExchange Code Review Q#125975, answer score: 5

Revisions (0)

No revisions yet.