Solve 3 order 2 variables function

Hi there,

I'm tring to solve folowing two equations, the two dependent variables are u2 and beta. My result doesn't make sense to me, don't know how to fix that. My file is attached.

In my results, I got symble z, not sure what is that stand for. Can someone help me solve this problem?

Thanks for the help.

-Hanrui

Here is my code:

clc
clear
syms x y
%Isentropic no fraction loss
P_2=10^5;
rho_g=1.225;
rho_f=997;
A_1g=0.1^2*pi/4;
A_2=0.5^2*pi/4;
A_1f=A_2-A_1g;
b=A_1g/A_2;
cp=1;
T1=323;
T2=293;
gamma=1.4;
R=287;
q=10000;
a=0.3;
B=5;
theta=1;
L=0.5;
M=1;
R=287;
u_1g=M*sqrt(gamma*R*T1);
p_tg=5*10^6;

% x and y are two unknowns that need to be solved 
eq1=p_tg*(1+M*(gamma-1)/2)^((gamma-1)/gamma)-P_2-rho_f*y^2*(1-a)*(1-(1-a)/(1-x))-rho_g*(y^2*a-u_1g^2*x)==0;
eq2=y^3*(1-a)*rho_f*(1-((1-a)/(1-x))^2)/2+y^3*a*rho_g*(1-(a/x)^2)/2+a*y*rho_g*cp*(T2-T1)==0;
%solve above equation
sol=solve(eq1,eq2,x,y);
%output u_2 and beta
beta=sol.x
u_2=sol.y

Results:

beta =

(767813726674340823155714050737263296411427490095261963310899673231481657405553608578345025458736310901807929501764900692329928923713489397556246227961511936*root(z^10 - (76245256816524517094787721664575*z^8)/5352866887287026539823104.......

.

.

.

.

.

u_2 =

root(z^10 - (76245256816524517094787721664575*z^8)/5352866887287026539823104 + (932018315641293413021053976330724993678573765595931929*z^6)/1232904261651835201667968301579019113988096 - (15149901962504148690356548006017832826970719790718281499456798538626475*z^4)/1334411437689202100645463103523731168093117650233720832 + (21560789318824412760547720965355965940321374664217347794108452577399262109977756586375*z^2)/641900390239748061475957626285817825407851726928630145741397950464 - 644430211103063219451233473971997200503671748328219055060769915199066015625/607305472086107978249312985781466984945489988373035614208, z, 1)

root(z^10 - (76245256816524517094787721664575*z^8)/5352866887287026539823104 + (932018315641293413021053976330724993678573765595931929*z^6)/1232904261651835201667968301579019113988096 - (15149901962504148690356548006017832826970719790718281499456798538626475*z^4)/1334411437689202100645463103523731168093117650233720832 + (21560789318824412760547720965355965940321374664217347794108452577399262109977756586375*z^2)/641900390239748061475957626285817825407851726928630145741397950464 - 644430211103063219451233473971997200503671748328219055060769915199066015625/607305472086107978249312985781466984945489988373035614208, z, 2)

.

.

.

.

.

Answer 1

[Walter Roberson]

Notice the root() around the expression that includes z.

root() is a Symbolic Toolbox placeholder in which you pass a symbolic polynomial as the first parameter, and pass the name of the variable that is the base of the polynomial in the second parameter. Optionally there may be an index number in the third parameter.

root() with two parameters stands for the complete set of roots of the given polynomial -- it stands in for all the values, z, such that the polynomial becomes 0. Well, except when it doesn't: if the polynomial does not have constant roots, then root() with two parameters stands in for the "first" root of the polynomial.

root() with three parameters stands in for one particular root of the given polynomial. The order of the roots is not documented: all we know is that the order of roots is consistent for any one fully-numeric set of coefficients.

root() is a placeholder: it does not try to evaluate the roots. However, the symbolic toolbox is able to reason about the placeholder. For example you can diff() a root() structure and it will return something appropriate.

Thus,

root(z^10 - (76245256816524517094787721664575*z^8)/5352866887287026539823104 + (932018315641293413021053976330724993678573765595931929*z^6)/1232904261651835201667968301579019113988096 - (15149901962504148690356548006017832826970719790718281499456798538626475*z^4)/1334411437689202100645463103523731168093117650233720832 + (21560789318824412760547720965355965940321374664217347794108452577399262109977756586375*z^2)/641900390239748061475957626285817825407851726928630145741397950464 - 644430211103063219451233473971997200503671748328219055060769915199066015625/607305472086107978249312985781466984945489988373035614208, z, 2)

is root #2 of the degree 10 polynomial given in the expression.

If you look carefully at the expression you can see that the coefficients are zero for the odd numbered roots, so this can be reduced to +/- the roots of a degree 5 polynomial.

There is no closed form solution for most polynonials of degree 5, but your polynomial happens to be one of the few that can be solved. You can get approximations of the numeric solutions by using vpa()

You will find that you have 10 solutions, all of them real-valued. Three of the solutions are positive for both u_2 and beta.

However ...

When you use solve() you are asking for indefinitely-precise solutions, closed form solutions when possible. And that does not make sense to ask for when any of your coefficients is not exact. For example,

gamma=1.4

In science, that would mean that gamma is some particular but not-precisely-known value that is between 135/100 (inclusive) and 145/100 (exclusive). It makes no sense to ask what the exact half-million-character solution is for an equation when you do not know exactly what the inputs are!

It is a categorical error to use solve() with floating point inputs.

If you have assumptions available, I recommend that you add them. Be careful on your assumptions: one of the potential u_2 values is exactly 0.

[Hanrui Qiu]

Thank you so much for the explanation, it really helps

[Hanrui Qiu]

Hi Walter,

My beta should be 0-1 and u_2 should be a positive number, is there any way I can get the exact solution?

Thanks!

-Hanrui

[Walter Roberson]

bv = vpa(beta) ;
uv = vpa(u_2);
mask = bv >= 0 & uv >= 0;
Beta = bv(mask) ;
U_2 = uv(mask) ;