0

I'm trying to solve the Machnumber-Area equation for isentropic flows. This equation has two real positive roots, but I'm not able to get any of these when trying to solve with the "solve" routine.

The code to solve the equation:

gamma = 1.4;
L     = 3;
n     = 30;

dx = L/n;
x  = 0;
for i=1:n+1
A(i)  = 1 + 2.2*( x - 1.5 )^2;
X(i)  = x;
x     = x + dx;
end
A0 = min(A);

for i=1:1 %n1
syms M    
eqn  = (A(i)/A0)^2 == ( 1/M^2 )*( ( 2/(gamma+1) )*( 1+0.5*(gamma-1)*M^2 ) )^( (gamma+1)/(gamma-1) )
solM = solve(eqn,M)
% Mach_is(i) = solM;
end

I commented out to loop through the entire vector since it crashes do the the output I get:

  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 1))^(1/2)
  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 2))^(1/2)
  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 3))^(1/2)
  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 4))^(1/2)
  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 5))^(1/2)
  (1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 6))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 1))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 2))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 3))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 4))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 5))^(1/2)
 -(1/root(z^6 - (40824726*z^5)/390625 + (3*z^4)/5 + (4*z^3)/25 + (3*z^2)/125 + (6*z)/3125 + 1/15625, z, 6))^(1/2)

Can anyone suggest why am I getting this and how to solve the equation?

Adriaan
  • 17,741
  • 7
  • 42
  • 75
timbjörn
  • 39
  • 5
  • 2
    The equation is solved, there is nothing wrong here. But if you would like to evaluate the symbolic solution in order to get a numerical value, use `vpa(solM)`. You should also use the function `assume` to add some constraints (non complex and positive solution). And your `for loop` is useless since you only use `A(1)`. – obchardon Jun 17 '20 at 08:26
  • @obchardon I was not aware of `vpa` and `assume`, but as you say assuming real and positive and taking `vpa` of that gives me the two roots. Thanks! (As I wrote above I took the first element of the vector only since I was not able to get and fill up a vector with the solution.) – timbjörn Jun 17 '20 at 08:36
  • @Adriaan No I only solves it once but I get the two real roots both positive and negative, which probably makes sense, and apparently 4 additional complex roots, both positive and negative. With the suggestions above I'm now able to sort out the correct roots though. – timbjörn Jun 17 '20 at 08:38

0 Answers0