Solving for a variable within a loop?

Situation: Rocket nozzle design.

Code:

 % code
gamma=1.4;
R=287; %Gas constant
P0=22313; %Pascals
A_star=0.2; %Throat area in m^2

Ae=1:10;
Me=1;
for i=Ae
    AeAs(i)=Ae(i)/A_star;
    Ae_As(i)=(1/Me(i))*((2/(gamma+1))*(1+((gamma-1)/2)*(Me(i)^2)))^((gamma+1)/2*(gamma-1))==AeAs(i);
    solMe(i)=solve(Ae_As(i),Me(i));
end

I have to vary the area ratio of the exit area of the nozzle (Ae) to the area of the throat (A_star). A_star is fixed, A_e varies. (I have it varying from 1:10). So there are 10 different values for Ae/A_star. Using the equation labeled as "Ae_As(i)", I can solve for Me (exit Mach number). I tried to do this using the solve function which can be seen in the solMe line. When I run this code, I get this error:

    % code
Error using solve>processString (line 337)
' 1 ' is not a valid expression or equation.
Error in solve>getEqns (line 267)
      eqns = processString(eqns, v, vc);
Error in solve (line 150)
[eqns,vars,options] = getEqns(varargin{:});
Error in Problem_2 (line 14)
    solMe(i)=solve(Ae_As(i)==AeAs(i),Me(i));

Is there an easier way to do this? I've tried messing with it for hours and I'm getting nowhere.

Thanks in advance for your help!

Answer 1

[Star Strider]

There are some errors in your code. I corrected them and got it to run using vpasolve:

syms Me
gamma=1.4;
R=287; %Gas constant
P0=22313; %Pascals
A_star=0.2; %Throat area in m^2

Ae=1:10;
% Me=1;
for i=Ae
    AeAs(i)=Ae(i)/A_star;
    Ae_As=(1/Me)*((2/(gamma+1))*(1+((gamma-1)/2)*(Me^2)))^((gamma+1)/2*(gamma-1))==AeAs(i);
    solMe(i)=vpasolve(Ae_As,Me);
end

I declared ‘Me’ a symbolic variable to solve for rather than a constant, and then solved for it in the last assignment in the loop. The code runs. I leave it to you to determine if it produces the correct results, since the Mach numbers seem low to me. (However, I am not an aerodynamicist, so they could be entirely appropriate.)

[Manderson]

I want to thank you for taking a look at my code and fixing any mistakes I had. I plugged the code in and ran it, and it does work; but I am not seeing any values for Me.

You can see in the workspace that the code recognizes solMe is there, but when I click on it there is no values. If I were to individually click on one of those 10 "1x1 sym"s, it is blank. Perhaps I'm trying to view the Me's the wrong way?

[Star Strider]

My pleasure.

When I add this line after the end of the loop:

MeResult = solMe.'                          % Transpose To Column Vector

I get:

MeResult =
  0.18383459985826265283395122362604
 0.091694492484796863225796246272871
 0.061102269630989349396263477739046
 0.045819518157215486942531767587339
 0.036652955258240984702254447872134
 0.030542925794733257980955355798478
 0.026179028688427669920476267221981
 0.022906296887477889226099469250772
 0.020360937541054753031554424415671
  0.01832470522166758734727608032376

[Manderson]

That did the trick! I don't know if these are correct but I'll find out soon enough I guess. Thank you again for all of your help!

[Star Strider]

As always, my pleasure!

[Manderson]

I hate to bother you again but I have one final question regarding this problem. You were correct in your assumptions, the Me results are too small, and I have found out why.

    % code
syms Me
gamma=1.4;
R=287; %Gas constant
k=1.5e-04;
h=10000; %Height in meters
P0=100e3; %Initial Pressure in Pascals
Pt=6e6; %Total Pressure in Pascals
Tt=2500; %Total Temperature in Kelvin
m_dot=970.042; %Mass flux in kg/s
A_star=0.2; %Throat area in m^2
Patm=P0*exp(-k*h);
Ae=1:10;
for i=Ae
    AeAs(i)=Ae(i)/A_star;
    Ae_As(i)=(1/Me)*(((2/(gamma+1))*(1+((gamma-1)/2)*(Me^2)))^((gamma+1)/2*(gamma-1)))==AeAs(i);
    solMe(i)=vpasolve(Ae_As(i),Me);

    Te(i)=Tt/(1+((gamma-1)/2)*solMe(i));
    Pe(i)=Pt/((Tt/Te(i))^(gamma/(gamma-1)));
    ue(i)=((gamma*R*Te(i))^0.5)*solMe(i);
    F(i)=m_dot*ue(i)+((Pe(i)-Patm)*Ae(i));
