Possibly spurious solutions - Matlab blocked with no answers

Good day,

I am trying to solve a set of equations for a problem of gas release using the Abel Noble eq.

Here the script

----------------------------------------------------------------------------------------

clear all;
close all;
clc;

R=8.314;
M=2.016;
T1=288;
b=2.65/M/100000;
A=33.066178;
B=-11.363417;
C=11.432816;
D=-2.772874;
E=-0.158558;
cp=(A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2)/M*1000;
cv=(((A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2))-(R/M))/M*1000;
gamma=cp/cv;

RH2=R/M;%gas constant for H2
K=2; %reduction area guess
F=0.02; %friction factor guess
p1=3.5*10000000; %tank pressure
p4=101325; %ambient pressure
d3=0.0095; %orifice diameter

syms  p2 p3 u2 u3 u4 T2 T3 T4 ro2 ro3 ro4 d4
eqn1=p2-p1+(ro2*u2^2*(K/4+1))==0;
eqn2=cp*T2+(K+1)*u2^2/2==cp*T1;
eqn3=ro2==b*p2/(p2+(RH2*T2));
eqn4=p3-p2+(ro2*u2^2*(F/4-1))+(ro3*u3^2*(F/4+1))==0;
eqn5=cp*T2+u2^2/2==(cp*T3)+((F/4+1)*u3^2/2);
eqn6=ro3==p3/((b*p3)+(RH2*T3));
eqn7=ro2*u2==ro3*u3;
eqn8=u3==(gamma*RH2*T3)^0.5/(1-(b*ro3));
eqn9=u4==(gamma*RH2*T4)^0.5;
eqn10=ro4==p4/(RH2*T4);
eqn11=(cp*T3)+(u3^2/2)==(cp*T4)+(u4^2/2);
eqn12=ro3*u3*d3^2==ro4*u4*d4^2;

eqns=[eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12];
vars=[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4];
[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4]=solve(eqns,vars);

--------------------------------------------------------------------------------------------

Then the RUNNING MATLAB i only received the message: Possibly spurious solutions

with nothing else!!!!

Answer 1

