在 R 中数值求解 Lotka-Volterra ODE

Question

免责声明：Stack Computational Science 上的交叉帖子

目标：我正在尝试使用 sde package 中的 de sde.sim() 函数对 R 中的 Lotka-Volterra ODE 进行数值求解。我想使用sde.sim() 函数，以便最终将此系统转换为 SDE。所以最初，我从一个没有噪声项的简单 ODE 系统（Lotka Volterra 模型）开始。

Lotka-Volterra ODE 系统：

x = 10 和 y = 10 的初始值。

alpha、beta、delta 和 gamma 的参数值分别为 1.1、0.4、0.1 和 0.4（模仿这个example）。

尝试解决问题：

library(sde)
d <- expression((1.1 * x[0] - 0.4 * x[0] * x[1]), (0.1 * x[0] * x[1] - 0.4 * x[1]))
s <- expression(0, 0)
X <- sde.sim(X0=c(10,10), T = 10, drift=d, sigma=s) 
plot(X)

但是，这似乎不会产生捕食者和猎物种群的良好循环行为。

预期输出

我使用 R 中的 deSolve 包来生成预期的输出。

library(deSolve)
alpha <-1.1
beta <- 0.4
gamma <- 0.1
delta <- 0.4

yini <- c(X = 10, Y = 10) 
Lot_Vol <- function (t, y, parms) {
  with(as.list(y), {
    dX <- alpha * X - beta * X * Y 
    dY <- 0.1 * X * Y - 0.4 * Y
    list(c(dX, dY))
  }) }
times <- seq(from = 0, to = 100, by = 0.01)
out   <- ode(y = yini, times = times, func = Lot_Vol, parms = NULL)
plot(y=out[, "X"], x = out[, "time"], type = 'l', col = "blue", xlab = "Time", ylab = "Animals (#)")
lines(y=out[, "Y"], x = out[, "time"], type = 'l', col = "red")

问题

我认为漂移功能可能有问题，但是，我不确定是什么。在尝试解决 sde.sim() 中的这个 ODE 系统时出了什么问题？

Answer 1

假设未指定方法采用列表中的第一个方法，并且所有其他未指定的参数采用默认值，则您正在执行步长 h=0.1 的 Euler 方法。

众所周知，对于具有凸同心轨迹的函数，欧拉方法将产生向外的螺旋。作为一阶方法，错误应该增长到大约T*h=10*0.1=1。或者，如果要进行更悲观的估计，则误差的大小为(exp(LT)-1)*h/L，而L=3 在一些适应的规范中给出了3.5e11 的规模。

---

探索欧拉法的实际误差e(t)=c(t)*h，得到如下图。左边是分量的误差，右边是欧拉方法中各种步长的轨迹。左图中函数c(t) 的误差系数按系数(exp(L*t)-1)/L 缩小，以获得较大时间间隔内的可比值，L=0.06 值达到最佳平衡。

可以看出实际的错误

abs(e(t))<30*h*(exp(L*t)-1)/L

介于线性和指数误差模型之间，但更接近线性模型。

---

要减少错误，您必须减小步长。在SDE.sim 的调用中，这是通过设置参数N=5000 或更大以获取步长h=10/5000=0.002 来实现的，这样您就可以希望前两位数字正确，错误范围为30*h*T=0.6。在 SDE 的情况下，您在每一步中都会累积大小为 sqrt(h) 的高斯噪声，因此 O(h^2) 的截断误差是随机数的一个相当小的扰动。