val = Empty sym: 0-by-1

Hello Guys!

I want to solve a series of dependent nonlinear equations in MATLAB. I have 31 equations and 31 unknown parameters. I entered these equations in MATLAB but when I use the 'solve' command, I get the following error in solved parameters:

val = Empty sym: 0-by-1

Does anyone know the reason for this problem? Is there a solution for it?

My code is:

clc 
clear 
close all

%DATA
di = 0.04; % m
alpha = 12.67; % degree
beta = 2.66; % degree
Ti = 291.65; % K
Pi = 90; % bar
Vi = 43.5; % m/s
k = 1.32;
R = 516.18; % j/kg.k
L = 0.2; %m
Mat=1;
Zi = 1;
Zt=1;
Zbs=1;
Zas=1;
Ze=1;


%inlet
Ai = pi*(di^2)/4; % m^2
Ci = sqrt(k*R*Zi*Ti); % Ci = 443.8; m/s
Mi = Vi/Ci; % Mi = 0.098; 
rhoi = (Pi*1e5)/(Zi*Ti*R); % rhoi = 60.17;  kg/m^3

%outlet
Pe = 70; %atm
de = 0.03; %m
Ae = pi*(de^2)/4; % m^2

%Equations:
syms Pt rhot Tt dt At Ct Vt xt Pbs rhobs Tbs dbs Abs Cbs Vbs Mabs xbs Pas rhoas Tas das Aas Cas Vas Maas xas rhoe Te Ce Ve Mae
eq1 = Ct == sqrt(k*R*Zt*Tt);
eq2 = (k*R*Zt*Tt) == (Vi^2+((2*k)/(k-1))*(R*Zi*Ti-R*Zt*Tt));
eq3 = ((Zi*Ti)/(Zt*Tt)) == (rhoi/rhot)^(k-1);
eq4 = (rhoi/rhot)^(k-1) == (Pi/Pt)^((k-1)/k);
eq5 = At == Ai*(rhoi/rhot)*(Vi/Vt);
eq6 = xt == (Ai^0.5 - At^0.5)/((pi^0.5)*tand(alpha)) ;
eq7 = At == pi*(dt^2)/4;
eq8 = Mat == Vt / Ct;

eq9 = Abs == At*(rhot/rhobs)*(Vt/Vbs);
eq10 = ((Zt*Tt)/(Zbs*Tbs)) == (rhot/rhobs)^(k-1);
eq11 = (rhot/rhobs)^(k-1) == (Pt/Pbs)^((k-1)/k);
eq12 = Cbs == sqrt(k*R*Zbs*Tbs);
eq13 = Vbs^2 == (Vt^2+((2*k)/(k-1))*(R*Zt*Tt-R*Zbs*Tbs));
eq14 = (Abs^0.5 - At^0.5)/((pi^0.5)*tand(beta)) == xbs;
eq15 = Mabs == Vbs / Cbs;
eq16 = Abs == pi*(dbs^2)/4;


eq17 = Ae == Aas*(rhoas/rhoe)*(Vas/Ve);
eq18 = ((Ze*Te)/(Zas*Tas)) == (rhoe/rhoas)^(k-1);
eq19 = (rhoe/rhoas)^(k-1) == (Pe/Pas)^((k-1)/k);
eq20 = Cas == sqrt(k*R*Zas*Tas);
eq21 = Ce == sqrt(k*R*Ze*Te);
eq22 = Vas^2 == (Ve^2+((2*k)/(k-1))*(R*Ze*Te-R*Zas*Tas));
eq23 = (Ae^0.5 - Aas^0.5)/((pi^0.5)*tand(beta)) == xas;
eq24 = Mae == Ve / Ce;
eq25 = Maas == Vas / Cas;
eq26 = Aas == pi*(das^2)/4;

eq27 = Pbs*(1+k*Mabs^2) == Pas*(1+k*Maas^2);
eq28 = Tbs*(1+((k-1)/2)*(Mabs^2)*Zbs) == Tas*(1+((k-1)/2)*(Maas^2)*Zas);
eq29 = ((Mabs)/(1+k*Mabs^2))*sqrt(1+((k-1)/2)*Mabs^2) == sqrt(Zbs/Zas)*(Maas)/(1+k*Maas^2)*sqrt(1+((k-1)/2)*Maas^2);
eq30 = Abs == Aas;
eq31 = L == xas+xbs+xt;



[Ptt,rhott,Ttt,dtt,Att,Ctt,Vtt,xtt,Pbst,rhobst,Tbst,dbst,Abst,Cbst,Vbst,Mabst,xbst,Past,rhoast,Tast,dast,Aast,Cast,Vast,Maast,xast,rhoet,Tet,Cet,Vet,Maet]=solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19,eq20,eq21,eq22,eq23,eq24,eq25,eq26,eq27,eq28,eq29,eq30,eq31],[Pt,rhot,Tt,dt,At,Ct,Vt,xt,Pbs,rhobs,Tbs,dbs,Abs,Cbs,Vbs,Mabs,xbs,Pas,rhoas,Tas,das,Aas,Cas,Vas,Maas,xas,rhoe,Te,Ce,Ve,Mae]);

