6 views (last 30 days)

Show older comments

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

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.

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)

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)

deriv(x0)

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))

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])

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))

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

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.

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

Start Hunting!