R: Vectorize Finite Difference Equations

Question

I'm trying to move some Fortran code to R for finite differences related to chemical kinetics.

Sample Fortran loop:

DOUBLE PRECISION, DIMENSION (2000,2) :: data=0.0
DOUBLE PRECISION :: k1=5.0, k2=20.0, dt=0.0005
DO i=2, 2000
  data(i,1) = data(i-1,1) + data(i-1,1)*(-k1)*dt
  data(i,2) = data(i-1,2) + ( data(i-1,1)*k1*dt - data(i-1,2)*k2*dt )
  ...
END DO

The analogous R code:

k1=5
k2=20
dt=0.0005
data=data.frame(cbind(c(500,rep(0,1999)),rep(0,2000)))
a.fun=function(y){
     y2=y-k1*y*dt
     return(y2)
 }
apply(data,2,a.fun)

This overwrites my first value in the dataframe and leaves zeros elsewhere. I'd like to run this vectorized and not using a for loop since they are so slow in R. Also, my function only calculates the first column so far. I can't get the second column working until I get the syntax right on the first.

Answer 1

Its not necessarily true that R is bad at loops. It very much depends on what you are doing. Using k1, k2, dt and data from the question (i.e. the four lines beginning with k1=5) and formulating the problem in terms of an iterated matrix, the loop in the last line below returns nearly instantaneously on my PC:

z <- as.matrix(data)
m <- matrix(c(1-k1*dt, k1*dt, 0, 1-k2*dt), 2)

for(i in 2:nrow(z)) z[i, ] <- m %*% z[i-1, ]

(You could also try storing the vectors in columns of z rather than rows since R stores matrices by column.)

Here is the first bit of the result:

> head(z)
           X1       X2
[1,] 500.0000 0.000000
[2,] 498.7500 1.250000
[3,] 497.5031 2.484375
[4,] 496.2594 3.703289
[5,] 495.0187 4.906905
[6,] 493.7812 6.095382

Answer 2

May be this can help.

I think you need to have the initial condition for data[1,2]. I assumed both data[1,1] as 500 and data[1,2 as 0 at the initial condition.

The code goes like this:

> ## Define two vectors x and y
> x <- seq(from=0,length=2000,by=0)
> y <- seq(from=0,length=2000,by=0)
> 
> ## Constants
> k1 = 5.0
> dt = 0.0005
> k2 = 20.0
> 
> ## Initialize x[1]=500 and y[1]=0
> x[1]=500
> y[1] = 0
> 
> for (i in 2:2000){
+   x[i]=x[i-1]+x[i-1]*-k1*dt
+   y[i] = y[i-1]+x[i-1]*k1*dt-y[i-1]*k2*dt
+ }
> 
> finaldata <- data.frame(x,y)
> head(finaldata)
         x        y
1 500.0000 0.000000
2 498.7500 1.250000
3 497.5031 2.484375
4 496.2594 3.703289
5 495.0187 4.906905
6 493.7812 6.095382

I hope this helps.