Error encountered using 'fzero'

Hello, I have a question regarding fzero. my code is as follows:

syms P
N=-P;
a=50; b=20;  h0=50; E=200;n=1.5;
 h1=40;
l=100;
I0=b*h0^3/12; I1=b*h1^3/12;
L0=sqrt(P/(E*I0-P*n));%L0=simplify(L00);
L1=sqrt(P/(E*I1-P*n));%L1=simplify(L11);
v=0.3;
s=0.5; 
f=1.93-3.07*s+14.53*s^2-25.11*s^3+25.8*s^4;
k=(E*I1)/(6*pi*h1*f*(1-v^2)); 
m=(E*I1+n*N)/N;
x=(E*I1+n*N)/(E*I0+n*N);
a11=(L0*a)-sin(L0*a);
a12=cos(L0*a)-1;
a13=(-tan(L1*l)*cos(L1*a))+sin(L1*a)-(a*(m*L1^3+L1)*((tan(L1*l)*sin(L1*l))+cos(L1*l)));
a14=1;
a21=(-L0*cos(L0*a))+L0+((E*I0+n*N)*L0^2*sin(L0*a));
a22=(-L0*cos(L0*a))-((E*I0+n*N)*L0^2*cos(L0*a));
a23=(k*L1*((tan(L1*l)*sin(L1*a))+cos(L1*a)))-((m*L1^3+L1)*((tan(L1*l)*sin(L1*l))+cos(L1*l)));
a24=0;
a31=L0^2*sin(L0*a);
a32=-L0^2*cos(L0*a);
a33=x*L1^2*((tan(L1*l)*cos(L1*a))-sin(L1*a));
a34=0;
a41=L0^3*cos(L0*a);
a42=L0^3*sin(L0*a);
a43=-x*L1^3*((tan(L1*l)*sin(L1*a))+cos(L1*a));
a44=0;
A=[a11 a12 a13 a14;a21 a22 a23 a24;a31 a32 a33 a34;a41 a42 a43 a44];
delta1=det(A);
delta2=subs(delta1);
delta=simplify(delta2);
% P=vpasolve(delta,P,1e3);
F=matlabFunction(delta);
P1=[];
for j = 0.01:10
    PP = fzero(F, j);
    P1 = [P1; PP];
end
P_cr = min(P1(P1 > 0)) 

It reports the following kind of errors,

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -2.5228 is 3.3961e-18-1.174e-23i.)
Check function or try again with a different starting value.

P_cr =

  0×1 empty double column vector

My question is that, how could i get positive real solutions. I shall be very thankful if somone can help me regarding this problem.

Answer 1

[Dana]

First of all, it looks like P is maybe supposed to be restricted to be a positive number, and when it's not, your function returns a complex number. However, fzero doesn't respect bounds on the search variable (P in your case): it looks for any value of P that zeros the function, regardless of whether P is postiive or negative. If it encounters a complex function value at some point, however (e.g., in your case if it evaluates F(P) for P<0), fzero doesn't know how to proceed, so it stops and reports an error. This is what you're seeing I believe.

There are various ways you can fix this problem, but the simplest to implement in your case is probably just to do a variable transformation like logP = log(P), P = exp(logP). For any real value of logP, P will be positive here, so you can use fzero to search among values of logP instead of directly for P. Here's how you might implement this:

% Replace the first line of your posted code with the following two:
syms logP
P = exp(logP);

. 
. % Then everything is the same until...
.

% replacing what you had, starting from the line P1=[];
logP1=[];
for j = log(0.01:10)
    logPP = fzero(F, j);
    logP1 = [logP1; logPP];
end
P1 = exp(logP1);
P_cr = min(P1(P1 > 0)) 

I think this should do what you want.

A couple of other comments/pointers:

• Do you really need to do this with symbolic variables? Couldn't you just write a regular MATLAB function called F.m that contains everything up to just before you start the fzero loop, but where logP is the argument passed into the function, and excluding the syms and the call to matlabFunction? Symbolic math is slow, so if you can avoid it you should.
• You have some crazy scaling on these values. For values of P between 0 and 10, the absolute value of F is on the order of 10^(-16) or smaller, but if you inspect F, these values are arrived at using some intermediate numbers that are on the order of like 10^(+15). That tells me that your F value is likely seriously compromised by numerical precision issues. Are you sure you haven't made a mistake somewhere in your code? It would be pretty unusual to encounter something like this in the real world.

[HINA]

Hello Dana,

Thank you very much for a detailed explaination to my problem.

About your comments: I am just a beginner at MATLAB and trying to learn it. I appreciate the suggestion provided in first comment and I 'll to implement it on my problem. The second comment is a very serious mistake made my me. Here I am looking for smallest nonzero positive root of P and as I understand i should recontruct my loop as follows:

logP1=[];
for j = log(0.01:10^15)
    logPP = fzero(F, j);
    logP1 = [logP1; logPP];
end
P1 = exp(logP1);
P_cr = min(P1(P1 > 0))

I will highly appreciate your further comments on this problem.

Answer 2

[Walter Roberson]

The only real root of your equation near 0 is 0. (There might be additional roots at large enough P.)

Near 0 the derivative is strictly negative for positive P.

Dana's explanation of why complex values are encountered is correct. In order to find the root at 0 fzero would have to bracket it from both sides but it cannot do that because of the complex values at negative P.