Getting a single value out of vpasolve with growing matrix

Question

Vidhan Malik 2016-4-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/279250-getting-a-single-value-out-of-vpasolve-with-growing-matrix

评论： Star Strider 2016-4-16

The equation "E" has two roots and I am only interested in the positive root. the function for To5 keeps growing to a predefined point, my trouble is that I'm getting this error occurring at the second line of code: "Assignment has more non-singleton rhs dimensions than non-singleton subscripts"

E = (1+R)*(Ym*Rg)*To5(i,2) + fo/(1+R)*Qr + (1+R+fo/(1+R) == 0.12*(1+R))*Yc*Rg*To6;
R(i,2) = vpasolve(E,R, [0 Inf]); %Step 6

I believe the problem is due to my not using vpasolve correctly in getting just the positive value out of the function and then assigning it to the growing matrix R, know how I can fix this?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2016-4-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/279250-getting-a-single-value-out-of-vpasolve-with-growing-matrix#answer_218112

在 MATLAB Online 中打开

You didn’t post the rest of your code, and it doesn’t work with the code you previously posted, so I can only post an example:

syms x
R = vpasolve(x^2 - x - 4 == 0)
ReR = R(real(R) > 0)
R =
 -1.561552812808830274910704927987
  2.561552812808830274910704927987
ReR =
2.561552812808830274910704927987

It’s relatively easy to select the solution that is greater than zero, or any number of other conditions.

38 个评论
显示 36更早的评论隐藏 36更早的评论

Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

Hello again!

clear; clc;
syms R
To7 = 1100;
To6 = 1800;
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
for i = 1:1:100
      To7 = To7+1;
      PrT = (To7/To6)^(Yc/0.92); %Step 1
      PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
      To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
      To5(i,2) = 0.92*(To7 - To4) + To4; %Step 4
      fostoic = 14/(1.5*(32+3.76*28)); 
      fo = 0.9*fostoic; %Step 5
      E = (1+R)*(Ym*Rg)*To5(i,2) + fo/(1+R)*Qr + (1+R+fo/(1+R) == 0.12*(1+R))*Yc*Rg*To6;
      R(i,2) = vpasolve(E,R, [0 Inf]); %Step 6
      Wnet(i,2) = (1+R(i,2)+fo/(1+R(i,2)) - 0.12*(1+R(i,2)))*Yc*Rg*(To6-To7) - (1+R(i,2))*Ym*Rg*(To4-T32); %Step 7
      Qin(i,2) = (1+R(i,2)+fo/(1+R(i,2))-0.12*(1+R(i,2)))*(Yc*Rg*To6-Ym*Rg*To5(i,2));
      Eff = Wnet(i,2)/Qin(i,2)*100; %Step 8
end
disp(Eff);

I tried using anonymous functions as you tried to teach me before but I honestly still don't fully grasp them and how to implement them properly. That is why I went back to doing it this way.

Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

That is very peculiar that you are getting non-real answers. I have a similar iteration of the previous code that gives me real values. I think for now I rather just use this method instead of using anonymous functions. Once I have gotten far enough I can always revisit the code and see what improvements I can make.

