solving nonlinear equation by fsolve

Hey guys,

I am having some problems with my function code. The main idea is to get parameter SI (main unknown) in which the equation Q_cal-Q=0. Can somebody help me?

Thanks a lot.

P=1.94;
Q=1.09;
P5=1.08;
fc=0;
lambda=0.2;
Ts=24;

[SI]=singhandyu(P,Q,P5,lambda,Ts,fc);


    function [SI] =singhandyu(P,Q,P5,lambda,Ts,fc)

    Fc=fc.*Ts; 

f=@(SI)((P5-0.2*SI)*SI)./(P5+0.8*SI);
M=@(SI)max(f(SI),0);
S=@(SI)(SI-M(SI));

Ia=@(SI)lambda.*S(SI);
Q_cal=@(SI)((P-Ia(SI)-Fc).*(P-Ia(SI)-Fc+M(SI)))./(P-Ia(SI)-Fc+M(SI)+S(SI));

H=@(SI)Q_cal(SI)-Q;
S0=0;
SI_sol=fsolve(H,S0)
end

Answer 1

[Star Strider]

You need to return ‘SI_sol’ from ‘singhandyu’, so declare it as:

function [SI_sol] =singhandyu(P,Q,P5,lambda,Ts,fc)

Also, 0 is usually not an acceptable initial parameter estimate. If you set ‘S0’ to 1 instead:

S0=1;
SI_sol=fsolve(H,S0)

you get:

SI =
     1.7414

Answer 2

[John D'Errico]

Using max inside an objective that will then be applied to fsolve is a bad idea. It creates a function that is at best non-differentiable, and may be discontinuous.

One assumption implicit in the use of fsolve? That your objective be both continuous AND differentiable.

I'm not sure what your problem is in that code. One issue is you seem to be wildly overusing functions. You define functions just to do some trivial operation, like multiply two numbers. That makes your code difficult to read and debug.

Anyway, this function is a function of one variable. First, let me get rid of all of that function crap you had, just writing a single function to evaluate it.

P=1.94;
Q=1.09;
P5=1.08;
fc=0;
lambda=0.2;
Ts=24;

function H = obj(SI,P,Q,P5,lambda,Ts,fc)
  Fc = fc.*Ts;
  f = (P5-0.2*SI).*SI./(P5+0.8*SI);
  M = max(f,0);
  S = SI - M;
  Ia = lambda.*S;
  Q_cal = ((P-Ia-Fc).*(P-Ia-Fc+M))./(P-Ia-Fc+M+S);
  H = Q_cal - Q;
end

See that I've care fully used .* where needed, something you missed in at least one spot. As well, eliminating the function crap,

Next, TEST IT!!!!!!!!!

obj(0:3,P,Q,P5,lambda,Ts,fc)
ans =
         0.85      0.39942     -0.12131      -0.4992

So, when I pass in a vector of inputs, I get out a vector of the same size and shape. I'm not worried yet about problems with the max. Anyway, I won't be using fsolve, but fzero. fzero is a little more tolerant of issues.

Before I even try that however, PLOT IT!!!!!!! A golden rule that is incredibly valuable is to PLOT EVERYTHING! Don't just throw something into an optimizer without doing so. That is pleading for garbage out the end.

fplot(@(SI) obj(SI,P,Q,P5,lambda,Ts,fc),[0,100])
grid on
yline(0);
xlabel SI

That looks to be quite reasonable. NOW you can use fzero. But which root would you look for? There are at least two of them. It looks like this relation seems to be asymptotically linear as SI grows large, so I expect only two positive roots.

[x1,fval] = fzero(@(SI) obj(SI,P,Q,P5,lambda,Ts,fc),[0,10])
x1 =
       1.7414
fval =
   2.2204e-16

[x2,fval] = fzero(@(SI) obj(SI,P,Q,P5,lambda,Ts,fc),[10,1000])
x2 =
       40.174
fval =
  -2.2204e-16

In the end, the max call inside there seems not to have been an issue. But again, fzero is more tolerant of stuff like that.