[Steven Lord]

It is easy to write a system of equations that has no solution.

Consider as one example the two functions y = x^2 and y = x^2+1. Is there a pair of values for x and y that satify both functions simultaneously? You should be able to easily convince yourself that there is not.

What reason do you have for believing that your system of equations has a solution? Have you double-checked that they have been translated from the text or notes where you obtained them into code accurately? Sometimes something as small as one typo (a + instead of a minus or vice versa) can have a big impact on what and how the code does.

[M.S.R]

You're right! but the equations I wrote are related to the design of a nozzle, which were copied from a senior thesis and entered correctly. In addition, I have the answer of these equations according to the input values. In fact, I am trying to write a code to solve these with MATLAB, so that I can easily get the desired results by changing the desired values.

[Torsten]

In addition, I have the answer of these equations according to the input values.

Analytically or numerical values ?

Are you sure the equations permit analytical expressions for the unknowns in terms of the input parameters ?

You must go through the equations step by step as I did in your previous post and see which unknowns can be eliminated in the next step.

"Solve" is a mathematical tool - not a creative human brain.

[John D'Errico]

Just because you think a solution should exist does not mean it can be found. In fact, it is trivial to write an equation that has NO analytically solvable solution, even though one does exist. For example:

syms a b c d e real
syms x
P5 = x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e == 0

P5 = 

Since the polynomial is a general polynomial of odd degree with real coefficients, it is mathematically provable that at least ONE real solution MUST exist.

At the same time, Abel-Ruffini provided a proof (long before any of us were born) that no algebraic solution can be found for exactly such a problem. So solve would NEVER return a result there, even for such a simply written problem.

solve can sometimes produce a solution for problems with completely numerical coefficients, but only because it uses a numerical solver there when necessary. And for large nonlinear problems, numerical solvers can often be subject to starting value issues.

As it turns out, your problem is wildly more complex than that simple polynomial. But solve was unable to solve your problem. No huge surprise there. We all sometimes want things that are not possible to achieve.

[M.S.R]

All the points you mentioned are correct, but what I am sure of is that there is an analytical solution to these equations.

I also know that these equations have already been coded and solved in MATLAB by someone else.

For the first 8 equations that relate to the nozzle inlet to its throat, I found the appropriate analytical solution code that perfectly matches the answers I have. (question of my previous post)

clc 
clear
close all

%DATA
di = 0.04; % m
alpha = 12.67; % degree
beta = 2.66; % degree
Ti = 291.65; % K
Pi = 90; % bar
Vi = 43.5; % m/s
k = 1.32;
R = 516.18; % j/kg.k

Zi = 1;
Zt=1;
Mat=1;

%inlet
Ai = pi*(di^2)/4; % m^2
Ci = sqrt(k*R*Zi*Ti); % Ci = 443.8; m/s
Mi = Vi/Ci; % Mi = 0.098; 
rhoi = (Pi*1e5)/(Zi*Ti*R); % rhoi = 60.17;  kg/m^3


%Equations:
syms Pt rhot Tt xt dt Vt At Ct
eq1 = sqrt(k*R*Zt*Tt) == Ct;
eq2 = (k*R*Zt*Tt) == (Vi^2+((2*k)/(k-1))*(Zi*R*Ti-Zt*Tt*R));
eq3 = ((Zi*Ti)/(Zt*Tt)) == (rhoi/rhot)^(k-1);
eq4 = (rhoi/rhot)^(k-1) == (Pi/Pt)^((k-1)/k);
eq5 = Ai*(rhoi/rhot)*(Vi/Vt) == At;
eq6 = (Ai^0.5 - At^0.5)/((pi^0.5)*tand(alpha)) == xt;
eq7 = At == pi*(dt^2)/4;
eq8 = Mat == Ct / Vt;


[PPt,rrhot,TTt,xxt,ddt,VVt,AAt,Ctt]=solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8],[Pt,rhot,Tt,xt,dt,Vt,At,Ct]);

From equation 9 onward, all equations are interdependent and therefore Salo's command cannot solve them easily, and this is exactly what I am looking for!

Perhaps simplifying the equations can help in this regard, and perhaps each equation should be solved separately and the obtained parameter replaced in another equation, but the point is that it cannot be done easily and the dependence of the equations prevents it !

Answer 1

[Walter Roberson]

%DATA
di = 0.04; % m
alpha = 12.67; % degree
beta = 2.66; % degree
Ti = 291.65; % K
Pi = 90; % bar
Vi = 43.5; % m/s
k = 1.32;
R = 516.18; % j/kg.k
L = 0.2; %m
Mat=1;
Zi = 1;
Zt=1;
Zbs=1;
Zas=1;
Ze=1;


