Help to speed up the code
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
Kevin Isakayoga
il 19 Ott 2020
Commentato: Durganshu
il 26 Ott 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 Commenti
Risposta accettata
Durganshu
il 19 Ott 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 Commenti
Durganshu
il 26 Ott 2020
Sorry for a late reply. I'm actually quite busy. i'll try to speed up your code when I'm free.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!