Finding nearest left and right root using Newton-Raphson

Question

The following is a problem in a Numerical Methods class.

Consider the following function:

$$f(x) = x^3 - \frac{31}{10}x^2+\frac{1}{10}x+\frac{21}{5}$$

Implement a method that will allow you to find the nearest left and right roots, using the Newton-Raphson method as your base root finding method.

You should find the nearest roots of the following start points: $x = 0.0161, x= 2.051, x=0.5$

This is to be programmed in Matlab, so I'm looking at making a method which can be iterated over.

So far, I find the first root of the function (be it left or right), by continually evaluating:

$$x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$$

Starting with $x = 0.0161$, I find the closest root lies to the left at $x= -1$.

Finding the other root, to the right, is where I get a bit lost.

What I've concluded is that if I find the first root to the left, then the result of the following $f(x)f'(x) > 0$ will be true. And for a root to the right, the inverse, $f(x)f'(x) < 0$ will be true.

However, since my start value of $x$ is equal to $0.0161$, the derivative evaluates at $9.5763×10^{-4}$ and the jump that I make with Newton-Raphson is very big because of the small denominator.

Mathematically, how can I arrive at the other root?

Answer 1

I am not sure if this is a very good solution, but it works at least for your problem. The standard Newton-Raphson method uses the linear approximation of the function. One could also use a quadratic approximation. This quadratic approximation can have two solutions, upon which one choose either of these solutions to further iterate using the standard Newton-Raphson method. There is, however, a possibility that the quadratic approximation has no solutions.

So let say that we start at $x_0$. The quadratic approximations has the following solution: $$ x = x_0 + \frac{-f'(x_0) \pm \sqrt{(f'(x_0))^2 - 4f(x_0)f''(x_0)}}{2f''(x_0)}. $$ From the two solutions found, one can use the standard Newton-Raphson method: $$ x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}. $$

I used the following code to check whether it works for your example:

% Define function and derivatives
f = @(x) x^3 - 31/10*x^2 + 1/10*x + 21/5;
fdot = @(x) 3*x^2 - 31/5*x + 1/10;
fddot = @(x) 6*x - 31/5;

% Starting point
x0 = 2.051;

% Quadratic approximation
a = fddot(x0);
b = fdot(x0);
c = f(x0);
D = b^2 - 4*a*c;
if D < 0
    disp('Quadratic function does not have any real roots')
    xx = x0;
else
    xx = [(-b-sign(a)*sqrt(D)) / (2*a), (-b+sign(a)*sqrt(D)) / (2*a)] + x0;
end;

% Newton-Raphson
for i = 1:numel(xx)
    for iter = 1:100
        xx(i) = xx(i) - f(xx(i)) / fdot(xx(i));
    end;
end;
if D < 0
    fprintf('Single solution: %.3f\n', xx);
else
    fprintf('Left solution: %.3f\n', xx(1));
    fprintf('Right solution: %.3f\n', xx(2));
end;

x0 = 0.0161 results in:

Left solution: -1.000
Right solution: 2.000

x0 = 2.051 results in:

Left solution: 2.000
Right solution: 2.100

x0 = 0.5 results in:

Left solution: -1.000
Right solution: 2.000

Answer 2

Following through with Tyberius's comment about factoring each successive root found out of the original polynomial led to the following code:

close all; clear *; clc;

x1 = 0.0161;
x2 = 2.051;
x3 = 0.5;
x4 = -2;
x5 = 3;

figureCount = 1;

Xs = [ x1 x2 x3 x4 x5 ]
for ( K = 1:size(Xs,2))


    syms x;


    funF = x^3 - (31/10)*x^2 + (1/10)*x + (21/5);
    funDF = 3*x^2 - 2*(31/10)*x + (1/10);

    x = Xs(K);
    f = subs(funF);
    df = subs(funDF);

    Xnext = x - (f/df);
    tolerance = 0.00000001;
    iterations = 100;
    i = 0;
    Xoriginal = x;
    solutions = [];
    numberOfSolutions = 0;
    flag = 0;
    originalLeft = 0;
    originalRight = 0;
    left = 0;
    right = 0;

    originalF = funF;
    originalDF = funDF;
    title('Newton Raphson Root Finding');
    xlabel('Iteration number (i)');
    ylabel('Value of Xi');


    while ( i <= iterations )
        if ( left == 1 || right == 1)
           flag = 0;
           right = 0;
           left = 0;
        end

        if (numberOfSolutions > 0 && flag == 0)             %% Flag controls that we only enter this block when we need to factor out a root
            flag = 1;
            syms x;
            currentSolution = solutions(end);
            funF = simplify(funF/(x - currentSolution));
            funDF = diff(funF);
            x = Xoriginal;
        end

        f = eval(subs(funF));
        df = eval(subs(funDF));
        if ( df == 0)
            df = tolerance;
        end

        ppd = f/df;
        Xnext = x - (ppd);
        Error = abs((Xnext - x)/Xnext);
        x = Xnext;

        if ( Error <= tolerance)
            if (x < Xoriginal ) %% we went left
                if (numberOfSolutions == 0)
                    originalLeft = 1;
                end
                left = 1;
                right = 0;
            else
                if ( x > Xoriginal ) %% we went right
                    if (numberOfSolutions == 0)
                        originalRight = 1;
                    end
                    right = 1;
                    left = 0;
                end
            end
            solutions = [solutions x]; %% solutions is equal to itself + next solution (x)
            numberOfSolutions=size(solutions,2); %% get size of solution set
            flag = 0;
            if( originalLeft == 1 && right == 1) 
                    %%%finish root finding
                break;
            else
                if ( originalRight == 1 && left == 1 )
                    %%%finish root finding
                    break;
                end

            end
        end
        i = i + 1;
        hold on
        stem(i, x)
    end
    if figureCount < size(Xs,2) %% make a figure for each value of X
        figureCount = figureCount + 1;
        figure;
    end

    leftSolutions = [];
    rightSolutions = [];

    %%% Build right solutions and left solutions matrix
    for  J = 1 : numberOfSolutions 
        if ( solutions(J) < Xoriginal)
            leftSolutions = [leftSolutions solutions(J)];
        else
            rightSolutions = [rightSolutions solutions(J)];
        end
    end

    leftSolution = max(leftSolutions);
    rightSolution = min(rightSolutions);
    if( ~isempty(leftSolution))  
        disp(strcat("Zero 1 found to the left of ", num2str(Xoriginal), " with value of ", num2str(leftSolution), " for function ", char(originalF))); 
    end
    if ( ~isempty(rightSolution))
        disp(strcat("Zero 2 found to the right of ", num2str(Xoriginal), " with value of ", num2str(rightSolution), " for function ", char(originalDF)));
    end

    disp('-------------------------------------');
end

If in the case no root is found to the left or right, the program has an iteration cap of 100.