end
MeResult=solMe.'
TeResult=Te.'
PeResult=Pe.'
ueResult=ue.'
FResult=F.'
end

This is my full code (never mind the Te, Pe, ue, and F lines). The Me values are so small because in the Ae_As(i) line, it's reading it as ((gamma+1)/2)*(gamma-1), which equals 12/25. If you unsuppress this line you can actually see it calculate it and it'll show that it's reading it as (12/25). No problem I thought, I'll just throw some parentheses in there so it reads it in the proper order. That would be changing it to:

    % code
(gamma+1)/(2*(gamma-1))
  end

Well when I do that MATLAB has a fit and I get the error:

    % code
Error using mupadmex
Error in MuPAD command: Subscripted assignment dimension mismatch
Error in sym/subsasgn (line 1450)
                C = mupadmex('symobj::subsasgn',A.s,B.s,inds{:});
Error in Problem_2 (line 19)
    solMe(i)=vpasolve(Ae_As(i),Me);
  end

I've tried re-writting it in I don't know how many ways, and even just hardcoding the actual value, 3 (2.4/.4=3) and I still get the same error. Any ideas?

Thank you so much again. This should be the last problem I have because everything else works accordingly.

[Star Strider]

As always, my pleasure.

The problem is that this line (and I believe I copied it correctly, since I had to insert your new (gamma+1)/(2*(gamma-1)) into it) now produces a (6x1) vector, so you have to subscript it and the subsequent assignments appropriately:

Ae_As=(1/Me)*(((2/(gamma+1))*(1+((gamma-1)/2)*(Me^2))))^((gamma+1)/(2*(gamma-1)))==AeAs(i);
solMe(:,i)=vpasolve(Ae_As,Me);
Te(:,i)=Tt/(1+((gamma-1)/2)*solMe(:,i));
Pe(:,i)=Pt/((Tt/Te(i))^(gamma/(gamma-1)));
ue(:,i)=((gamma*R*Te(i))^0.5)*solMe(:,i);
F(:,i)=m_dot*ue(i)+((Pe(i)-Patm)*Ae(:,i));

The new subscripting (:,i) now stores it. It is not a single scalar.

You don’t need to store ‘AeAs’ or ‘Ae_As’ (unless you want to) since they are one-offs in each iteration and you don’t use them outside the loop. I removed the subscripts from ‘AeAS’ because I wanted to see what it was and why it wasn’t working in the following line.

The code now throws a number of warnings about rank-deficient matrices. I can’t interpret the outputs, but they are all numeric.

I get the impression this is not the direction you want your code to go. I know nothing about rocket engines (except that they look cool when they work correctly), so I can’t help you troubleshoot your ‘Ae_As’ assignment. I would look at it very carefully to be certain you coded it correctly.

I tend to avoid long lines of code for the very reason that they are difficult to troubleshoot, especially when they involve large numbers of nested parentheses. Consider breaking the ‘Ae_As’ line into smaller segments, and combining them in the final ‘Ae_As’ assignment. (There’s nothing wrong with that, and it’s not inefficient since the computer has to do all the parsing and calculations anyway. Breaking it up simply makes it easier for you to parse it, and to be sure it’s parsed correctly.) Then you can be sure the individual segments are correct, and you can avoid a lot of the nested parentheses that make it difficult to work with.

I’ll be glad to continue to help you on the MATLAB code.

[Manderson]

Here's a picture of the equations written out (I understand that they aren't the easiest to read when written in MATLAB form). While I don't know whether or not the code is going in the wrong direction, I get the feeling I approached it wrong initially and that there is possibly an easier way to go about doing it.

