在求解一些诸如动力学方程的时候，可以利用R包 deSolve 来完成

这里举两个实用性的例子：

普通的微分方程组：

library(deSolve)
library(scatterplot3d)
## Chaos in the atmosphere
Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c = 28)
state <- c(X = 1, Y = 1, Z = 1)
times <- seq(0, 100, by = 0.01)

out <- ode(y = state, times = times, func = Lorenz, parms = parameters)

利用函数 ode() 可以计算出三个函数 X，Y，Z 的轨迹函数图像

带有误差延时的微分方程组

首先是单个方程的情况：

其中 (t-1)代表延迟，也就是这一部分中，函数 y() 在 t 时刻 前移了1个单位

library(deSolve)

# 传入方程参数 y 和 t;函数y写在state里面
derivs <- function(t, state, parms) {
  with(as.list(c(state, parms)), {
    if (t < 1){
      dy <- 1
    }
    else{
      dy <- -1*lagvalue(t - 1) + t^2+1
    }
    list(c(dy))
  })
}

state <- c(y = 0)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = derivs, parms = NULL)
plot(yout)

### 或者这样写
library(deSolve)

# 传入方程参数 y 和 t
derivs <- function(t,  y, parms) {
    if (t < 1){
      dy <- 1
    }
    else{
      dy <- -1*lagvalue(t - 1) + t^2+1
    }
    list(c(dy))
}

state <- c(y = 0)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = derivs, parms = NULL)
plot(yout)

其中：

1. dy <- -1*lagvalue(t - 1) + t^2+1代表方程形式为 dy / dx = -y( t - 1 ) + t² + 1
2. 该方程的其中一个解为：当 t < 1 时，y = t；当t ≥ 1时，y = t²

第二种是方程组的形式：

其中(t-1)代表延迟，也就是在这一部分中，函数y() 在 t 时刻前移了1个单位；
而(t-0.5)也代表延迟，也就是在这一部分中，函数x() 在 t 时刻前移了0.5个单位

library(deSolve)

# 注意如果是方程组需要传入待求解函数(本例中为 x,y;因此传入 x,y并存入state)
derivs <- function(t,state, parms) {
  with(as.list(c(state, parms)), {
    if (t < 1){
      dx <- 1
      dy <- 1
    }
    else{
      dx <- -1*lagvalue(t - 0.5,1) + 2*lagvalue(t - 1,2) + y^2 + 1
      dy <- -1*lagvalue(t - 1,2) + lagvalue(t - 0.5,1) + x
    }
    list(c(dy,dx))
  })
}

state <- c(x = 10,y = 0)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = derivs, parms = NULL)
plot(yout)


### 或者这样写
library(deSolve)

# 注意如果是方程组需要传入待求解函数(本例中为 x,y;因此传入 x,y)
derivs <- function(t, x , y, parms) {
    if (t < 1){
      dx <- 1
      dy <- 1
    }
    else{
      dx <- -1*lagvalue(t - 0.5,1) + 2*lagvalue(t - 1,2) + y^2 + 1
      dy <- -1*lagvalue(t - 1,2) + lagvalue(t - 0.5,1) + x
    }
    list(c(dy,dx))
}

state <- c(x = 10,y = 0)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = derivs, parms = NULL)
plot(yout)

其中：

1. lagvalue(t - 0.5,1) 代表函数 x( t - 0.5 )；lagvalue(t - 1,2) 代表函数 y( t - 1 )；1 代表 state里面的第一个函数 x()，2 代表 state里面的第二个函数 y()
2. state 代表初始值函数，对应 x 和 y 的初始值
3. derivs <- function(t, x,y, parms) {} 如果待求解函数 （dx，dy），则需要传入 x 和 y两个参数

以不同时间段设置不同的延时:

微分方程形式：

library(deSolve)

# lamb为可变的延时系数
func_dirct <- function(t, state, parms) {
  with(as.list(c(state, parms)), {
    if (t < 1){
      dy <- 1
    }
    else if (t < 2 & t > 1) {
      lamb = 1
      dy <- -1*lagvalue(t - lamb,1)
    }
    else {
      lamb = 0.1
      dy <- -1*lagvalue(t - lamb,1) + y
    }
    list(c(dy))
  })
}

