在求解一些诸如动力学方程的时候，可以利用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)
其中：
dy <- -1*lagvalue(t - 1) + t^2+1代表方程形式为 dy / dx = -y( t - 1 ) + t2 + 1- 该方程的其中一个解为：当 t < 1 时，y = t；当t ≥ 1时，y = t2
第二种是方程组的形式：
其中(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)
其中：
lagvalue(t - 0.5,1)代表函数 x( t - 0.5 )；lagvalue(t - 1,2)代表函数 y( t - 1 )；1 代表 state里面的第一个函数 x()，2 代表 state里面的第二个函数 y()- state 代表初始值函数，对应 x 和 y 的初始值
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)), {}) 语句了
其中
- 在函数derivs中
print(x)打印出来的是
分别对应 x 和 y 的值，因此 x[1] 代表的是 x(t)函数随时间 t 变化的值；x[2] 代表的是 y(t)函数随时间 t 变化的值；因此可以依据 x[1] 和 x[2] 进行判断- 这种情况就不能用
with(as.list(c(state, parms)), {})语句了
网友评论