How do I store Iteration Results
    3 views (last 30 days)
  
       Show older comments
    
Hi, Please I need to get values for Z, Za, V, Va at every value of T. The code is as follows:
z=[0.2 0.2 0.2 0.4];
W_PR=0.245;
C=4;
omega=[0.344 0.467 0.578 0.789];
Tc=[600 700 500 570];
Pc=[50 70 58 76];
R=8.314; 
P=20; 
Z_store = ones(C,C,5);
V_store = ones(5,C);
K_PR=[1.546e-2, 0.456, 1.432e2, 14.32];
iter = 0;
for k=40:10:80
T(k)=273.15+k;
  for c=1:C
    x_PR(1,c)=z(c)/(1+W_PR*(K_PR(c)-1));
    x_PR(2,c)=K_PR(c)*x_PR(1,c);  
    kappa_PR=0.37464+1.54226.*omega(c)-0.26992.*omega(c).^2;
    alpha_PR=(1+kappa_PR.*(1-sqrt(T(k)./Tc(c)))).^2;
    a_PR(c,c)=0.45724.*R.^2.*Tc(c).^2./Pc(c).*alpha_PR;
    b_PR(c)=0.07780*R.*Tc(c)./Pc(c);
end
for c=2:C
    for n=1:(c-1)
        a_PR(c,n)=sqrt(a_PR(c,c).*a_PR(n,n));
        a_PR(n,c)=a_PR(c,n);
    end
end
for c=1:C
    A_PR(c,c)=a_PR(c,c).*P./(R.*T(k)).^2;
    B_PR(c)=b_PR(c).*P./(R.*T(k)); 
    Z(c,c)=A_PR(c,c)./5;
    V(c)=B_PR(c).*6;
end
for c=1:C
      Aa_PR(c,c)=a_PR(c,c).*P./(R.*T(k)).^3;
      Ba_PR(c)=b_PR(c).*P./(R.*T(k));
      Za(c,c)=A_PR(c,c)./8;
      Va(c)=B_PR(c).*9;
   end       
end
Thanks
Isa
[EDITED, Jan, code formatted]
3 Comments
  Image Analyst
      
      
 on 3 Sep 2012
				You should have done that in your first post, not copied it as a comment.
Answers (3)
  Azzi Abdelmalek
      
      
 on 3 Sep 2012
         z=[0.2 0.2 0.2 0.4]; W_PR=0.245; C=4; 
omega=[0.344 0.467 0.578 0.789]; 
Tc=[600 700 500 570]; Pc=[50 70 58 76]; R=8.314; P=20;
Z_store = ones(C,C,5); V_store = ones(5,C);
K_PR=[1.546e-2, 0.456, 1.432e2, 14.32]; iter = 0;
 for k=40:10:80
T(k)=273.15+k; 
  for c=1:C 
      x_PR(1,c)=z(c)/(1+W_PR*(K_PR(c)-1));
      x_PR(2,c)=K_PR(c)*x_PR(1,c); 
      kappa_PR=0.37464+1.54226.*omega(c)-0.26992.*omega(c).^2; 
      alpha_PR=(1+kappa_PR.*(1-sqrt(T(k)./Tc(c)))).^2;
      a_PR(c,c)=0.45724.*R.^2.*Tc(c).^2./Pc(c).*alpha_PR;
      b_PR(c)=0.07780*R.*Tc(c)./Pc(c);
  end
  for c=2:C
      for n=1:(c-1)
          a_PR(c,n)=sqrt(a_PR(c,c).*a_PR(n,n)); 
          a_PR(n,c)=a_PR(c,n);
      end
  end
  for c=1:C 
      A_PR(c,c)=a_PR(c,c).*P./(R.*T(k)).^2;
      B_PR(c)=b_PR(c).*P./(R.*T(k));
      Z(c,c)=A_PR(c,c)./5; V(c)=B_PR(c).*6; 
  end
  for c=1:C 
      Aa_PR(c,c)=a_PR(c,c).*P./(R.*T(k)).^3;
      Ba_PR(c)=b_PR(c).*P./(R.*T(k)); 
      Za(c,c)=A_PR(c,c)./8; Va(c)=B_PR(c).*9;
  end
  result{k,1}=Z;result{k,2}= Za; result{k,3}= V;result{k,4}=  Va
end
5 Comments
  Jan
      
      
 on 11 Sep 2012
				@Isa Isa, I still do not get the problem. Does the code you have posted calculate Z, Za, V and Va correctly? If not, why? Did you set some breakpoints in the code to inspect what's going on? What kind of help do you expect from us now?
  Jan
      
      
 on 11 Sep 2012
        This looks strange:
for c=2:C
  for n=1:(c-1)
    a_PR(c,n)=sqrt(a_PR(c,c).*a_PR(n,n));
    a_PR(n,c)=a_PR(c,n);
  end
end
The values of a_Pr are overwritten by its own transposed. Please check if this is really doing what you want.
I do not see a chance that we can find the bugs of your program, because we cannot know what you are wanting to achieve. I suggest to use the debugger to find out, what's going on.
3 Comments
  Jan
      
      
 on 11 Sep 2012
				Hi Isa, and what I ask is, if the posted code is correct, because it is a frequently implemented bug (FIB) to overwrite values in a matrix at an inplace-transpose.
  Jan
      
      
 on 11 Sep 2012
        
      Edited: Jan
      
      
 on 12 Sep 2012
  
      Your code has this problem:
for k = 40:10:80
  T(k) = 273.15+k; 
  ...
  result{k,1}=Z;result{k,2}= Za; result{k,3}= V;result{k,4}=  Va
end
Now the first value written to the index k=40! Improvement:
kList = 40:10:80;
for ik = 1:length(kList)
  k = kList(ik);    % [EDITED, "iK" -> "ik"]
  T(ik) = 273.15 + k;   % Not T(k)!
  ...
  result{ik, 1} = Z; 
  result{ik, 2} = Za;
  result{ik, 3} = V;
  result{ik, 4} = Va
end
BTW. your code would be cleaner and faster, if it is vectorized:
% Instead of: for c=1:C 
x_PR       = z ./ (1 + W_PR .* (K_PR.' - 1));
x_PR(2, :) = K_PR .* x_PR.'; 
kappa_PR   = 0.37464+1.54226.*omega-0.26992.*omega.^2; 
alpha_PR   = (1+kappa_PR.*(1-sqrt(T(ik)./Tc))).^2;
a_PR       = diag(0.45724.*R.^2.*Tc .^ 2 ./ Pc .* alpha_PR);
b_PR       = 0.07780 * R .* Tc ./ Pc;
3 Comments
  Jan
      
      
 on 12 Sep 2012
				The error message contain "iK" with an uppercase "K", while the variable is called "ik" with a lowercase "k". Is it really too complicated to fix this by your own?
See Also
Categories
				Find more on Data Type Identification in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


