MATLAB Answers

Newton Raphson Method not solving

6 views (last 30 days)
Max House
Max House on 13 Apr 2021
Commented: Max House on 13 Apr 2021
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

Accepted Answer

Alan Stevens
Alan Stevens on 13 Apr 2021
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.
  4 Comments
Max House
Max House on 13 Apr 2021
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.

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 13 Apr 2021
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.
  1 Comment
Max House
Max House on 13 Apr 2021
Thanks, that's been helpful! I will go back and check for errors.
Max

Sign in to comment.

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by