Getting a single value out of vpasolve with growing matrix

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?

Answer 1

[Star Strider]

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.

[Vidhan Malik]

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.

[Star Strider]

‘The equation "E" has two roots and I am only interested in the positive root.’

E =
(14992954385710700973*R)/4398046511104000 + 1653351/(572*(R + 1)) + 14992954385710700973/4398046511104000 == (6328943496324176301*R)/4398046511104000 + 393750/(143*(R + 1)) + 6328943496324176301/4398046511104000

R =
 - 1.0 - 0.26369061414332790360591577294225i
 - 1.0 + 0.26369061414332790360591577294225i

‘Houston, we have a problem ...’

There are no positive real roots, at least in this iteration, and possibly others.

Anonymous functions are relatively straightforward, at least for simple problems. I will be glad to help you understand them if they will work in your intended application.

I find it a bit difficult to follow your code.

What do you mean by ‘positive’, (real or imaginary parts) and do you want to do?

[Vidhan Malik]

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]

The first code I posted had a mistyped value of To7, it should have been 1200 not 1100.

[Star Strider]

I still get complex values for ‘R’ with that change.

[Vidhan Malik]

Have you tried using the second code without the loop?

[Star Strider]

O.K. With that, I get one positive root of 49.123...

[Vidhan Malik]

49 is actually the efficiency but at the end of the day it means that we are getting a positive value of R. So with the second code, I want to vary the value of To6 and graph the impact it has on the Eff term.

It seems that I may have made some error in between when trying to make a loop.

[Star Strider]

How do you want to vary ‘To6’?

It would be much easier to do this without the Symbolic Math Toolbox, since it’s intended for one-offs, not iterations, and is slow.

You might even be able to vectorise it, although considering the number of computations, that doesn’t seem likely.

Getting the positive root would be relatively easy with the fzero function.

It would help if I knew what you want to do.

[Vidhan Malik]

I want to vary To6 from the value of 1000 to 1800 by 1 while To7 = 2/3*To6. So with that, what would be the most efficient way of plotting and calculating "Eff" by To6?

[Star Strider]

...Taking a break...Back in the morning...

[Vidhan Malik]

haha I followed suite!

[Vidhan Malik]

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]

Sorry for such an unclear code

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

That line of code is in line 38 and this is the equation that we need to solve for the positive root of R. The rest of the equations before this line are just intermediate steps.

[Vidhan Malik]

Sorry I made another mistake, To6 has to start at 1700, any value below that gives a value of R greater than 1. Which in context to the problem, doesn't make sense.

[Vidhan Malik]

I get a value now, it just happens to be way too high to be correct. But is this code correct in terms of setting up and solving an anonymous function? Also how would you solve for the positive root in this case?

      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(1), To6(2),R);
      R = fsolve(@(R) Eqn(R), 0.9);

[Star Strider]

I’m lost as to what changes you’ve made.

Please post your complete current code.

[Vidhan Malik]

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,2) = 1699;

for i = 2:1:202

      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 = @(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(1), To6(2),R);
      R = fsolve(@(R) Eqn(R), 0.9);

      Wnet(i,2) = (1+R)+fo/(1+R) - 0.12*(1+R)*Yc*Rg*(To6(i,2)-To7) - (1+R)*Ym*Rg*(To4-T32); %Step 7

      Qin(i,2) = (1+R+fo/(1+R)-0.12*(1+R))*(Yc*Rg*To6(i,2)-Ym*Rg*To5(i,2));

      Eff(i,2) = Wnet(i,2)/Qin(i,2)*100; %Step 8

end

disp(R);

[Star Strider]

I don’t understand what you’re doing here:

Eqn = @(R) E(To5(1), To6(2),R);

You’re not changing either ‘To5’ or ‘To6’, so ‘R’ will never change.

See if this does what you want:

        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(i), To6(i),R);
        R = fsolve(Eqn, 0.9);
        Rv(i) = R;                                        % Record ‘R’

(I added ‘Rv’ because I want to see what it is and how it changes between iterations. Delete it if you don’t want to.)

[Vidhan Malik]

