How do I create a loop that will repeat my Newton Raphson algorithm until all the roots are found, in matlab?

Question

I wrote the following code which uses the Newton Raphson method to find roots. It works to find 1 root, but then it breaks. How should a modify the code to repeat the algorithm until all the roots are found in the desired range?

I know I should use a for or while loop I just cant figure it out.

'''

function [xn nrfail] = newraphson(fun,xg,xl,xr,tol)
% fun is the function of interest = (0.01-1) + x*cot(x), xg is the initial root guess, xl and xr are the bounds of the range I am using
% Initializing
i=0;
nrfail=0;
check=1;
h=1e-4;

% Loop
    while tol<check % Tolerence check
        i=i+1;
        fp=derivative(fun,xg,h); %derivative is another function I wrote that finds the derivative of the function of interest
        y=feval(fun,xg);
        xn=xg-y/fp; % NR method

        if xn<xl || xn>xr || i>10 || fp==0 % Method check
            nrfail=1;
            break
        end

        if abs(xn) <= 1
            check=abs(xg-xn); % Absolute error
        else
            check=abs(1-xg/xn); % Relative error
        end

        xg=xn; % New guess
        end
    end

'''

Initial conditions are crucial, and its not trivial how to choose them, to the point that this algorithm creates a fractal. Check Newtons fractal. You can indeed fall into an infinite loop of solutions here. — Ander Biguri, Oct 13 '22 at 22:21
So to avoid an infinite loop would it be more efficient to add a for loop outside of the code ? — cschlum, Oct 13 '22 at 22:35
So, as far as I understand, for a given initial condition, your code works, right? but then your polynomial may have more roots. So, the only thing you can really do here is go trying different initial conditioins, and just hope that you will find a new root. But the Newton fractals show you that there is no trivial way of chosing these initial conditions. It gets inifnitely complex. So unless you have more information on the specifics of your equations, the answer is that you can't find a general solution to chose starting conditions such that it ensures you will find all roots. — Ander Biguri, Oct 14 '22 at 13:18

score 0 · Answer 1 · edited Dec 31 '22 at 01:14

First of all, there is no way to be sure that the Newton-Raphson method will converge for an initial input. The algorithm converges once the derivative of the function reaches abs (f') <= errThreshold. This means that the algorithm fails for

Very 'flat' functions or
Regions in the function domain where the derivative gets flat, i.e. functions that are locally flat.

A second issue is that for functions with multiple maximas or minimas, there will be multiple points where the function flattens. So the outcome will actually depend on the starting point. I give these initial warnings because unless you are totally 100% sure that the function is very well defined, you could well end up in an infinite loop. So the best approach would be to, instead of setting an infinite loop, allow a maximumNumberOfIterations, after which the algorithm 'breaks down' and returns 'could not find some of the roots' ...

That said depending on what you are calculating, consider using the bissection method instead. It tends to have a slower convergence, but then again, the convergence is less dependent on the derivative, so it will tend to converge pretty much in all cases.

The algorithm you are looking for is the following (on Python):

from numpy.random import rand

#
# The maximum number of iterations before the 
# algorithm breaks.
##
maximumNumberOfIterations : int

#
# The function that we are evaluating.
##
function : Callable

#
# The lower bound.
##
lowerBound : float;

#
# The upper bound.
##
upperBound : float;

#
# The tolerance.
##
tolerance : float;

#
# An array to store the outcomes 
# for each iteration
##
roots: list[float] = [] 

#
# For each iteration DO:
##
for iteration in range (maximumNumberOfIterations):

    #
    # Generate a 'random' initial guess
    # specified within the allowed boundary
    ##
    initialGuess : float = (upperBound - lowerBound) * rand()

    #
    # Actually try to run the algorithm and append
    # the root to the list.
    ##
    try:
        roots.append (
            newraphson(function, guess, lowerBound,
            upperBound, tolerance)
        )

    #
    # It is likely that the algorithm will fail and
    # return a RuntimeFault for 'some' cases. We
    # must handle these... Not a big deal, just do
    # nothing and continue.
    ##
    except RuntimeError:
        pass

#

# Now we just need to handle the repeated roots. Being
# lazy I can just use a set and display the results...
##
roots = set (roots)

for root in roots:
    print (roots)
    
#
# If you were expecting a certain number of roots ...
##

expectedNumberOfRoots : int

if len (roots) < expectedNumberOfRoots:
    print ({n} roots could not be found"
    .format (n=expectedNumberOfRoots-len (roots)))

How do I create a loop that will repeat my Newton Raphson algorithm until all the roots are found, in matlab?

1 Answers1