vectorization this for loop

2017 10 24
1 个回答
更新时间：2017 10 24
8 次查看（30 天）

0 个投票

Hi all I have some problem when I run this script. It has two scripts, and the main script is running very slow. I tried to vectorize all the loop I can and the for loop outside is the only loop I can't vectorize.
script one:function [Vstep]=iteratingvol(~)
%--------------engine speces---------------- b=3.7; %bore s=3.9; %stroke l=5.7; %connecting rod length a=(1/2)*s; %crank radius
%x is piston location; at TDC, x=0; at BDC, x=s theta=-10:1:90;
Vstep=(a+l-(a.*cosd(theta)+(l^2-a^2.*(sind(theta)).^2).^0.5))*(b/2)^2*pi;
Vstep=fliplr(Vstep);
main script:P2=177.8034; T2=1185.4; V2=6; n2=0.000048544;
kcomb=1.4; kcomp=1.4; mrc=1.05; dT=5035; %-----------------------slice model set up--------------------------------- slicenm=100;
P=repmat(P2,1,slicenm); T=repmat(T2,1,slicenm);
V=repmat((V2/slicenm),1,slicenm);
n=repmat((n2/slicenm),1,slicenm);
Vstep=iteratingvol();
%------------Iteration terms----------- deltap=1; %incremental pressure for iterations
for i=1:slicenm
%selected slice combust and increase temperature T(1,i)=T2+dT; n(1,i)=(n2/slicenm)*mrc; P(1,i)=n(1,i)*T(1,i)*R*12/V(1,i);
%expand to the same pressure as unburn slices if i==slicenm V(1,i)=(P(1,i)*V(1,i)^kcomb/P(1,i-1))^(1/kcomb); P(1,i)=P(1,i-1); else V(1,i)=(P(1,i)*V(1,i)^kcomb/P(1,i+1))^(1/kcomb); P(1,i)=P(1,i+1); end
T(1,i)=P(1,i)*V(1,i)/(n(1,i)*R*T(1,i)*12);
%compress slices back to orginal volume before combustion
while sum(V(1,:)) > Vstep1(1,i)
V=((P.*V.^kcomp)./(P+deltap)).^(1/kcomp);
P=P+deltap;
T=P.*V./(n.*R.*12);
end
Tr=T;
Vr=V;
Pr=P;
%compress slices into smaller volume due to piston motion
while sum(V(1,:)) > Vstep1(1,i+1)
V=((P.*V.^kcomp)./(P+deltap)).^(1/kcomp);
P=P+deltap;
T=P.*V./(n.*R.*12);
end
end
The script one is used to define term 'Vstep' which is going to use in main script. If 'Vstep' is defined purely by typing each element one by one, then the main script is running fine. But if term 'Vstep' is calculated by another script, then the main script is running super slow.
Let me know if you guys have any thought, I will appreciate that very much.

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

回答（1 个）

shanghua chen
shanghua chen 2017-10-24

0 个投票

I may retype the scripts.
script one:
function [Vstep]=iteratingvol(~)
%--------------engine speces---------------- b=3.7; %bore s=3.9; %stroke l=5.7; %connecting rod length a=(1/2)*s; %crank radius
%x is piston location; at TDC, x=0; at BDC, x=s theta=-10:1:90;
Vstep=(a+l-(a.*cosd(theta)+(l^2-a^2.*(sind(theta)).^2).^0.5))*(b/2)^2*pi;
Vstep=fliplr(Vstep);
main script:
P2=177.8034; T2=1185.4; V2=6; n2=0.000048544;
kcomb=1.4; kcomp=1.4; mrc=1.05; dT=5035; %-----------------------slice model set up--------------------------------- slicenm=100;
P=repmat(P2,1,slicenm); T=repmat(T2,1,slicenm);
V=repmat((V2/slicenm),1,slicenm);
n=repmat((n2/slicenm),1,slicenm);
Vstep=iteratingvol();
%------------Iteration terms----------- deltap=1; %incremental pressure for iterations
for i=1:slicenm
%selected slice combust and increase temperature T(1,i)=T2+dT; n(1,i)=(n2/slicenm)*mrc; P(1,i)=n(1,i)*T(1,i)*R*12/V(1,i);
%expand to the same pressure as unburn slices if i==slicenm V(1,i)=(P(1,i)*V(1,i)^kcomb/P(1,i-1))^(1/kcomb); P(1,i)=P(1,i-1); else V(1,i)=(P(1,i)*V(1,i)^kcomb/P(1,i+1))^(1/kcomb); P(1,i)=P(1,i+1); end
T(1,i)=P(1,i)*V(1,i)/(n(1,i)*R*T(1,i)*12);
%compress slices back to orginal volume before combustion
while sum(V(1,:)) > Vstep1(1,i)
V=((P.*V.^kcomp)./(P+deltap)).^(1/kcomp);
P=P+deltap;
T=P.*V./(n.*R.*12);
end
Tr=T;
Vr=V;
Pr=P;
%compress slices into smaller volume due to piston motion
while sum(V(1,:)) > Vstep1(1,i+1)
V=((P.*V.^kcomp)./(P+deltap)).^(1/kcomp);
P=P+deltap;
T=P.*V./(n.*R.*12);
end
end

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Jan
Jan 2017-10-24
No, do not "retype"it and post in again in the section for answers. Remove this pseudo answer and use the "{} Code" button to format your code. See http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup#answer_18099

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Spline Postprocessing 的更多信息

标签

评论：

Jan
Jan
2017-10-24

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by