For the equation that you pointed out that didn't make sense to you, I was just trying out using anonymous functions, I don't know what I'm doing.

I get the same value as I did before interestingly. So in the loop, the value of To5 and To6 are changing, and so the value of R should also be changing.

[Star Strider]

The value of ‘R’ (saved in ‘Rv’) did change when I ran it. I don’t know what you’re doing, so I don’t know how much it should change.

Plotting ‘Eff’ as a function of ‘To6’ shows an initial sharp drop, then a linear or close-to-linear increase from about -39 to about -36 with ‘To6’ going from 1700 to 1900:

figure(1)
plot(To6, Eff)
grid

[Vidhan Malik]

In the solve equation, how can we specify which root should be used?

[Star Strider]

The roots fsolve returns depends on the intital guess you give it. It will return the ‘nearest’ root to that estimate, ‘near’ taking into account the shape (derivative) of the curve at the initial estimate. This also applies to complex roots. If you give it a complex initial estimate, it will search the complex numbers for a root. If you give it a real initial estimate, it will search the real numbers for a root. It is frequently necessary to experiment.

One way to see where the potential roots may be is to create a vector of values for ‘R’ (here), and then plot the function, here ‘E(To5(i), To6(i),R)’ to see where to zero-crossings are. Unfortunately, with three independent variables, it would be difficult to analyse it graphically with respect to all three.

[Vidhan Malik]

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.

[Vidhan Malik]

For some reason, the values of R are coming out to be different using the code with the code, as opposed to the code without the code. And as far as I can tell, the only difference is that in the code without the loop I use vpasolve and in the code with the loop I am using fsolve.

[Vidhan Malik]

I'm not really too sure of what I did but its really starting to look like what it's supposed to!

Just have to figure out why there is that initial steep increase, but I think that is mainly due to the different value of R that I get as compared to using vpasolve

[Star Strider]

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]

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!

[Vidhan Malik]

To put all this in perspective, here is the system I am analyzing!

And this code just deals with the microturbine subsystem, now onwards to the rest of the system! :) As always, your time and effort are much appreciated! (They should really be paying you for all of your help)

[Star Strider]

My pleasure!

I actually do this for fun!

I had no idea what you were doing. It looks like a jet engine.

Did you see my one-line solution for the positive root for ‘R’ in my previous Comment? You no longer need any symbolic or numeric solution step for it any more, since the root is already solved for symbolically. Your code should also be faster.

[Vidhan Malik]

no I didn't use it, but now that I am going through it again, it definitely seems like the better way to do this. I once again actually didn't know that you could solve for a single variable in a equation via matlab, which definitely makes things easier!

How did you use the symbolic toolbox to solve equation E for R?

Also the system is actually for terrestrial applications such as in a powerplant!

[Vidhan Malik]

Just saw that you already showed me, thank you as always!

[Vidhan Malik]

How did you find what the R_fcn function in it's variable form and it's numeric form?

[Star Strider]

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]

Since R has To6 in it, the value of it actually changes with each iteration so you do need to solve for it in each iteration. And then it gets complicated I would think when you would have to assign the matrix value of To6(i) to the variable To6 and use that in the "Rss" equation. But I definitely prefer this way of doing it, makes more sense to do it this way and is more efficient!

In that, how do you make R a function where inputs are To6 and To5?

[Star Strider]

Looking at the ‘R’ assignment I posted, it will pick up ‘Qr’,‘Rg’,‘To6’,‘To7’,‘Yc’,‘Ym’ from the workspace. It is not a direct function of ‘To5’, according to the simplification. It inherits that value through ‘To7’.

The root calculation for ‘R’ will be the same regardless of the values in your workspace, just as it would be for any quadratic equation with varying coefficients. The value will of course change depending on the values it uses. That is the reason I chose to do a symbolic calculation that you can then use directly rather than having to recalculate the roots either symbolically or numerically each time. The positive root of ‘E’ with respect to ‘R’ is what you want to calculate, so it is best to do it once, since the structure of ‘E’ does not change, and so the solution with respect to ‘R’ does not change.

[Vidhan Malik]

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]

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.