Initially I thought I could just create a loop for Ae ranging from 1 to 10, then solve for Me for each of those ten cases, and then use that Me to solve for the other variables (Te, Pe, ue, and F). I don't have time at the moment to go back and carefully comb through the code but I'll continue to fiddle with it. If I get it working I'll let you know! Then you can tell people you know how to optimize a rocket engines thrust for given parameters, haha.

[Star Strider]

When I code your top equation:

syms Me 
Ae = 1;
gamma = 1.4;
A_star = 0.2;
Ae_As = (1/Me) * ( 2/(gamma+1) * (1 + (gamma-1)/2 * Me^2) )^(gamma+1)/(2*(gamma-1)) == Ae/A_star;
solMe = vpasolve(Ae_As, Me)

solMe =
3.6525399777641596497340539217213

Does that look correct? If it does, put my equation in your loop, making the appropriate substitutions. Use a scalar subscript ‘Ae_As(i)’ since the output is a scalar. I put significant ‘whitespace’ in my MATLAB expression so I could easily read the nested parentheses and be certain they are correct.

[Manderson]

Yes, I think that is correct! That Mach numbers look appropriate except for Ae=3 and and Ae=10 which is weird, but I will check those by hand later. I will start adding more whitespace from now on when dealing with long lines of code. Thank you again for all of your help!

[Star Strider]

My pleasure, as always!

Other possibilities to consider are the core MATLAB functions numerical solvers fzero and fminsearch, and fsolve in the Optimization Toolbox. Those actually might be easier to use in your application than the Symbolic Toolbox. They would definitely be faster, because the Symbolic Math Toolbox is not designed for such iterative calculations. They require initial parameter estimates, so if you know approximately what they should be, using those functions should be straightforward.

[Manderson]

I found out why the numbers are so low for Ae(3) and Ae(10). I tried solving for Me by hand by setting Ae/A_star=15 and then solving for Me using Wolfram Alpha. It turns out there are two solutions for Me when Ae/A_star=15, one is 0.0539 (the number that MATLAB was outputting), and 5.395 (the correct Mach number). If there are multiple solutions to a function, is there any way I can tell MATLAB to choose one solution as opposed to another?

I was thinking along the lines of putting a statement towards the beginning of the loop saying something like: 'if Me>2 store, else discard', something maybe like:

    % code
for i=Ae
    AeAs(i)=Ae(i)/A_star;
    Ae_As(i)=(1/Me) * ( 2/(gamma+1) * (1 + (gamma-1)/2 * Me^2) )^(gamma+1)/(2*(gamma-1)) == Ae(i)/A_star;
    solMe(i)=vpasolve(Ae_As(i),Me);
        if solMe(i)>2==solMe(i)*1;
        else solMe(i)==solMe(i)*0;
    Te(i)=Tt/(1+((gamma-1)/2)*solMe(i));
    Pe(i)=Pt/((Tt/Te(i))^(gamma/(gamma-1)));
    ue(i)=((gamma*R*Te(i))^0.5)*solMe(i);
    F(i)=m_dot*ue(i)+((Pe(i)-Patm)*Ae(i));
        end
end

[Star Strider]

If you know the maximum is always going to be the correct solution, I would simply put (in place of the if block):

solMev = vpasolve(Ae_As(i),Me);
solMe(i) = max(solMev);

With ‘solMev’ being the vector of solutions. You have to allow for the possibility of multiple solutions (that will be a column vector) so I would not subscript the first ‘solMe’ assignment. Let it be whatever it is, and choose the maximum in the array assignment for ‘solMe(i)’ in the following line. That does not affect scalar values of ‘solMe’, since if there is only one value, it is the maximum.

I would also not store ‘AeAs’ or ‘Ae_As’, since your never refer to them again. (If you need to store them for other reasons, ignore this suggestion.)

See how this works. I’ll continue to help.

[Manderson]

Okay, here's my code:

    % code
for i=Ae
    AeAs=Ae(i)/A_star;
    Ae_As=(1/Me) * ( 2/(gamma+1) * (1 + (gamma-1)/2 * Me^2) )^(gamma+1)/(2*(gamma-1)) == Ae(i)/A_star;
    solMe=vpasolve(Ae_As(i),Me);
    solMev=vpasolve(Ae_As(i),Me);
    solMe(i) = max(solMev);
    Te(i)=Tt/(1+((gamma-1)/2)*solMe(i));
    Pe(i)=Pt/((Tt/Te(i))^(gamma/(gamma-1)));
    ue(i)=((gamma*R*Te(i))^0.5)*solMe(i);
    F(i)=m_dot*ue(i)+((Pe(i)-Patm)*Ae(i));