state <- c(y = 10)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = func_dirct, parms = NULL)
plot(yout)

### 或者这样写
library(deSolve)

# lamb为可变的延时系数
func_dirct <- function(t,  y, parms) {
    if (t < 1){
      dy <- 1
    }
    else if (t < 2 & t > 1) {
      lamb = 1
      dy <- -1*lagvalue(t - lamb,1)
    }
    else {
      lamb = 0.1
      dy <- -1*lagvalue(t - lamb,1) + y
    }
    list(c(dy))
}

state <- c(y = 10)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = func_dirct, parms = NULL)
plot(yout)

其中：
lamb 为可变的延时变量，用if语句进行控制

如果是想让延时系数与时间 t 有关，则需要进行如下改下：

library(deSolve)

# 连续变量型的延时系数书写方式,传入参数 parms <- c(lamb = 1)
func_dirct <- function(t, state,parms) {
  with(as.list(c(state, parms)), {
    # 对延时系数lamb进行条件筛选
    if (t < 1){
      lamb <- 1
    }
    else if (t < 2 & t > 1) {
      lamb <- t/100000
    }
    else {
      lamb <- 3*t/100000
    }
    
    print(lamb)
    if (t < 1){
      dy <- 1
    }
    else if (t < 2 & t > 1) {
      dy <- -1*lagvalue(t - lamb,1)
    }
    else {
      dy <- -1*lagvalue(t - lamb,1) + y
    }
    list(c(dy))
  })
}

parms <- c(lamb = 1)
state <- c(y = 10)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = func_dirct, parms = parms)
plot(yout)

### 或者这样写
library(deSolve)

# 连续变量型的延时系数书写方式, 把参数 lamb 写在 function() 里面
func_dirct <- function(t, y, lamb,parms) {
  # 对延时系数lamb进行条件筛选
  if (t < 1){
    lamb <- 1
  }
  else if (t < 2 & t > 1) {
    lamb <- t/100000
  }
  else {
    lamb <- 3*t/100000
  }
  
  print(lamb)
  if (t < 1){
    dy <- 1
  }
  else if (t < 2 & t > 1) {
    dy <- -1*lagvalue(t - lamb,1)
  }
  else {
    dy <- -1*lagvalue(t - lamb,1) + y
  }
  list(c(dy))
}

parms <- c(lamb = 1)
state <- c(y = 10)
times <- seq(0, 2, 0.1)

yout <- dede(y = state, times = times, func = func_dirct, parms = parms)
plot(yout)

对延时系数 lamb 进行条件判断以及设置 lamb 和 t 的函数关系式，可以进行可变延时系数的计算

对函数值进行判断，对应不同的微分方程：

library(deSolve)

# 利用函数值做判断,原理是迭代计算每一个time对应的函数值
derivs <- function(t,x, y, parms) {
  print(x)
  if (x[1] == 10 & x[2] == 0){
    dx <- 1
    dy <- 3
  }
  else if (x[1] < 7 & x[2] < 1.2) {
    dx <- -1*x + y^2
    dy <- -5*y + x*y^2 + x^2*y
  }
  else if (x[1] < x[2]) {
    dx <- -1*x + y^2
    dy <- -5*y + x*y^2 + x^2*y
  }
  else {
    dx <- 2*x
    dy <- 5*y
  }
  list(c(dy,dx))
}

state <- c(x = 9,y = 0)
times <- seq(0, 0.5, 0.05)

yout <- dede(y = state, times = times, func = derivs, parms = NULL)
plot(yout)

## 这种情况就不能用 with(as.list(c(state, parms)), {}) 语句了

其中

1. 在函数derivs中print(x)打印出来的是
   
   分别对应 x 和 y 的值，因此 x[1] 代表的是 x(t)函数随时间 t 变化的值；x[2] 代表的是 y(t)函数随时间 t 变化的值；因此可以依据 x[1] 和 x[2] 进行判断
2. 这种情况就不能用 with(as.list(c(state, parms)), {}) 语句了