在 R 中实现快速数值计算

Question

我试图在 R 中进行大量计算。十八个小时过去了，但我的 RStudio 似乎继续工作。我不确定我是否可以用不同的方式编写脚本以使其更快。我试图在 50000 x 350 矩阵上实现 Crank–Nicolson type method，如下所示：

#defining the 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 the big matrix (all filled with zeros)
m <- as.data.frame(matrix(0, t/dt+1, z/dz+2)) #extra columns/rows for boundary conditions

#fill the first and last columns with constant boundary values
m[,1]<-400
m[,length(m)]<-0

#implement the 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)
  }}

有没有办法知道 R 实现它需要多长时间？有没有更好的方法来构建数值计算？在这一点上，我觉得 excel 本来可以更快！！

Answer 1

只需进行一些简单的优化就很有帮助。您代码的原始版本代码在我的笔记本电脑上需要大约 5 天的时间。使用矩阵并仅计算一次在循环中重复使用的值，我们将其缩短到大约 7 分钟

想想像

这样的混乱结构

m[i,][2:length(m)-1][[j]]

这相当于

m[[i, j]]

哪个会更快（也更容易理解）。进行此更改进一步将运行时间减少了 2 倍以上，大约 3 分钟

综合起来我们有

dt<-1
t<-50000
dz<-0.0075
z<-350*dz

#velocity & diffusion 
v<-2/(24*60*60)
D<-0.02475/(24*60*60) 

#make the big matrix (all filled with zeros)
m <- (matrix(0, t/dt+1, z/dz+2)) #extra columns/rows for boundary conditions

# cache a few values that get reused many times 
NC = NCOL(m)
NR = NROW(m)
C1 = D*dt / dz^2
C2 = v*dt / (2*dz)

#fill the first and last columns with constant boundary values
m[,1]<-400
m[,NC]<-0

#implement the 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)
    }
  }

如果你需要比这更快，你可以尝试更多的优化。例如，请参阅 here 以了解对同一元素进行索引的不同方式如何具有非常不同的执行时间。一般来说，最好先引用列，然后再引用行。

如果您可以在 R 中进行的所有优化都不足以满足您的速度要求，那么您可以改为在 RCpp 中实现循环。

在 R 中实现快速数值计算

Implementing fast numerical calculations in R

r

numerical-methods

fluid-dynamics