Help to speed up the code
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Kevin Isakayoga
am 19 Okt. 2020
Kommentiert: Durganshu
am 26 Okt. 2020
Hello everyone, I need to speed up this simulation code because it took around several hours, can somebody help me how to speed it up?
pc=3.16;
T=293;
rwc3s=0.47;
rwc2s=0.27;
rwc3a=0.10;
rwc4af=0.09;
mc=0.4;
mw=0.157;
ms=0.658;
mg=1.129;
wc=mw/mc;
vc=1/((wc*pc)+1);
ro=(3*vc/(4*pi))^1/3;
%% Proposed Model
yg=0.25;
yw=0.15;
RH=0.88;
b1=1231;
b2=7579;
B293=0.3*10^-10;
C=5*10^7;
ER=5364;
De293=((rwc3s*100)^2.024)*(3.2*10^-14);
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);
v=2.06;
cw=((RH-0.75)/0.25)^3;
pw=1;
%De=De293*exp(-b2*(1/T-1/293));
B=B293*exp(-b1*(1/T-1/293));
kr=kr293*exp(-ER*(1/T-1/293));
%% Determine the day
hr=15;
da=24*hr;
dt=1;
Nn=(da)*(1/dt);
time=1:Nn;
time2=1:Nn+1;
%% Simulation MC
n=100000; %this is the number of trial
rlower=0.5;
rupper=125;
r_initial=rlower+(rupper-rlower)*rand(1,n);
m=length(r_initial);
n1=length(time);
%% Pre-allocation speed
n2=length(time2);
Sr=zeros(1,n1);
cst=zeros(1,n1);
alfa=zeros(1,n1);
kd=zeros(1,n1);
De=zeros(1,n1);
rt_inc=zeros(1,n1);
Rt_inc=zeros(1,n1);
Alfa=zeros(m,n1);
Rt_inc1=zeros(m,n1);
rt_inc1=zeros(m,n1);
ro_init1=zeros(m,1);
alfag1=zeros(m,n1);
Vdp1=zeros(m,1);
vdp1=zeros(m,1);
for i = 1:m
ro_init=r_initial(i)*0.001;
L=((4*pi*(wc*pc/pw+1)/3)^(1/3))*ro_init;
for AA=1:Nn
if (AA==1)
rt_inc(1)=ro_init(1)-0.00001; %boundary condition
Rt_inc(1)=ro_init(1)+0.00001; %boundary condition
else
rt_inc(1)=ro_init(1)-0.00001;
Rt_inc(1)=ro_init(1)+0.00001;
end
end
for BB=1:Nn
if Rt_inc(BB)/L<=ro
Sr(BB)=0;
elseif (Rt_inc(BB)/L>=ro)&&(Rt_inc(BB)/L<0.5)
Sr(BB)=4*pi*(Rt_inc(BB)/L)^2;
elseif (0.5<=Rt_inc(BB)/L)&&(Rt_inc(BB)/L<0.5*(2^0.5))
Sr(BB)=(4*pi*(Rt_inc(BB)/L)^2)-(6*pi*(1-(0.5/(Rt_inc(BB)/L))));
elseif (Rt_inc(BB)/L>=0.5*(2^0.5)) && (Rt_inc(BB)/L<0.5*(3^0.5))
syms y x
fun = @(y,x) 8*(Rt_inc(BB)/L)./(sqrt((Rt_inc(BB)/L)^2-(x.^2)-(y.^2)));
ymin=sqrt((Rt_inc(BB)/L)^2-0.5);
xmin=@(x) sqrt((Rt_inc(BB)/L)^2-0.25-x.^2);
Sr(BB)=integral2(fun,ymin,0.5,xmin,0.5);
else
Sr(BB)=0;
end
cst(BB)=Sr(BB)/(4*pi*Rt_inc(BB)^2);
alfa(BB)=1-(rt_inc(BB)/ro_init)^3;
kd(BB)=(B/(alfa(BB)^1.5))+C*(Rt_inc(BB)-ro_init)^4;
De(BB)=De293*(log(1/alfa(BB)))^1.5;
rt_inc(BB+1)=rt_inc(BB)-(dt*((pw*cst(BB)*cw)/((yw+yg)*pc*rt_inc(BB)^2))*1/((1/(kd(BB)*rt_inc(BB)^2))+(((1/rt_inc(BB))-(1/Rt_inc(BB)))/De(BB))+(1/(kr*rt_inc(BB)^2))));
Rt_inc(BB+1)=((v-1)*((rt_inc(BB)^2)/(Rt_inc(BB)^2))*((rt_inc(BB)-rt_inc(BB+1))/dt))*dt+Rt_inc(BB);
end
Sr1(i,:)=Sr;
Alfa(i,:) = alfa ;
Rt_inc1(i,:) = Rt_inc(1:n1);
rt_inc1(i,:) = rt_inc(1:n1);
rt_inc2(i,:) = rt_inc;
ro_init1(i,:)= ro_init;
b=0.0517;
n=1.0145;% 0.6;
Vdp=(1-exp(-b*((2*ro_init*1000)^n)));
vdp=2*b*n*exp(-b*(2*ro_init*1000)^n)*(2*ro_init*1000)^(n-1);
N=6*(10^12)*mc*vdp/(pc*pi*(2*ro_init*1000)^3);
alfag=alfa*vdp;
alfag1(i,:)=alfag;
Vdp1(i,:)=Vdp;
vdp1(i,:)=vdp;
N1(i,:)=N;
end
Thankyou in advance.
Best Regads,
Kevin
0 Kommentare
Akzeptierte Antwort
Durganshu
am 19 Okt. 2020
Well, the code inside the loops can be edited and vectorized to obtain faster results. I would suggest you go through this documentation on vectorization that may help you:
2 Kommentare
Durganshu
am 26 Okt. 2020
Sorry for a late reply. I'm actually quite busy. i'll try to speed up your code when I'm free.
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!