[John D'Errico]

You are attempting to solve for an analytical solution of 12 simultaneous nonlinear equations, in 12 unknowns. This, even if it is possible to solve, will likely take an immense amount of time.

But what happens is solve does some work, and sees that the solution approach will very possibly introduce a solution that is not a valid solution. That often happens when you have powers of a variable or expression, since in order to resolve things, it may need to raise an expression to a power, which can introduce apparent solutions that are not really solutions.

The problem is, solve apparently saw this was going to happen, so it wrote out the warning message. But then it got bogged down. I tested your code, and it does exactly that. You might try the use of vpasolve, which I did on my computer for your problem, but it returns:

vpasolve(eqns,vars)
ans = 
struct with fields:

p2: [0×1 sym]
p3: [0×1 sym]
u2: [0×1 sym]
u3: [0×1 sym]
u4: [0×1 sym]
T2: [0×1 sym]
T3: [0×1 sym]
T4: [0×1 sym]
ro2: [0×1 sym]
ro3: [0×1 sym]
ro4: [0×1 sym]
d4: [0×1 sym]

So vpasolve found no solution, but I gave it no starting values since I had not a clue what to choose. Without any inteligent starting values, it is possible that vpasolve will start in a poor place and fail to converge. Or it may mean that no solution exists. You may have an error in your equations. I cannot know.

[davide papurello]

Really thanks, I will try to add some initial guess. Are you suggesting to force some values?

[John D'Errico]

No. Don't force any variable to some value, as that will merely insure that no solution is found at all.

The metaphor I often use for any optimization tool is that of a blind person, tasked with finding the point of lowest elevation on earth. Given only a cane and an altimeter, can they walk to the solution? Ok, give them scuba gear too.

The point is, you need to provide intelligent starting values, else you get complete crapola much of the time.

You might decide to try a tool like lsqnonlin. It is sometimes possible that it can find a near solution, even when no exact solution can be found. Even there though, you need to provide starting values. Good starting values are better than bad ones of course.

Answer 2

[Walter Roberson]

Turns out there are some solutions

In the below, fullsubs is not actually valid for all entries -- only for entries [5 8 13 16 25 27 30 32 41 43 46 48]

Note that it is a category error to use floating point coefficients with solve().

Floating point coefficients are inherently uncertain, expressing numbers that are "somewhere in the range" of roughly 1 part in 2^53. Numbers expressed to 6 or so floating point digits by convention express the entire range of numbers that round to the number that is 6 significant digits. -1.1363417 by convention expresses -113634175/10^7 to -113634165/10^7.

When you use floating point coefficients, you express uncertainty about what the number actually is. And uncertainty about what the number actually is, is antagonistic to the purpose of solve() -- which is to find indefinitely precise solutions whenever possible.

If you must use floating point coefficients, then don't use solve() -- use vpasolve()

Q = @(v) sym(v);

R = Q(8314)/Q(10)^3;
M = Q(2016)/Q(10)^3;
T1 = Q(288);
b = Q(265)/Q(10)^3/M/100000;
A = Q(33066178)/Q(10)^6;
B = Q(-11363417)/Q(10)^6;
C= Q(11432816)/Q(10)^6;
D = Q(-2772874)/Q(10)^6;
E= Q(-0158558)/Q(10)^6;
cp = (A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2)/M*1000;
cv = (((A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2))-(R/M))/M*1000;
gamma = cp/cv;

RH2 = R/M;%gas constant for H2
K = Q(2); %reduction area guess
F = Q(0.02); %friction factor guess
p1 = Q(3.5)*10000000; %tank pressure
p4 = Q(101325); %ambient pressure
d3 = Q(95)/Q(10)^4; %orifice diameter

syms  p2 p3 u2 u3 u4 T2 T3 T4 ro2 ro3 ro4 d4
eqn1=p2-p1+(ro2*u2^2*(K/4+1))==0;
eqn2=cp*T2+(K+1)*u2^2/2==cp*T1;
eqn3=ro2==b*p2/(p2+(RH2*T2));
eqn4=p3-p2+(ro2*u2^2*(F/4-1))+(ro3*u3^2*(F/4+1))==0;
eqn5=cp*T2+u2^2/2==(cp*T3)+((F/4+1)*u3^2/2);
eqn6=ro3==p3/((b*p3)+(RH2*T3));
eqn7=ro2*u2==ro3*u3;
eqn8=u3==(gamma*RH2*T3)^0.5/(1-(b*ro3));
eqn9=u4==(gamma*RH2*T4)^0.5;
eqn10=ro4==p4/(RH2*T4);
eqn11=(cp*T3)+(u3^2/2)==(cp*T4)+(u4^2/2);
eqn12=ro3*u3*d3^2==ro4*u4*d4^2;

eqns=[eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12].'

eqns = 

vars=[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4];
%[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4]=solve(eqns,vars)
sol1 = solve(eqns([1 2 3 4 5 6 8 9 10 11 12]), [u2 T2 ro2 p3 T3 ro3 u3 u4 ro4 T4 d4])

Warning: Possibly spurious solutions.

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

sol1 = struct with fields:
     u2: [48×1 sym]
     T2: [48×1 sym]
    ro2: [48×1 sym]
     p3: [48×1 sym]
     T3: [48×1 sym]
    ro3: [48×1 sym]
     u3: [48×1 sym]
     u4: [48×1 sym]
    ro4: [48×1 sym]
     T4: [48×1 sym]
     d4: [48×1 sym]

eqns2 = subs(eqns(7), sol1);

p2sol = arrayfun(@vpasolve, eqns2(5), 'uniform', 0)

p2sol = 1×1 cell array
    {[-1185.5858461081479176182578584473]}

backsubs = subs(sol1, p2, p2sol);
backsubs.p2 = repmat(p2sol, 48, 1);
fullsubs = subs(vars, backsubs)

fullsubs = 

[davide papurello]

Thanks for the answer, but I can not get the point of the last two lines. Is it reachable the final solution?

Thanks

[Walter Roberson]

p2sol = arrayfun(@vpasolve, eqns2(5), 'uniform', 0)

It happens that p2sol is the same for all of the equations that can be solved at all.

But to get the "real" answer, you need to

p2solutions = arrayfun(@vpasolve, eqns2, 'uniform', 0);

this will give you a mix of empty (no solution) and -1185.5858461081479176182578584473 (solution). Where you get empty, you have to discard the rows of fullsubs.

[Walter Roberson]

Really you should check all of the return conditions

Q = @(v) sym(v);

R = Q(8314)/Q(10)^3;
M = Q(2016)/Q(10)^3;
T1 = Q(288);
b = Q(265)/Q(10)^3/M/100000;
A = Q(33066178)/Q(10)^6;
B = Q(-11363417)/Q(10)^6;
C= Q(11432816)/Q(10)^6;
D = Q(-2772874)/Q(10)^6;
E= Q(-0158558)/Q(10)^6;
cp = (A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2)/M*1000;
cv = (((A+B*(T1/1000)+C*(T1/1000)^2+D*(T1/1000)^3+E/(T1/1000)^2))-(R/M))/M*1000;
gamma = cp/cv;

RH2 = R/M;%gas constant for H2
K = Q(2); %reduction area guess
F = Q(0.02); %friction factor guess
p1 = Q(3.5)*10000000; %tank pressure
p4 = Q(101325); %ambient pressure
d3 = Q(95)/Q(10)^4; %orifice diameter

syms  p2 p3 u2 u3 u4 T2 T3 T4 ro2 ro3 ro4 d4
eqn1=p2-p1+(ro2*u2^2*(K/4+1))==0;
eqn2=cp*T2+(K+1)*u2^2/2==cp*T1;
eqn3=ro2==b*p2/(p2+(RH2*T2));
eqn4=p3-p2+(ro2*u2^2*(F/4-1))+(ro3*u3^2*(F/4+1))==0;
eqn5=cp*T2+u2^2/2==(cp*T3)+((F/4+1)*u3^2/2);
eqn6=ro3==p3/((b*p3)+(RH2*T3));
eqn7=ro2*u2==ro3*u3;
eqn8=u3==(gamma*RH2*T3)^0.5/(1-(b*ro3));
eqn9=u4==(gamma*RH2*T4)^0.5;
eqn10=ro4==p4/(RH2*T4);
eqn11=(cp*T3)+(u3^2/2)==(cp*T4)+(u4^2/2);
eqn12=ro3*u3*d3^2==ro4*u4*d4^2;

eqns=[eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12].'

eqns = 

vars=[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4];
%[p2,p3,u2,u3,u4,T2,T3,T4,ro2,ro3,ro4,d4]=solve(eqns,vars)
sol1 = solve(eqns([1 2 3 4 5 6 8 9 10 11 12]), [u2 T2 ro2 p3 T3 ro3 u3 u4 ro4 T4 d4], 'returnconditions', true)

Warning: Possibly spurious solutions.

sol1 = struct with fields:
            u2: [48×1 sym]
            T2: [48×1 sym]
           ro2: [48×1 sym]
            p3: [48×1 sym]
            T3: [48×1 sym]
           ro3: [48×1 sym]
            u3: [48×1 sym]
            u4: [48×1 sym]
           ro4: [48×1 sym]
            T4: [48×1 sym]
            d4: [48×1 sym]
    parameters: z
    conditions: [48×1 sym]

sol1.conditions

ans = 

Answer 3

[Alex Sha]

There is a approximate solution for original equations：

p2: 35001350.3279785
p3: 35113789.3799513
u2: -8.90075607998193
u3: -1115.00289883981
u4: 37.2181947223662
t2: 287.991671106106
t3: 244.208945707645
t4: 287.728066504045
ro2: -11.3630314304234
ro3: -0.0907078996826766
ro4: 85.3915491921384
d4: -0.00169472452853458

Fevel:
2.18617515201913E-8
9.77888703346252E-9
-5.38904565416942E-9
6.85395207256079E-9
-4.19095158576965E-9
-0.00270597877844363
-2.48826381721301E-9
5.77870196138974E-9
7.08126890458516E-11
-1.25750432289351E-9
1.2572854757309E-8
-2.28636101197444E-9

If making some reformulation like (Turning division into multiplication):

ro2=b*p2/(p2+(RH2*T2)) -> (p2+(RH2*T2))*ro2=b*p2;

ro3=p3/((b*p3)+(RH2*T3)); -> ((b*p3)+(RH2*T3))*ro3=p3;

u3=(gamma1*RH2*T3)^0.5/(1-(b*ro3)); -> (1-(b*ro3))*u3=(gamma1*RH2*T3)^0.5;

ro4=p4/(RH2*T4); -> (RH2*T4)*ro4=p4;

Then the result will be a little better:

p2: 35008358.7526556
p3: 35725949.1997961
u2: -22.145118119983
u3: -2859.44111961317
u4: 37.1376872183248
t2: 287.948442775181
t3: 7.70387476147959E-22
t4: 286.484630895113
ro2: -11.3630315600834
ro3: -0.0880016987847933
ro4: 85.762176031689
d4: -0.00267026541160297

Feval:
-2.29192664846778E-8
3.25962901115417E-9
1.19209289550781E-7
3.59723344445229E-8
-6.98491930961609E-9
-1.49011611938477E-8
-1.90871674021764E-9
-6.2170057916112E-11
-2.00515870574236E-9
-2.95403879135847E-9
9.77888703346252E-9
-4.28637500840545E-9