R 中关于欧拉方法(微分方程)的警告消息

Warning message on Euler Method (Differential Equations) in R

我正在研究微分方程和欧拉方法,但在我的程序中找不到错误。

想法是使用欧拉方法模拟微分方程的估计。

这是我们在 class 中尝试过的一些代码。

# Define differential equation
y.prime.diff.eqn <- function(p, y) {return(5*y-p)}


initial.condition.x1 <- 0
initial.condition.y1 <- 0

# Define Euler Method's estimations
euler <-  function(x1 = initial.condition.x1,
                   y1 = initial.condition.y1,
                   y.prime = y.prime.diff.eqn(x1, y1),
                   iter = 5,
                   step.size = 1) {

  for (i in 2:iter)

  {

    x1[i]     <-  x1[i-1] + step.size
    y1[i]     <-  y1[i-1] + step.size * (y.prime)
    y.prime   <-  y.prime.diff.eqn(x1[i], y1[i]) 

  }

  return(data.frame(cbind(x1,y1)))

}

output <- euler()

output

它输出正确的结果,但有一条警告消息:

Warning message:
In y1[i] <- y1[i - 1] + step.size * (y.prime) :
number of items to replace is not a multiple of replacement length

为什么我会收到此警告?

你遇到了奇怪的范围问题,因为函数参数中的函数在函数体内被调用....而你认为它是在单个值上调用的,例如x1[i] 但它看起来像是在增长向量 x1 上调用的。只需像这样从迭代函数中取出起始参数:

y.prime.diff.eqn <- function(p, y) {return(5*y-p)}


initial.condition.x1 <- 0
initial.condition.y1 <- 0
y.prime = y.prime.diff.eqn(initial.condition.x1, initial.condition.y1)

# Define Euler Method's estimations
euler <-  function(x1 = initial.condition.x1,
                   y1 = initial.condition.y1,
                   iter = 5,
                   step.size = 1) {

  for (i in 2:iter)

  {

    x1[i]     <-  x1[i-1] + step.size
    y1[i]     <-  y1[i-1] + step.size * (y.prime)
    y.prime   <-  y.prime.diff.eqn(x1[i], y1[i]) 

  }

  return(data.frame(cbind(x1,y1)))

}

output <- euler()

output

同样的答案,没有警告。

我的一位同学想出了一个没有警告的不同解决方案。

我很想知道什么在性能方面更好。据我了解,不在每次迭代中计算下一个 y.prime 将是一种改进。但与此同时,我发现斯蒂芬的回答更简洁。

# Define differential equation
y.prime.diff.eqn <- function(p, y) {return(5*y-p)}


initial.condition.x1 <- 0
initial.condition.y1 <- 0

# Define Euler Method's estimations
euler <-  function(x1 = initial.condition.x1,
                   y1 = initial.condition.y1,
                   y.prime = y.prime.diff.eqn(x1, y1),
                   iter = 5,
                   step.size = 1) {

  for (i in 2:iter)

  {

    x1[i]     <-  x1[i-1] + step.size
    y1[i]     <-  y1[i-1] + step.size * y.prime.diff.eqn(x1[i-1], y1[i-1)

  }

  return(data.frame(cbind(x1,y1)))

}

output <- euler()

output