# Finding nearest left and right root using Newton-Raphson

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?

• Once you have found a root, you can factor it out of your polynomial and start again from the same initial point. I'd be interested to hear a better way to force Newton's method to go in one direction, because this way would be a bit of a hassle. – Tyberius Mar 3 '18 at 22:32
• @Tyberius This method actually works quite well, but in some cases we won't find the nearest root on the right, but some other root. I wonder if there is a way to control it to only find the closest roots. – Michael McQuade Mar 4 '18 at 7:28
• – Tyberius Mar 4 '18 at 19:46
• the idea offered by the related question is to combine Newton's method with something like bisection. – Tyberius Mar 4 '18 at 19:58

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;

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


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.