7
\$\begingroup\$

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', itercount);
    print(filename)

    % Append the zero to the array of zeroes calculated so far
    xvals = [xvals; x1];

    % Calculate the epsilon value
    eps = abs(x1-x0);

    x0 = x1;
    itercount = itercount+1;
end


% Print the result of the iteration
xvals
f_zero = f(xvals(end))
eps
itercount

This program produces a series of images which I have combined to create an animation:

Simulation of Newton's Raphson method on a function not having any zeroes

The command I used to convert the images is:

convert -delay 100 -loop 0 *.png animation.gif

This code is also stored here.

I would like to get it reviewed on the following:

  • Octave conventions and coding patterns (This is my first programming attempt using Octave)
  • Any improvement on the algorithm itself
  • Readability and performance of the code
  • Method and quality of producing animations
\$\endgroup\$
5
\$\begingroup\$

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 < itermax instead of <= itermax in the loop condition. Otherwise you'll end up with itercount = itermax + 1.


Note: Always try to keep compatibility with MATLAB when programming in Octave.

  • Don't use itercount ++ instead of itercount = itercount + 1, because it's easier.
  • 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!

\$\endgroup\$
3
\$\begingroup\$

Although I'm with Stewie Griffin in many points I want to highlight where I disagree.

Octave tries to complete the Matlab language and therefore constructs like itercount++, a = b = c = 0 (which should be familiar to every C rogrammer) were added. If you don't have to be Matlab compatible just use them. They make your code easier to read and more fun to write. In my opinion the endfor, endwhile, endfunction keywords also increase the readability.

Btw, both graphics toolkits qt and fltk can use OSMesa for offscreen rendering if you use a Mesa based GPU driver on GNU/Linux. (your "bug" comment)

Octave conventions and coding patterns:

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");

  • Create the output directory before blindly trying to write to it: mkdir ("output");
  • I would keep the last 3..5 steps in the plots and/or let the last steps vanish (decrease saturation)
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.