Solving iteratively a system of equations by vpasolve while there are "if statements" in the system

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)

Answer 1

[Walter Roberson]

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

[Maryam S]

Thanks a lot Walter or your answer, you are right. But, I apply exactly the literature formulas for this problem (see the attached image) and it works well in Excel circular reference formulas, but not in MATLAB!

[Walter Roberson]

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

[Maryam S]

Thank you so much! I learnt a lot!