clear; clc;
syms R
To6 = 1800;
To7 = To6*(2/3);
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
PrT = (To7/To6)^(Yc/0.92); %Step 1
PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
To5 = 0.92*(To7 - To4) + To4; %Step 4
fostoic = 14/(1.5*(32+3.76*28)); 
fo = 0.9*fostoic; %Step 5
E = (1+R)*(Ym*Rg)*To5 + fo/(1+R)*Qr == (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
R = vpasolve(E,R,[0 Inf]); %Step 6
Wnet = (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*(To6-To7) - (1+R)*Ym*Rg*(To4-T32); %Step 7
Qin = (1+R+fo/(1+R)-0.12*(1+R))*(Yc*Rg*To6-Ym*Rg*To5);
Eff = Wnet/Qin*100; %Step 8
disp(Eff);

The the end of the day I want to plot Eff(efficiency) vs To6 where:

To6 = 1800;
To7 = To6*(2/3);

I would love to explain to you the equations but they all have technical derivations of gas turbines equations so it may take a while :). But essentially, I have commented each of the steps that define an important parameter which is to be used in a later line of code.

Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

Using the previous second code I posted, I have tried doing it through using anonymous functions as follows:

      E = @(R) -(1+R)*(Ym*Rg)*To5 - fo/(1+R)*Qr + (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
      R = fzero(E, 0);

In the code, everything besides R is a variable predefined so we just need to find the positive root of R.

Here is the mess of the code I have currently.

clear; clc;
syms R
To6 = 999;
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
To6(1,2) = 999
for i = 2:1:801
      To6(i,2) = To6(i-1,2) +1;
      To7 = To6(i,2)*(2/3);
      PrT = (To7/To6(i,2))^(Yc/0.92); %Step 1
      PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
      To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
      To5(i,2) = 0.92*(To7 - To4) + To4; %Step 4
      fostoic = 14/(1.5*(32+3.76*28)); 
      fo = 0.9*fostoic; %Step 5
      %E = (1+R)*(Ym*Rg)*To5(i,2) + fo/(1+R)*Qr + (1+R+fo/(1+R) == 0.12*(1+R))*Yc*Rg*To6;
      %R(i,2) = fsolve(E,R, [0 Inf]); %Step 6
      E = @(R) -(1+R)*(Ym*Rg)*To5 - fo/(1+R)*Qr + (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
      R = fzero(E, 0);
      Wnet(i,2) = (1+R(i,2)+fo/(1+R(i,2)) - 0.12*(1+R(i,2)))*Yc*Rg*(To6(i,2)-To7) - (1+R(i,2))*Ym*Rg*(To4-T32); %Step 7
      Qin(i,2) = (1+R(i,2)+fo/(1+R(i,2))-0.12*(1+R(i,2)))*(Yc*Rg*To6(i,2)-Ym*Rg*To5(i,2));
      Eff = Wnet(i,2)/Qin(i,2)*100; %Step 8
end
disp(Eff);

I have tried setting up the loop so that the value of To6 changes from 1000 to 1800 by 1. Also, To6 is the only independent variable.

Vidhan Malik 2016-4-16

编辑：Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

The reason why I ask how to specify which root that the "fsolve" function is outputting is because I have a similar code to the one with the loop, where there is no loop. Through that I can check if my answers from the loop are correct or not, and so far they are not.

Here is the code without the loop. In this code, I am using vpasolve and using the positive root for 'R' in step 6

clear; clc;
syms R
To6 = 1700;
To7 = To6*(2/3);
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
PrT = (To7/To6)^(Yc/0.92); %Step 1
PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
To5 = 0.92*(To7 - To4) + To4; %Step 4
fostoic = 14/(1.5*(32+3.76*28)); 
fo = 0.9*fostoic; %Step 5
E = (1+R)*(Ym*Rg)*To5 + fo/(1+R)*Qr == (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
R = vpasolve(E,R,[0 Inf]); %Step 6
Wnet = (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*(To6-To7) - (1+R)*Ym*Rg*(To4-T32); %Step 7
Qin = (1+R+fo/(1+R)-0.12*(1+R))*(Yc*Rg*To6-Ym*Rg*To5);
Eff = Wnet/Qin*100; %Step 8
disp(R);

By varying the value of To6 initially, I can find the impact that it has on the value of "R" and "Eff". So I just want to turn this code into a loop where To6 is varied from 1700 to 1900.

Star Strider 2016-4-16

在 MATLAB Online 中打开

In that example, with on positive root at about 1 and a negative root at about -3, giving it a positive initial estimate would find the positive root.

However, it seems easier than all that. I looked at your ‘E’ expression and noted that it is quadratic in ‘R’, making identification of the positive root straightforward, since all the constants are positive. After commenting-out the constant assignments and adding them to the syms call, the symbolic expressions for the two roots were the same except for sign. The positive root is:

E = (1+R)*(Ym*Rg)*To5 + fo/(1+R)*Qr == (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
Rs = solve(E, R);
Rspos = simplify(Rs(2), 'steps',10);
R_fcn = matlabFunction(Rspos)
R_fcn = @(Qr,Rg,To6,To7,Yc,Ym) (sqrt(Rg.*(Qr-Rg.*To6.*Yc).*(To6.*Yc.*-2.2e1+To7.*Ym.*2.3e1+Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2).*-5.005e3).*(5.0./2.86e2)+Rg.*To6.*Yc.*2.2e1-Rg.*To7.*Ym.*2.3e1-Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2)./(Rg.*To6.*Yc.*-2.2e1+Rg.*To7.*Ym.*2.3e1+Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2);

So, just substitute in the appropriate arguments in ‘R_fcn’ and the result will always be the positive root. You can then use it in the rest of your code. I checked the value of ‘Rss’, and it was 0.9900534780448165117718294123644, so it will be positive with these constants.

If you don’t want (or need) to use the function form of ‘R’, you can just assign it directly as:

R = ((5*(-5005*Rg*(Qr - Rg*To6*Yc)*(23*To7*Ym - 22*To6*Yc + 552*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym))))^(1/2))/286 + 22*Rg*To6*Yc - 23*Rg*To7*Ym - 552*Rg*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym)))/(23*Rg*To7*Ym - 22*Rg*To6*Yc + 552*Rg*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym)));

