3
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Tyberius Mar 3 '18 at 22:32
  • $\begingroup$ @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. $\endgroup$ – Michael McQuade Mar 4 '18 at 7:28
  • $\begingroup$ Related: math.stackexchange.com/q/1644377/418145 $\endgroup$ – Tyberius Mar 4 '18 at 19:46
  • $\begingroup$ the idea offered by the related question is to combine Newton's method with something like bisection. $\endgroup$ – Tyberius Mar 4 '18 at 19:58
2
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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.

$\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.