val = Empty sym: 0-by-1
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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]);
5 commentaires
John D'Errico
le 30 Nov 2022
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
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.
Réponses (1)
Walter Roberson
le 30 Nov 2022
%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.
0 commentaires
Voir également
Catégories
En savoir plus sur Linear Least Squares dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!