how to write a for loop for the following task?

Known:

Kp(1)=a;

Then:

x = roots([Kp-1 7.56*Kp -18.12*Kp 9.56*Kp]); (we only need 0<x<1). 

T = roots([-0.001273*x+0.00365 0.544*x+44.3191 283338.4*x-407295]); 
mu_CO2 = -394088+44.3191*(T-298-T*ln(T/298))-0.0073/2*(T-298)^2-213.984*T; 
mu_CO = -110700+29.6127*(T-298-T*ln(T/298))-0.00301/2*(T-298)^2-197.81*T; 
mu_O2 = 30.5041*(T-298-T*ln(T/298))-0.00349/2*(T-298)^2-205.31*T; 
Kp_new = exp(-(mu_O2+2*mu_CO-2*mu_CO2)/(8.314*T));

Finally, I need to use Kp_new to calculate new x, new T, new mu. The iteration number is supposed to be 20.

I will appreciate if someone can help me on this problem!

[Kaushik Lakshminarasimhan]

Instead of calling your variable as Kp_new, just call it Kp. And put the whole code starting from ( x = ...) to (Kp = ...) inside a for loop.

Look here for how to create 20 loops: https://www.mathworks.com/help/matlab/matlab_prog/loop-control-statements.html

[Ivy Shen]

Thank you! Shall I use i for all these variables everywhere? such as x(i),T(i),mu_co2(i)....

[Ivy Shen]

I tried, but it said that the subscripted assignment dimension mismatch for x.

[Kaushik Lakshminarasimhan]

You can just overwrite the variables. No need to store intermediate values as x(i), T(i) etc. unless you want to look at the values in intermediate iterations.

You're probably getting the error because you have multiple roots. You can retain the root between 0 and 1 using:

x = x(x>0 & x<1)

You many have to do the same thing for T if there are multiple roots.

[Ivy Shen]

I am still not sure how can I make the for loop work. Where should I put the iteration index? If I just wrote from "x=.." to "Kp=..," I only get one value per time. How can I make it automatically calculate 20 times with the new Kp values?

[KSSV]

copy the whole code here...so that we can help you.

[Ivy Shen]

Kp(1) = 7.75*10^(-4);

for i = 2:21

      x_root = roots([Kp(i-1)-1 7.56*Kp(i-1) -18.12*Kp(i-1) 9.56*Kp(i-1)]);
      x = x_root(x_root>0 & x_root<1);
      T_root = roots([-0.001273*x+0.00365 0.544*x+44.3191 283338.4*x-407295]);
      T = T_root(T_root>0 & T_root<1);
      mu_CO2 = -394088+44.3191*(T-298-T.*log(T/298))-0.0073/2*(T-298).^2-213.984*T;
      mu_CO = -110700+29.6127*(T-298-T.*log(T/298))-0.00301/2*(T-298).^2-197.81*T;
      mu_O2 = 30.5041*(T-298-T.*log(T/298))-0.00349/2*(T-298).^2-205.31*T;
      Kp = exp(-(mu_O2+2*mu_CO-2*mu_CO2)/(8.314*T));
  end

[KSSV]

   T = T_root(T_root>0 & T_root<1);

This line giving you empty output....

[KALYAN ACHARJYA]

It seems problem arises in the following:

x_root = roots([Kp(i-1)-1 7.56*Kp(i-1) -18.12*Kp(i-1) 9.56*Kp(i-1)]);

x = x_root(x_root>0 & x_root<1);

Check the following what does this statement mean-

[Walter Roberson]

On the first iteration, when i = 2, then x_root generates 3 roots, two of which are complex. The selection of only some of those narrows it down to one value for x. The root() for T_root then results in two roots. Neither of them is in the desired range, so T comes out empty. That gives you empty variables for the next several lines. You then write that empty content over the entire Kp array instead of just Kp(i) . But you cannot just assign into Kp(i) because you have only emptiness to write in.

[Ivy Shen]

Thank you! I revised the code again as following. Actually, what I want to get from the code is the x value when Kp does not change much (in other words, when the solution converges). Do you think this code works correctly?

Kp(1) = 7.75*10^(-4);

for i = 2:21

      x_root = roots([Kp(i-1)-1 7.56*Kp(i-1) -18.12*Kp(i-1) 9.56*Kp(i-1)]);
      x(i-1) = x_root(x_root>0 & x_root<1);
      T_root = roots([-0.001273*x(i-1)+0.00365 0.544*x(i-1)+44.3191 283338.4*x(i-1)-407295]);
      T = T_root(T_root>0);
      mu_CO2 = -394088+44.3191*(T-298-T.*log(T/298))-0.0073/2*(T-298).^2-213.984*T;
      mu_CO = -110700+29.6127*(T-298-T.*log(T/298))-0.00301/2*(T-298).^2-197.81*T;
      mu_O2 = 30.5041*(T-298-T.*log(T/298))-0.00349/2*(T-298).^2-205.31*T;
      Kp(i) = exp(-(mu_O2+2*mu_CO-2*mu_CO2)/(8.314*T));
  end

[Walter Roberson]

You should change

x(i-1) = x_root(x_root>0 & x_root<1);

to

x(i-1) = x_root(imag(x_root) == 0 & x_root>0 & x_root<1);

[Walter Roberson]

Last night I tried finding the range of Kp values that left x_root in the range 0 to 1. It turned out that for 0 exactly, the solutions were 0 and 3775/7744 with that second value being an isolated real-valued spot in the middle of complex-valued Kp solutions. It also turned out that due to discontinuities, there were no finite values of Kp that made any of the x_root values exactly 1, but infinite Kp made it one (that is, the values were greater than one for finite values but converged to one at infinity.)