Matlab fsolve function not returning correct answer.

Hello,

I am trying to solve a system of 8 non-linear equations with 8 unknowns using the Matlab fsolve function and the below function file. The point of the exercise is to find the solutions for V=0,1,2,...,10. The resultant data for x(6) is supposed to be linear up to a point, but then level off abruptly. Now with x = [1 1 1 1 1 1 1 1] and executing "x=fsolve(@myfun,x0)", I am getting the correct values for V=0-6, but at V=7, I can't seem to get Matlab to spit out the correct answer. I should be getting something above 23, but I keep getting an imaginary number like 5+0.26i. Can you please tell me why this is happening and how I fix it? Thank you.

function G = myfun(x)
V=7;
cb=100;
omega=94.25;
T=298;
H=0.1;
Dp=0.72*10^(-9);
Dm=1.065*10^(-9);
iO=10;
R=8.314;
F=96487;
dA=0.06;
dC=0.04;
alphaA=0.84;
alphaC=1.16;
nu=0.89*10^(-6);
gamma=0.42;
zp=2;
zm=-2;
D=(Dp*Dm*(zp-zm))/(zp*Dp-zm*Dm);
kmta = (D/dC)*0.0791*((omega*dC^3)/(2*nu*dA))^0.7*(nu/D)^0.356;
kmtc = (D/dC)*0.0791*((omega*dC^2)/(2*nu))^0.7*(nu/D)^0.356;
tp = (zp*Dp)/(zp*Dp-zm*Dm);
kappa = (F^2/(R*T))*((zp^2*Dp)+(zm^2*Dm))*cb;

G(1) = -V+x(1)+x(2)+x(3)-x(4)-x(5);
G(2) = -x(6)/(pi*H*dA)+iO*(x(7)/cb)^gamma*(exp(alphaA*F*x(2)/(R*T))-exp(-alphaC*F*x(2)/(R*T)));
G(3) = x(6)/(pi*H*dC)+iO*(x(8)/cb)^gamma*(exp(alphaA*F*x(4)/(R*T))-exp(-alphaC*F*x(4)/(R*T)));
G(4) = -x(6)/(zp*F)+(pi*dA*H)*(-kmta*(cb-x(7))/(1-tp));
G(5) = x(6)/(zp*F)+(pi*dC*H)*(-kmtc*(cb-x(8))/(1-tp));
G(6) = -x(3)+R*T/F*log(x(7)/cb)+R*T/F*tp*(1-x(7)/cb);
G(7) = -x(5)+R*T/F*log(x(8)/cb)+R*T/F*tp*(1-x(8)/cb);
G(8) = -x(1)+(x(6)/(2*pi*H*kappa)*log(dA/dC));
end

[Alex Sha]

An approximate real-value solution:

x1: 6.43175857461548
x2: 0.142388505395668
x3: 0.0075073193074162
x4: -0.201073218358661
x5: -0.217272383006374
x6: 26.7402614677563
x7: 202.768554734273
x8: 0.0141275553155712

Fevl:
6.83599149509106E-10
-6.08679329161532E-10
1.07775122160092E-10
-8.00298647320652E-8
1.92619586341643E-5
7.62215120880816E-10
-4.01468596214483E-10
1.78754788748847E-11

Answer 1

[John D'Errico]

You need to understand that the equations in here have some regions where complex numbers will result from some of those computations, depending on the parameters. Once a complex result happens, fsolve gets confused.

How could that possibly happen?

What is gamma? (gamma is a terrible name to use for a variable, by the way, because one day you will want to use the gamma function, a terribly useful tool. But if you have variables named gamma, then forget about being able to use the FUNCTION named gamma.

gamma=0.42;

So gamma is 0.42. How are you using gamma?

(x(7)/cb)^gamma

As an exponent. What is a negative number raised to the 0.42 power?

(-1) ^ 0.42
ans =
          0.248689887164855 +     0.968583161128631i

Similarly, you have places where you compute the log of a number.

log(x(8)/cb)

Logs also don't like negative numbers.

So if your parameters ever go negative, fsolve will have problems, as then complex numbers will propagate through the computations.

Fsolve does not offer constraints. So even if one trial ends up with a negative value for one of the elements, fsolve will have problems.

The classic and simple solution is to use better starting values.

[Casey]

Thank you for your answer. Unfortunately, that solution does not apply to this case. x(7) and x(8) represent ionic concentrations, and so their values will never be negative, because you can't never a negative concentration of something. If the solution to this system gave negative values, then I would have setup the system correctly, which I haven't, because I've checked the system for lower values of V.

However, could you please elaborate on what you mean when you say "use better starting values"?

Thank you

[Walter Roberson]

It is not a matter of the system "giving" negative values: it is a matter that when fsolve() is trying to calculate the roots, it probes at negative values to try to predict where the root is. Then, having reached a complex result, it gets confused.

One of the options you can specify to fsolve is 'FunValCheck', 'on', which will cause the computation to abort if it encounters a complex result, or an inf, or a NaN.

[monty singh]

Then, how to solve this issue? Same problem is happening with me, its not possible to have negative or imaginary solution to the system. But fsolve gives imaginory solution and if we give the fulvalcheck condition it terminates in between. Error: Error using lsqfcnchk/checkfun (line 135) User function 'methane20/abc' returned a complex value when evaluated; FSOLVE cannot continue.

[Walter Roberson]

You need switch to a bounded solver or you need arrange for a noncomplex value to be returned if the probe values are outside of the useful area.

If you have the symbolic toolbox you could consider vpasolve as it accepts bounds if I recall correctly.