Solving iteratively a system of equations by vpasolve while there are "if statements" in the system
5 次查看(过去 30 天)
显示 更早的评论
Hi,
I have a system of three unknown variables that I would like to solve by vpasolve iteratively, but there are "if statements" in the Eq2 of the system and vpasolve does not accept them. I was wondering if there is any solution for that or there are any other solver functions in MATLAB that works with "if statements". I put my MATLAB code here
Thank
g = 9.818; rhoP = 2700; rho = 1060; mu = 0.0063; dp = 88e-4; Cv = 11; Fc = 0.0999
syms Nre Cd v
Eq1 = Nre == dp * rho * v * Cv / mu
if Nre <= 1
Eq2 = Cd == 24 ./ Nre;
elseif Nre < 1000
Eq2 = Cd == (1 + 0.15 * Nre .^0.687) * 24 ./ Nre;
else Eq2 = Cd == 0.44;
end
Eq3 = v == sqrt(4 / 3 * dp * g * Fc * (rhoP / rho - 1) / Cd);
[Nre,Cd,v] = vpasolve(Eq1,Eq2,Eq3,Nre,Cd,v)
0 个评论
采纳的回答
Walter Roberson
2019-3-1
If you run the attached code you will find that v is 0 for the solution. If you go through the three piecewise conditions, you will find that v = 0 generates Nre = 0 which passes Nre <= 1 test. If you take the three piecewise parts individually and solve for v for each of them, the first one works with v = 0. The second one gives a v > 2 which fails the Nre < 1000 test. The third one gives the constant v = 11/25 but that branch cannot be reached because when v = 11/25 that would be caught by the Nre <= 1 branch. So there are no alternative solutions for v, only 0.
However... when v = 0 becaue Nre = 0, and you define Cd = 24/Nre that would be division by 0, giving Cd = infinity. Then in your definition of Eq3 on the right hand side, you have sqrt() of a positive quantity divided by that infinity, giving 0. This is consistent, 0 == 0, but it is icky.
Final solution: Nre = 0, v = 0, Cd = infinity
3 个评论
Walter Roberson
2019-3-1
Your revised equations gain a second solution beyond 0: v = 16000885^(1/2)/19875 does not match the first two conditions, so the third conditional case is consistent. This is Cd = 0.44 and Nre = (3872*160008855^(1/2))/4725
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!