Solve Equation Using a Loop?

Hello, I am trying to solve an equation for a specific variable using a loop. The equation is a bit long, so I've shortened it using a, b, c, d, so bear with me. The equation is:

if true 
Tau = 1.8378
a = 1+gamma*M1.^2;
b = 1+gamma*M2.^2;
c = 1+((gamma-1)/2)*M2.^2;
d = 1+((gamma-1)/2)*M1.^2;
M2sol = (a./b).*((M2./M1).^2).*(c./d) == Tau;
M2 = solve(M2sol,M2)
end

M1 is an array starting at 0.0, increasing in increments of 0.01, all the way up to 0.4.

I want to solve the M2sol equation for M2 using the various values for M1. How do I accomplish this using a loop? So far I have:

if true
M1 = 0:0.01:0.4;
Tau = 1.8378
for i=1:41
  syms M2
  a = 1+gamma*M1.^2;
  b = 1+gamma*M2.^2;
  c = 1+((gamma-1)/2)*M2.^2;
  d = 1+((gamma-1)/2)*M1.^2;
  M2sol = (a./b).*((M2./M1).^2).*(c./d) == Tau;
  M2(i) = solve(M2sol,M2)
end
end

Clearly this isn't right; how do I fix it? I'm guessing I'm screwing up something in with the indexing? And probably far too many periods in that Msol equation...

Either way, any help would be greatly appreciated!

Answer 1

[Roger Stafford]

1. In computing a, d, and M2sol you need to use M1(i) instead of M1.

2. Using the ‘solve’ function here has the problem that there will be four possible solutions for M2 at each step and we have no idea which one matlab will select to place in the M2 array. You can see this from the fact that the equation can easily be converted to a quadratic equation in M2^2 which would have two roots for M2^2 and therefore four roots for M2. If you solve the equation this way you then have the option of selecting which of the four roots you want placed in M2 array. For this reason I recommend that you use this method instead of ‘solve’.

[Manderson]

Would it be easier if I initialed M1 as 0.0 and then wrote in the loop:

 if true
M1=0.0
  for i=1:41
   syms M2
   a = 1+gamma*M1.^2;
   b = 1+gamma*M2.^2;
   c = 1+((gamma-1)/2)*M2.^2;
   d = 1+((gamma-1)/2)*M1.^2;
   M2sol = (a./b).*((M2./M1).^2).*(c./d) == Tau;
   M2(i) = solve(M2sol,M2)
   M1(i) = M1+0.01
  end

Otherwise I will have to use the quadratic way. Speaking of, I'm a little confused on how to get a quadratic out of that equation. I tried expanding M2sol out and got a fraction that doesn't seem very helpful. How would you go about doing that?

[Roger Stafford]

As to your question of how to get a quadratic equation, your original equation is

   (1+gamma*M1^2)/(1+gamma*M2^2)*(M2/M1)^2*(1+(gamma-1)/2*M2^2)/(1+(gamma-1)/2*M1^2) = tau

If you first transpose tau to the left side and then multiply by

   (1+gamma*M2^2)*M1^2*(1+(gamma-1)/2*M1^2),

after collecting like powers of M2 and simplifying, you will get an equation of the form

   A*(M2^2)^2 + B*(M2^2) + C = 0

where

   A = (gamma-1)/2*(1+gamma*M1^2)
   B = 1+gamma*M1^2*(1-(1+(gamma-1)/2*M1^2)*tau)
   C = -M1^2*(1+M1^2*(gamma-1)/2)*tau

This is a quadratic equation in the unknown M2^2 so its two roots are

   M2^2 = (-B+sqrt(B^2-4*A*C))/(2*A)

and

   M2^2 = (-B-sqrt(B^2-4*A*C))/(2*A)

Taking the plus and minus square roots of these gives you four possible roots for M2, though some of them may be complex-valued. It will be up to you to select the one out of the four possible which you regard as the appropriate solution to be placed in the M2 array.