since there is no longer a need to solve for it.

This all comes directly from your code. All I did was use the Symbolic Math Toolbox to solve it once symbolically, so solving it either symbolically or numerically in each iteration with different variables is no longer necessary.

Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

I wish I could cook something for you! haha, thank you so much for all your help once again! Here is my final code that I think is giving me the proper results, if you fancy a look!

clear; clc;
syms R
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
To6(1) = 1600;
i = 1;
Eff(1) = 45.844052030815569938093624771305;
while i <300
      i = i +1;
      To6(i) = To6(i-1) +1;
      To7 = To6(i)*(2/3);
      PrT = (To7/To6(i))^(Yc/0.92); %Step 1
      PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
      To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
      To5 = 0.92*(To7 - To4) + To4; %Step 4
      fostoic = 14/(1.5*(32+3.76*28)); 
      fo = 0.9*fostoic; %Step 5
      E = @(To5, To6,R) -(1+R)*(Ym*Rg)*To5 - fo/(1+R)*Qr + (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
      Eqn = @(R) E(To5, To6(i),R);
      R = fsolve(Eqn, 0.55);
      Wnet(i) = (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*(To6(i)-To7) - (1+R)*Ym*Rg*(To4-T32); %Step 7
      Qin(i) = (1+R+fo/(1+R)-0.12*(1+R))*(Yc*Rg*To6(i)-Ym*Rg*To5);
      Eff(i) = Wnet(i)/Qin(i)*100; %Step 8
      disp(R);
end
%figure(1)
plot(To6(:), Eff(:))
grid on
xlabel('To6 (K)')
ylabel('Thermal Efficiency (%)')

And the code gives me this beautiful graph!

Star Strider 2016-4-16

在 MATLAB Online 中打开

As always, my pleasure!

I used this version of your code:

syms R To6 To7 Ya Ym Yc Rg Qr T32
% To6 = 1700;
% To7 = To6*(2/3);
% 
% Ya = 1.4/0.4;
% Ym = 1.35/0.35;
% Yc = 1.3/0.3;
% Rg = 0.287;
% Qr = 45e3;
% T32 = 276;
PrT = (To7/To6)^(Yc/0.92); %Step 1
PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
To5 = 0.92*(To7 - To4) + To4; %Step 4
fostoic = 14/(1.5*(32+3.76*28)); 
fo = 0.9*fostoic; %Step 5
E = (1+R)*(Ym*Rg)*To5 + fo/(1+R)*Qr == (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
Es = simplify(E, 'steps',10);
Rs = solve(E, R);
Rss = simplify(Rs(2), 'steps',10)
Rspos = vpa(Rss)
R_fcn = matlabFunction(Rspos)
% R = vpasolve(E,R,[0 Inf]); %Step 6
Wnet = (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*(To6-To7) - (1+R)*Ym*Rg*(To4-T32); %Step 7
Qin = (1+R+fo/(1+R)-0.12*(1+R))*(Yc*Rg*To6-Ym*Rg*To5);
Eff = Wnet/Qin*100; %Step 8
% disp(R);
% R_fcn = @(Qr,Rg,To6,To7,Yc,Ym) (sqrt(Rg.*(Qr-Rg.*To6.*Yc).*(To6.*Yc.*-2.2e1+To7.*Ym.*2.3e1+Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2).*-5.005e3).*(5.0./2.86e2)+Rg.*To6.*Yc.*2.2e1-Rg.*To7.*Ym.*2.3e1-Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2)./(Rg.*To6.*Yc.*-2.2e1+Rg.*To7.*Ym.*2.3e1+Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2);

I commented-out all the constants and added their assignments to the syms call so they would remain symbolic for my analysis. Then I did a symbolic solution. I looked at the results of ‘Rs’ (‘R solved’), saw that the second one was positive, and addressed it as such, since the solution is a vector of the two solutions. (It’s not possible to test for it, because unless you list a number of assumptions, the Symbolic Math Toolbox has no way of knowing which root will be positive, depending on the values of the variables involved.) I then simplified it to create ‘Rspos’ (‘R solved positive’). I used the matlabFunction function to create ‘R_fcn’, the anonymous function representation of it, but then it occurred to me that you were using the anonymous function to solve for the root, and since ‘Rspos’ and ‘R_fcn’ were the solutions, the anonymous funciton wasn’t necessary. You only need to use ‘Rspos’ that I renamed ‘R’ and copied to my earlier Comment as:

R = ((5*(-5005*Rg*(Qr - Rg*To6*Yc)*(23*To7*Ym - 22*To6*Yc + 552*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym))))^(1/2))/286 + 22*Rg*To6*Yc - 23*Rg*To7*Ym - 552*Rg*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym)))/(23*Rg*To7*Ym - 22*Rg*To6*Yc + 552*Rg*Ym*(4503599627370496000/(4050119001849345243*(To7/To6)^((25*Yc)/23)))^(10/(9*Ym)));

