Implementing fast numerical calculations in R -
i trying extensive computation in r. eighteen hours have passed rstudio seems continue work. i'm not sure if have written script in different way make faster. trying implement crank–nicolson type method on 50000 350 matrix shown below:
#defining discretization of cells dt<-1 t<-50000 dz<-0.0075 z<-350*dz #velocity & diffusion v<-2/(24*60*60) d<-0.02475/(24*60*60) #make big matrix (all filled zeros) m <- as.data.frame(matrix(0, t/dt+1, z/dz+2)) #extra columns/rows boundary conditions #fill first , last columns constant boundary values m[,1]<-400 m[,length(m)]<-0 #implement calculation for(j in 2:(length(m[1,])-1)){ for(i in 2:length(m[[1]])){ m[i,][2:length(m)-1][[j]]<-m[i-1,][[j]]+ d*dt*(m[i-1,][[j+1]]-2*m[i-1,][[j]]+m[i-1,][[j-1]])/(dz^2)- v*dt*(m[i-1,][[j+1]]-m[i-1,][[j-1]])/(2*dz) }}
is there way know how long take r implement it? there better way of constructing numerical calculation? @ point, feel excel have been faster!!
just making few simple optimisations helps here. original version code of code take ~ 5 days on laptop. using matrix , calculating once values reused in loop, bring down around 7 minutes
and think messy constructions
m[i,][2:length(m)-1][[j]]
this equivalent to
m[[i, j]]
which faster (as easier understand). making change further reduces runtime factor of on 2, around 3 minutes
putting have
dt<-1 t<-50000 dz<-0.0075 z<-350*dz #velocity & diffusion v<-2/(24*60*60) d<-0.02475/(24*60*60) #make big matrix (all filled zeros) m <- (matrix(0, t/dt+1, z/dz+2)) #extra columns/rows boundary conditions # cache few values reused many times nc = ncol(m) nr = nrow(m) c1 = d*dt / dz^2 c2 = v*dt / (2*dz) #fill first , last columns constant boundary values m[,1]<-400 m[,nc]<-0 #implement calculation for(j in 2:(nc-1)){ for(i in 2:nr){ ma = m[i-1,] ma.1 = ma[[j+1]] ma.2 = ma[[j-1]] m[[i,j]] <- ma[[j]] + c1*(ma.1 - 2*ma[[j]] + ma.2) - c2*(ma.1 - ma.2) } }
if need go faster this, can try out more optimisations. example see here how different ways of indexing same element can have different execution times. in general better refer column first, row.
if optimisations can in r not enough speed requirements, might implement loop in rcpp instead.
Comments
Post a Comment