R:向量化有限差分方程

Question

我试图把一些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)

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

Answer 1

这不一定是正确的，R是坏的循环。这在很大程度上取决于你在做什么。使用问题中的k1、k2、dt和data (即以k1=5开头的四行)并以迭代矩阵的形式构造问题，下面最后一行中的循环几乎立即返回到我的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, ]

(您也可以尝试将向量存储在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

Answer 2

也许这能帮上忙。

我认为您需要具备data[1,2]的初始条件。在初始条件下，我假设data[1,1]为500，data[1,2为0。

代码是这样的：

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

我希望这能帮到你。