%inlet
Ai = pi*(di^2)/4; % m^2
Ci = sqrt(k*R*Zi*Ti); % Ci = 443.8; m/s
Mi = Vi/Ci; % Mi = 0.098; 
rhoi = (Pi*1e5)/(Zi*Ti*R); % rhoi = 60.17;  kg/m^3

%outlet
Pe = 70; %atm
de = 0.03; %m
Ae = pi*(de^2)/4; % m^2

%Equations:
syms Pt rhot Tt dt At Ct Vt xt Pbs rhobs Tbs dbs Abs Cbs Vbs Mabs xbs Pas rhoas Tas das Aas Cas Vas Maas xas rhoe Te Ce Ve Mae
eq1 = Ct == sqrt(k*R*Zt*Tt);
eq2 = (k*R*Zt*Tt) == (Vi^2+((2*k)/(k-1))*(R*Zi*Ti-R*Zt*Tt));
eq3 = ((Zi*Ti)/(Zt*Tt)) == (rhoi/rhot)^(k-1);
eq4 = (rhoi/rhot)^(k-1) == (Pi/Pt)^((k-1)/k);
eq5 = At == Ai*(rhoi/rhot)*(Vi/Vt);
eq6 = xt == (Ai^0.5 - At^0.5)/((pi^0.5)*tand(alpha)) ;
eq7 = At == pi*(dt^2)/4;
eq8 = Mat == Vt / Ct;

eq9 = Abs == At*(rhot/rhobs)*(Vt/Vbs);
eq10 = ((Zt*Tt)/(Zbs*Tbs)) == (rhot/rhobs)^(k-1);
eq11 = (rhot/rhobs)^(k-1) == (Pt/Pbs)^((k-1)/k);
eq12 = Cbs == sqrt(k*R*Zbs*Tbs);
eq13 = Vbs^2 == (Vt^2+((2*k)/(k-1))*(R*Zt*Tt-R*Zbs*Tbs));
eq14 = (Abs^0.5 - At^0.5)/((pi^0.5)*tand(beta)) == xbs;
eq15 = Mabs == Vbs / Cbs;
eq16 = Abs == pi*(dbs^2)/4;


eq17 = Ae == Aas*(rhoas/rhoe)*(Vas/Ve);
eq18 = ((Ze*Te)/(Zas*Tas)) == (rhoe/rhoas)^(k-1);
eq19 = (rhoe/rhoas)^(k-1) == (Pe/Pas)^((k-1)/k);
eq20 = Cas == sqrt(k*R*Zas*Tas);
eq21 = Ce == sqrt(k*R*Ze*Te);
eq22 = Vas^2 == (Ve^2+((2*k)/(k-1))*(R*Ze*Te-R*Zas*Tas));
eq23 = (Ae^0.5 - Aas^0.5)/((pi^0.5)*tand(beta)) == xas;
eq24 = Mae == Ve / Ce;
eq25 = Maas == Vas / Cas;
eq26 = Aas == pi*(das^2)/4;

eq27 = Pbs*(1+k*Mabs^2) == Pas*(1+k*Maas^2);
eq28 = Tbs*(1+((k-1)/2)*(Mabs^2)*Zbs) == Tas*(1+((k-1)/2)*(Maas^2)*Zas);
eq29 = ((Mabs)/(1+k*Mabs^2))*sqrt(1+((k-1)/2)*Mabs^2) == sqrt(Zbs/Zas)*(Maas)/(1+k*Maas^2)*sqrt(1+((k-1)/2)*Maas^2);
eq30 = Abs == Aas;
eq31 = L == xas+xbs+xt;


eqns = [eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19,eq20,eq21,eq22,eq23,eq24,eq25,eq26,eq27,eq28,eq29,eq30,eq31];

current = eqns;
part1to8 = solve(eqns(1:8));
current2 = subs(eqns(9:end), part1to8);
vars1to8 = sym(fieldnames(part1to8));
Nv1to8 = length(vars1to8);

part9to12 = solve(current2(1,1:4));
vars9to12 = sym(fieldnames(part9to12));
Nv9to12 = length(vars9to12);

current3 = subs(current2(1,5:end), part9to12);
part13to15 = solve(current3(1:2), [Abs, xbs, Mabs]);
vars13to15 = sym(fieldnames(part13to15));
Nv13to15 = length(vars13to15);

current4 = subs(current3(1,3:end), part13to15);

But at this point you get stuck. current4(1,1) is only in terms of the variable Cbs. The left hand side is 0; the right hand side is a ratio of two polynomials, with a discontinuity at Cbs == 0, and it does not cross 0, so there is no solution.

If you follow through the code you will see that at each place I generate a new "current" variable, we are dealing with multiple solutions. In the above code, in each case I select the "first" of the solutions (the second solution generates rows 2 and onwards). So it is potentially possible that if you were to select different branches at those points that you might be able to get further.