end

I believe I am no longer storing AeAs and Ae_As, and after I inserted the two lines of code you had suggested, I get an error saying: "Undefined function 'max' for input arguments of type 'sym'. Error in Problem_2 (line 21) solMe(i) = max(solMev);". Again, I thank you for your dedication.

[Star Strider]

As always, my pleasure.

Delete this line because it could be setting ‘solMe’ as a symboic variable (and you don’t need it anyway):

solMe=vpasolve(Ae_As(i),Me);

If necessary (after deleting that line), change the ‘solMe(i)’ assignment to:

solMe(i) = max(double(solMev));

That should work.

[Manderson]

Well, the code runs now but it's still choosing the smaller of the two values for Ae=3 and Ae=10. Here is the full code:

    % code
syms Me
gamma=1.4;
R=287; %Gas constant
k=1.5e-04;
h=10000; %Height in meters
P0=100e3; %Initial Pressure in Pascals
Pt=6e6; %Total Pressure in Pascals
Tt=2500; %Total Temperature in Kelvin
m_dot=970.042; %Mass flux in kg/s
A_star=0.2; %Throat area in m^2
Patm=P0*exp(-k*h);
Ae=1:10;
for i=Ae
    AeAs=Ae(i)/A_star;
    Ae_As=(1/Me) * ( 2/(gamma+1) * (1 + (gamma-1)/2 * Me^2) )^(gamma+1)/(2*(gamma-1)) == Ae(i)/A_star
    solMev=vpasolve(Ae_As,Me);
    solMe(i) = max(double(solMev));
    Te(i)=Tt/(1+((gamma-1)/2)*solMe(i));
    Pe(i)=Pt/((Tt/Te(i))^(gamma/(gamma-1)));
    ue(i)=((gamma*R*Te(i))^0.5)*solMe(i);
    F(i)=m_dot*ue(i)+((Pe(i)-Patm)*Ae(i));
end
MeResult=solMe.'
TeResult=Te.'
PeResult=Pe.'
ueResult=ue.'
FResult=F.'
  end

[Star Strider]

I’m only getting one (wrong) value for those. In desperation, I turned to fzero with this code:

Me3 = fzero(@(Me) (5*(Me^2/6 + 5/6)^(12/5))/(4*Me) - 15, 10)
Me10 = fzero(@(Me) (5*(Me^2/6 + 5/6)^(12/5))/(4*Me) - 50, 10)

Me3 =
     5.3948

Me10 =
     7.7861

I will let you experiment with those.

If the correct values in the symbolic solutions are always the second values, instead of using max, try this:

solMe(i) = solMev(end);

That will take the last (or only) value of the vector and return it to ‘solMe(i)’.

[Manderson]

Still no luck using solMe(i) = solMev(end);, that is, the code is only outputting the smaller of the two values. I'll keep messing with it but if I don't get anywhere I might just abandon it. I have so many other things to do that I can't afford to get caught up on a small thing like this. Again, I'll let you know if I do end up getting it to work.

You've been a great help and I have learned a few things along the way. I cannot thank you enough for all you have done. I don't know one person who would have put so much time and effort into a code they had no attachment to, haha. Take care :)

[Star Strider]

My pleasure, as always!

Rather than use the Symbolic Toolbox, use fzero. I would. It works! It’s also much faster and more efficient.

Solving problems like this and explaining them are fun for me (I’m one of the unpaid volunteers here), and I always learn things. I encounter problems I probably would never see in my own code or my own areas of research and expertise.

Meanwhile, submit a Bug Report to MathWorks about vpasolve returning a lower of two values, and your not being able to get the correct value, while fzero produces the correct result. Include the URL to this thread: http://www.mathworks.com/matlabcentral/answers/254029-solving-for-a-variable-within-a-loop. This thread tells the entire story, so you need do nothing more than briefly describe the problem and include the URL.