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); endSr1(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
Best Answer