1

我正在尝试将一些 Fortran 代码移至 R，以获得与化学动力学相关的有限差异。

示例 Fortran 循环：

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

类似的 R 代码：

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)

这会覆盖我在数据框中的第一个值，并在其他地方留下零。我想运行这个矢量化而不使用 for 循环，因为它们在 R 中太慢了。此外，我的函数到目前为止只计算第一列。在第一列的语法正确之前，我无法让第二列正常工作。

4

2 回答 2

1

可能这会有所帮助。

我认为你需要有初始条件data[1,2]。我在初始条件下都假设data[1,1]为 500 和0。data[1,2

代码如下：

> ## 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

我希望这有帮助。

于 2013-04-19T02:45:44.543 回答
1

R 不擅长循环并不一定是真的。这在很大程度上取决于你在做什么。使用k1、k2和问题（即以 开头的四行）dt并根据迭代矩阵制定问题，下面最后一行中的循环几乎立即在我的 PC 上返回：datak1=5

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, ]

（您也可以尝试将向量存储在 z 的列中而不是行中，因为 R 按列存储矩阵。）

这是结果的第一位：

> 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
于 2013-04-19T03:19:48.693 回答