Newton's Method returns complex value.

So, I've written a program that carries out Newton's method. The root of the equation which I am trying to find is approximately 13.1. This is fine and my program returns the correct value when my initial guess is around this value (up to about x = 30 as my initial guess), however when I start using values such as 100+ it returns a complex root. The real part will be approximately 13.1 and the imaginary part will be VERY close to 0. Why is this and is there any way that I can fix/safeguard against this? Thanks.

function xnew = Newton (f, df, xi, tol)
xold = xi; %x(old) is assigned the value of the initial guess
xnew = xold - f(xold)/df(xold); %Implement newtons method to find x(new)
k = 0; %Assigns k(the counter) an initial value
fprintf('\nTable of Iteration No.(k) and Depth(h)\n')
fprintf('\nIteration No.\tDepth\n')
fprintf('%5u\t\t%2.6e\n',k,xi)

while ((abs(xnew - xold))/(abs(xnew)) > tol) %Running condition
    if (k <= 100) %Max number of iterations
        xold = xnew; %x(old) get's x(new)'s calculated value as per Newton Method's
        xnew = xold - f(xold)/df(xold); 
        k = k + 1; %Increment k
    else
        warning('Maximum number of iterations has been reached') 
        break;
    end
    fprintf('%5u\t\t%2.6e\n',k,xnew)
end

[A Jenkins]

Is it possible your f and df are returning imaginary values? Try printing them out.

Remember with this algorithm it is possible to have "weird" values of x (negative, 0, etc), so you may also want to take a look at the xnew values it has printed out.

Functions like abs() or real() can hide the fact that you are getting imaginary results, but there are drawbacks in terms of linearity if you choose to do this.

[David]

Yeah, so f and df are both printing out complex values once the initial guess is about 45 and above(remember the actual value of the root is about 13.1) and once the initial guess is about 1000, NaN is returned as the root. This is due to f being Inf + some tiny imaginary number and df being -Inf + some tiny imaginary number. Thus, f/fd will return NaN. First of all, is this due to a problem with my algorithm or is this just a shortcoming of Newton's Method? And second, how can I safeguard against getting imaginary values and values that are too large(i.e. Inf) for f and df? Like, is there a way to check before Newton's Method is carried out to see if the initial guess is valid? Thanks.

[Walter Roberson]

Imagine a function that looks nice and smooth and gives every indication that you can go ahead and extrapolate a value. But at some place a little before the location that would be extrapolated at, put in a discontinuity, or make the function go non-real. Newton's Method will not be able to deal with that function.

Answer 1

[A Jenkins]

If this is a school project requiring that you use Newton's method on some non-smooth function, than your professor probably is looking for you to notice this limitation of Newton's Method, as Walter described above.

MATLAB has lots of other built in functions that allow you to set bounds, or to try to solve functions that are non-differentiable, so if you are interested, you can take a look at those:

The regular Optimization Toolbox

http://www.mathworks.com/help/matlab/optimization.html

fminbnd()
fminsearch()

or the Global Optimization Toolbox

http://www.mathworks.com/products/global-optimization/examples.html?file=/products/demos/shipping/globaloptim/nonSmoothOpt.html

[David]

It's a school project. The function itself isn't defined for certain values of x, it contains a fraction with x in the denominator, so there's definitely values where it isn't defined properly. As well as that, there's an exponential aspect to it too. As x approaches +Inf, f(x) approaches -6 or so and as x approaches -Inf, f(x) approaches +Inf. So, if you choose an initial guess large enough, the tangent to f(x) may never hit the x-axis (or if it does it'll be at a very large negative value). So, Newton's method won't really work if these values of x are chosen as the initial guess. See what I'm saying? So, that clears that up, in my head anyway. I still have absolutely no clue where the imaginary numbers are coming from. Could you explain how I could use the above functions you stated to help with the problem of choosing too large an x? Fairly new to Matlab, so I'm not all that competent in it.

[Walter Roberson]

What are you passing in for f and df ?

[David]

Is it possible to post mathematical equations on here in LaTeX form? Some sites allow that. Otherwise, they're fairly long, and may be tough to follow.

[Walter Roberson]

LaTeX is not directly supported here. There is at least one site that allows you to paste in a LateX expression, and it gives you back a URL that when viewed gives a rendered version of it. Other than that, you can take a snapshot and post that.

With Newton's Method, as soon as f(x) or df(x) returns a complex value, all subsequent calculations will also return complex values (rather it is very very likely.) So we would need to examine f(x) and df(x) in order to determine wear the complex values are coming from.

[A Jenkins]

Per your question about how to use the toolboxes (assuming you have them):

Intead of calling:

xnew = Newton (f, df, xi, tol)

Try calling

xnew = fminsearch(f, xi, optimset('Display','iter'))

(f is the handle to your function, xi is the initial condition, display makes it print the iterations)

Or, if you like GUIs, type

optimtool

(You'll see there are lots of options. Sometimes the defaults work, but if not, then you have to start digging into your function "f" to understand a suitable way to set it up.)

[David]

@Walter Robertson. I'll just post the Matlab code iself so, with all the defined constants. Should be easier to analyse in Matlab that way.

rhow = 1025;
g = 9.81;
P0 = 101325;
gamma = 1.4;
m = 67;
Vbody = 0.06;
Vlungs0 = 0.004;
tol = 0.01;
%Anonymous Functions:
P = @(h) P0 + rhow .* g .* h; 
Vlungs = @(h) ((P0)./(P(h))).^(1/gamma) .* Vlungs0;
Vsuit = @(h) 0.005 .* exp(-0.035 .* h);
V = @(h) Vbody + Vlungs(h) + Vsuit(h);
f = @(h) (-m) + rhow .* V(h);
%Derivatives of Anonymous Functions:
Pd = @(h) rhow * g;
Vlungsd = @(h) -(P0 .* Vlungs0 .* g .* rhow .*(P0./P(h)).^(1./gamma - 1))./(gamma.*(P(h)).^2);
Vsuitd = @(h) -1.7500e-04 .* exp(-0.035 .* h);
Vd = @(h) Vlungsd(h) + Vsuitd(h);
fd = @(h) rhow .* Vd(h); 

h = 0:25; %Chosen domain for f(h).
%Graph:
figure
hold on
plot(h,f(h),'ro', 'MarkerSize', 5)
plot([0,25], [0,0] , 'k-'); %Plots x = 0 line to better display root.
grid minor
title('h vs. f(h)')
xlabel('h')
ylabel('f(h)')

[David]

Also, if that's somewhat hard to follow, let me know and I'll try format it in another way.

[David]

Any luck Walter Robertson?

[Walter Roberson]

If h goes negative then P(h) could go negative and then P0/P(h) would be negative, and you would be raising that negative value to a power, which is going to give you a complex result.

You do not show the line invoking your Newton function. Keep in mind that your Newton function is not going to be restricted to invoking the function handles on h in the range you assign in your script (h = 0:25)

Answer 2

[David]

@A Jenkins. Thanks for that. I'll have a look into that now.