How do I store Iteration Results
古いコメントを表示
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 件のコメント
Jan
2012 年 9 月 3 日
Please format your code. The topic code formatting has been discussed in this forum such frequently, that I'm starting to get strange and wild emotions when I see the next ugly posting.
Image Analyst
2012 年 9 月 3 日
You should have done that in your first post, not copied it as a comment.
回答 (3 件)
Azzi Abdelmalek
2012 年 9 月 3 日
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 件のコメント
Isa Isa
2012 年 9 月 3 日
Azzi Abdelmalek
2012 年 9 月 3 日
i agree, there are many []. but that's the result of your code
Azzi Abdelmalek
2012 年 9 月 3 日
but there are some results non empty
result{end,1}
result{end,2}
result{end,3}
result{end,3}
Isa Isa
2012 年 9 月 11 日
Jan
2012 年 9 月 11 日
@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
2012 年 9 月 11 日
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.
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;
カテゴリ
ヘルプ センター および File Exchange で Argument Definitions についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!