# Newton Raphson Method not solving

28 views (last 30 days)
Max House on 13 Apr 2021 at 11:55
Commented: Max House on 13 Apr 2021 at 15:39
I am a mechanical enginering student trying to solve an equation frequiring a newton raphson iteration. However, when I run the programme, I either get the initial guess or NaN + InfI as my solution, depending on the tolerance chosen. The expected answer is a real number. Below is my script:
function x = NewtonRaphson(f,df,x0,tol)
%Using Newton-Raphson with tolerance
x = x0;
y = f(x);
while abs(y)>tol
x = x - (f(x)/df(x));
y = f(x);
end
end
And command:
>> c3=-0.56;
>> c4=0.182;
>> c5=-0.09;
>> c6=0.0102;
>> Et=0.021;
>> fn = @(x) c4*(2*x).^(c3) + c6*(2*x).^(c5) - Et;
>> deriv = @(x) (c3)*(c4)*(2*x).^(c3-1) + (c5)*(c6)*(2*x).^(c5-1);
>> NewtonRaphson(fn,deriv,1000,0.01);
(The equation I'm trying to solve is the Coffin-Manson-Basquin equation, where 2x is number of reversals, if that helps)
Many thanks,
Max

Alan Stevens on 13 Apr 2021 at 12:45
Try
c3=-0.56;
c4=0.182;
c5=-0.09;
c6=0.0102;
Et=0.021;
fn = @(x) c4*(2*x).^(c3) + c6*(2*x).^(c5) - Et;
deriv = @(x) (c3)*(c4)*2^c3*x.^(c3-1) + (c5)*(c6)*2^c5*x.^(c5-1); %%%%%
x = NewtonRaphson(fn,deriv,10,0.01); %%%%%
disp(x)
function x = NewtonRaphson(f,df,x0,tol)
%Using Newton-Raphson with tolerance
x = x0;
y = f(x);
while abs(y)>tol
x = x - (f(x)/df(x));
y = f(x);
end
end
Note that when the function is differentiated, only the x's drop a power, not the 2's.
Also, the shape of the function is such that it's best to give an initial guess on the low side of the root.
Max House on 13 Apr 2021 at 15:34
Thank you @John D'Errico and @Alan Stevens , your answers have been very helpful. I suspect I must have made a mistake somewhere! Thanks again.

John D'Errico on 13 Apr 2021 at 15:04
This answer is here just to add some information. I try to teach that one should always plot everything, then look at what they see.
c3=-0.56;
c4=0.182;
c5=-0.09;
c6=0.0102;
Et=0.021;
fn = @(x) c4*(2*x).^(c3) + c6*(2*x).^(c5) - Et;
fplot(fn,[0,2000])
yline(0);
grid on Here, the start point is 1000, but the root to be found appears to be much closer to zero, probably near 50 I'd guess from the plot. And that is a clear problem, because there is a singularity at x==0.
In fact, below x == 0, if you try to evaluate fn, you will find it produces complex numbers. That is because of the negative fractional powers of x.
fn(-1)
ans = -0.0349 - 0.1239i
If we view Newton's method as approximating the function by a tangent line at the iteration point, then where will that first step send us?
x0 = 1000;
deriv = @(x) (c3)*(c4)*2^c3*x.^(c3-1) + (c5)*(c6)*2^c5*x.^(c5-1); %%%%%
fn(x0)
ans = -0.0133
deriv(x0)
ans = -1.9076e-06
So the locally tangent approximation to fn has a VERY tiny and negative slope at that point. What ends up happening is that line now forces the Newton step to be a bad idea.
xnew = x0 - (fn(x0)/deriv(x0))
xnew = -5.9588e+03
So after the very first step, we are in negative territory for x. And once you go into that complex rabbit hole, Newton's method will not escape alive.
This is why you should ALWAYS be using a better tool to solve your problems. Newton's meth0d is a poor choice here. Instead, use fzero. For example...
[xfinal,fval,exitflag] = fzero(fn,[eps,1e6])
xfinal = 47.3820
fval = 3.4694e-18
exitflag = 1
So fzero is entirely happy. See that I gave it a nice wide bracket, but one that forces the method to stay above zero. It found the expected root at around 50, so my wild guess was quite a good one. But this brings up a point that as I said, Newton is a terrible choice of method for this shape of a function with a root near a problematic singularity as you have here. As long as you start above the root, the first iteration often will blast well below zero. For example, had you started with an x0 of 200, it still wil fail.
x0 = 200;
xnew = x0 - (fn(x0)/deriv(x0))
xnew = -225.1279
Only if I start sufficiently close to the root does it now converge. 60 seems to work.
x0 = 60;
for i = 1:5
x0 = x0 - (fn(x0)/deriv(x0))
end
x0 = 44.9199
x0 = 47.2835
x0 = 47.3819
x0 = 47.3820
x0 = 47.3820
You stated in a comment that you expected a root in the vicinity of 1000 for this function. While the plot does not look as if that is true, we might try expanding the region of the plot a bit, to see if a root will ever exist.
fplot(fn,[2000,100000]) And just looking at that figure, my gut tells me it will never turn around, in order to eventually cross zero. As you can see, going all the way out to 1e5 sees no root beyond the one at x==47.
Can we prove that a root will never exist above a certain point, that this curve will never turn around? That is probably doable, though much easier with positive integer powers of x. Then I could use a tool like Descarte's rule of signs.
My gut tells me not to bother, and this answer has now gotten long, so I won't.
Max House on 13 Apr 2021 at 15:39
Thanks, that's been helpful! I will go back and check for errors.
Max

R2020a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!