빠른 수치 계산을 구현 R

Question

R에서 광범위한 계산을 시도했습니다. 18 시간이 지났지 만 RStudio는 계속 작동하는 것 같습니다. 스크립트를 다른 방식으로 작성하여 더 빠르게 만들 수 있는지 확신 할 수 없습니다.

#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 그것을 구현하는 것이 얼마나 걸릴지 알 수있는 방법이 있나요 : 나는 아래와 같이 Crank–Nicolson type method (350)에 의해 50,000 이상 매트릭스를 구현하려고했다? 수치 계산을 구성하는 더 좋은 방법이 있습니까? 이 시점에서, 나는 더 빠른 것처럼 느껴졌습니다!

Answer 1

간단한 최적화를하면 도움이됩니다. 귀하의 코드의 원래 버전 코드는 내 노트북에서 ~ 5 일이 걸릴 것입니다. 매트릭스를 사용하여 루프에서 재사용 한 번만 값을 계산, 우리는 주변 칠분

을이 쓰러 뜨리 그리고이

m[[i, j]] 

에 해당

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

처럼 지저분한 구조에 대해 생각

은 더 빠릅니다 (이해하기가 훨씬 쉽습니다). 2 이상의 다른 요인에 의해 런타임이 약 삼분

에 더이 변경을 감소 만들기 당신이, 당신은 좀 더 시도 할 수도 이보다 빨리 이동해야하는 경우 우리가

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에 루프를 구현할 수 있습니다.