Newton Rapson method code question

I am creating a MATLAB code that takes a function, differential function, initial value, and error tolerance and solves it using Newton-Raphson's method, but it repeats infinitely. Where is the problem?

function [root,iter] = newton_raphson(f,df,x0,tol)
iter=0;
err = 1;
x=x0;
while(err>=tol)
    iter=1+iter;
    xm=x-(f(x)/df(x));
    err=abs(xm-x);
    x=xm;
end
root=xm;
end

[Torsten]

Where is the problem?

We don't know because you didn't supply the inputs to the function.

Answer 1

[John D'Errico]

Actually, this is surprisingly nice looking code. (A true compliment from me, given most of the student codes we see on here!) You should consider adding comments in your work though. My goal is often to have at least one comment line per line of code.

There are some subtle problems though. You start out with err = 1. BAD. What if you supplied a tolerance of 2? Make your code robust to crazy things!

Initially, set

err = inf;

Now your solver will always take at least one step!

Next, it would be better to set up a maximum number of iterations, so have the user pass in a parameter, probably called itermax, that will have the code stop when iter exceeds that max. Good code would probably issue a warning when that happens. I won't add that feature for you, but it should be here.

function [root,iter] = newton_raphson(f,df,x0,tol)
% 1-d Newton-Raphson solver
% [root,iter] = newton_raphson(f,df,x0,tol)
%  f = provided function handle
%  df = provided derivative function handle
%  x0 = start point
%  tol = tolerance

% initialize iterations
iter = 0;
err = inf;
x = x0;

% looping until convergence
while(err>=tol)
    % update iteration counter
    iter = 1 + iter;
    
    % Newton update
    xm=x-(f(x)/df(x));
    
    % set err as the step length taken
    err=abs(xm-x);
    
    % new estimate for x replacing the old
    x=xm;
end
root=xm;
end

We can try it out.

f = @(x) sin(x) - 0.3;
df = @(x) cos(x);
x0 = 0.1;
tol = 1e-12;
disp("Known solution: " + asin(0.3))

Known solution: 0.30469

[root,iter] = newton_raphson(f,df,x0,tol)

root = 0.3047

iter = 4

Now I made no changes at all to your code, beyond adding some comments to make it more readable, and some white space. As you can see, it worked with no problems.

So, what did you do wrong? Perhaps you gave it a function that has no root at all. For example:

f = x^2 + 1;

df = 2*x;

Or perhaps you tried to solve a problem where you gave it a starting point that diverges, or something equally silly. (There are actually subtle ways your teacher could have caused this to happen, that I won't get into.) But we cannot know. Your code was actually correct. All I did was change the initial value for err to inf.

[kjw]

Thank you so much for your answers.

f = @(x) x^3-2*x+2

df = @(x) 3*x^2-2

In the above function, we need to find the root when the initial value is 0/-2x/2. But if I put 0, it doesn't stop running, and if I put -2x, it doesn't run. Can't we fix this? OR If the root can't be found, can't we just make it output an error message?

[John D'Errico]

f = @(x) x^3-2*x+2;
df = @(x) 3*x^2-2;

Nothing is impossible. Well, as long as you don't ask for my favorite sports team to win a game.

fplot(f,[-2,2])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

grid on

Sigh. This is one of those cases where I said your teacher might have decided to throw in a monkey wrench. Teachers can be a pain in the.....

THINK about what happens. That is, think about what Newton's method does, how it works. This is a classic problem with a Newton method.

First, where is the root? It clearly lies around x == -1.75. We can see that from the plot.

Where did you start it? At x==0.

Where is the start point, in terms of the shape of your function? It lives in that little belly of the curve.

Newton's method approximates your function locally with a straight line, tangent to your curve. You need to understand this. Then it extrapolates that line down to y==0.

x0 = 0;
f0 = f(x0);
df0 = df(x0);
f1hat = @(x) (x - x0)*df0 + f0; % The tangent line

figure
fplot(f,[-2,2])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
fplot(f1hat,[-2,2])
yline(0)
grid on

Do you see how it works? So at the first iteration, xm will be

x1 = x0 - f0/df0

x1 = 1

Ok. that works. But now we are on the other side of that "belly". The next stop of Newton's method takes you BACKWARDS. But not far enough.

f1 = f(x1);
df1 = df(x1);
x2 = x1 - f1/df1

x2 = 0

f2hat = @(x) (x - x1)*df1 + f1; % The NEW tangent line

figure
fplot(f,[-2,2])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
fplot(f1hat,[-2,2])
fplot(f2hat,[-2,2])
yline(0)
grid on

Do you see what happens? x at the second iteration is EXACTLY 0. Newton's method will bounce back and forth forever, from 0 to 1 to 0 to 1 to 0 ... Thank your teacher for this.

Instead, think about what happens when you start the iterations below about x == -1. YUP. Now it will probably converge.

All of this MIGHT be a homework problem designed to teach you about a concept called a basin of attraction, or at least introduce the general idea.

When you show up in class, look into the eyes of your teacher. Do those eyes look just a little scary? Would you buy a used car from them? ;-) Seriously, this is a classic problem, one you should learn from, and why it was assigned to you.

By way of comparison...

function [root,iter] = newton_raphson(f,df,x0,tol)
% 1-d Newton-Raphson solver
% [root,iter] = newton_raphson(f,df,x0,tol)
%  f = provided function handle
%  df = provided derivative function handle
%  x0 = start point
%  tol = tolerance
% initialize iterations
iter = 0;
err = inf;
x = x0;
% looping until convergence
while(err>=tol)
    % update iteration counter
    iter = 1 + iter;
    
    % Newton update
    xm=x-(f(x)/df(x));
    
    % set err as the step length taken
    err=abs(xm-x);
    
    % new estimate for x replacing the old
    x=xm;
end
root=xm;
end

[root,iter] = newton_raphson(f,df,-1,1e-12)

root = -1.7693

iter = 9

Yup. No problem, when started from there.

[kjw]

Thank you very much. You're kind!

[kjw]

If I want to draw like this, how should I write the code?

For f(x)=x3−6x2+11x−6, the error tolerances are ϵ=10^-1, 10^-2, 10^-3, 10^-4, 10^-5, 10^-6 , 10^-7, 10^-8, 10^-9, 10^-10, find each solution using the bisection method and Newton's method and draw a graph.

a) x-axis: error range (Log scale)/ ϵ y-axis: number of repetitions iter (Linear scale)

b) x-axis: error range (Log scale)/ ϵ y-axis: convergence solution (root)

[John D'Errico]

I'm sorry. I won't do your assignment for you. But you know exactly what to do. Break large problems down into small ones that you can handle.

1. Write a bisection function (if you don't already have one) that will take an error tolerance, just like the one you have for Newton. You already have a Newton code.
2. Set up a loop, where you call the codes repeatedly, one for each error tolerance. Save the results in vectors.
3. Plot the results as needed. You will need to use the functions loglog and semilogy for those plots as appropriate.

Again. make a big probem into a set of managable small ones. When all the pieces are in place, put it all together.