So, one symbolic solution and you’re done! You don’t need to solve it in each iteration. I wish I’d thought of that earlier, but it took me a while to understand what you want to do.

Vidhan Malik 2016-4-16

在 MATLAB Online 中打开

I tried updating my code using that method but can't tell if it is more efficient or not.

clear; clc;
syms R To5 fo To6
Ya = 1.4/0.4;
Ym = 1.35/0.35;
Yc = 1.3/0.3;
Rg = 0.287;
Qr = 45e3;
T32 = 276;
To6(1) = 1600;
i = 1;
Eff(1) = 45.844052030815569938093624771305;
E = (1+R)*(Ym*Rg)*To5 + fo/(1+R)*Qr == (1+R+fo/(1+R) - 0.12*(1+R))*Yc*Rg*To6;
Rs = solve(E, R);
Rspos = simplify(Rs, 'steps',10);
Rr = matlabFunction(Rspos);
Rr = @(Qr,Rg,To6,To7,Yc,Ym) (sqrt(Rg.*(Qr-Rg.*To6.*Yc).*(To6.*Yc.*-2.2e1+To7.*Ym.*2.3e1+Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2).*-5.005e3).*(5.0./2.86e2)+Rg.*To6.*Yc.*2.2e1-Rg.*To7.*Ym.*2.3e1-Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2)./(Rg.*To6.*Yc.*-2.2e1+Rg.*To7.*Ym.*2.3e1+Rg.*Ym.*((To7./To6).^(Yc.*(-2.5e1./2.3e1)).*1.111967234867441).^((1.0e1./9.0)./Ym).*5.52e2);
while i <300
      i = i +1;
      To6(i) = To6(i-1) + 1;
      To6 = To6(i);
      To7 = To6*(2/3);
      PrT = (To7/To6)^(Yc/0.92); %Step 1
      PrC = 1/(0.97*0.99^3*PrT*0.975*0.98); %Step 2
      To4 = 276*PrC^(1/(0.9*Ym)); %Step 3
      To5 = 0.92*(To7 - To4) + To4; %Step 4
      fostoic = 14/(1.5*(32+3.76*28)); 
      fo = 0.9*fostoic; %Step 5
      Wnet(i) = (1+Rr(Qr,Rg,To6,To7,Yc,Ym)+fo/(1+Rr(Qr,Rg,To6,To7,Yc,Ym)) - 0.12*(1+Rr(Qr,Rg,To6,To7,Yc,Ym)))*Yc*Rg*(To6-To7) - (1+Rr(Qr,Rg,To6,To7,Yc,Ym))*Ym*Rg*(To4-T32); %Step 7
      Qin(i) = (1+Rr(Qr,Rg,To6,To7,Yc,Ym)+fo/(1+Rr(Qr,Rg,To6,To7,Yc,Ym))-0.12*(1+Rr(Qr,Rg,To6,To7,Yc,Ym)))*(Yc*Rg*To6-Ym*Rg*To5);
      Eff(i) = Wnet(i)/Qin(i)*100; %Step 8
      disp(vpa(Rr(Qr,Rg,To6,To7,Yc,Ym)));
end
figure(1)
plot(To6(:,2), Eff(:,2))
grid on

Either way I am getting the error of Index exceeds matrix dimensions. in line 27:

To6(i) = To6(i-1) + 1;

I don't see where it exceeds the matrix dimension.

Star Strider 2016-4-16

在 MATLAB Online 中打开

I do!

Your ‘To6’ value is a scalar:

To6(1) = 1600;
. . .
        To6(i) = To6(i-1) + 1;                      % <— DEFINED AS A VECTOR ...
        To6 = To6(i);                               % <— ... THEN REDEFINED AS A SCALAR

For scalars, any subscript greater than 1 will throw the error you saw.

One solution is to define two different variables before the loop:

To6(1) = 1600;
To6v = To6;

... then within the loop:

        To6v(i) = To6v(i-1) + 1;
        To6 = To6v(i);                             % <— ADD THIS LINE

That will run. (I tested it.) I will leave you to determine if it produces the correct result.

请先登录，再进行评论。

Getting a single value out of vpasolve with growing matrix

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

38 个评论
显示 36更早的评论隐藏 36更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Getting a single value out of vpasolve with growing matrix

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

38 个评论 显示 36更早的评论隐藏 36更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

38 个评论
显示 36更早的评论隐藏 36更早的评论