Matlab fsolve function not returning correct answer.

5 次查看（过去 30 天）
Casey
Casey 2016-9-12
评论： Alex Sha 2020-2-2
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
  1 个评论
Alex Sha
Alex Sha 2020-2-2
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2016-9-12
编辑：John D'Errico 2016-9-12
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.
  4 个评论
显示 2更早的评论
monty singh
monty singh 2017-1-20
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
Walter Roberson 2017-1-20
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by