r - 如何使用 R 求解和绘制 ODE 系统？

Question

我有一个颂歌系统，

我想根据时间 t 绘制 V1 和 V2。我的代码是

library("deSolve")
library("reshape")
library("tidyverse")

parameters <- c(tau = 0.005, tau_r = 0.0025, mui=0, Ve=0.06, Vi=-0.01, s=0.015, mue=10000)

state <- c(X = 0.015, Y = 0)

Odesolver <-function(t, state, parameters) {
 with(as.list(c(state, parameters)),{
 # rate of change
 dX <- -(1/tau + mue - mui)*X + (Y-X)/tau_r + mue*Ve - mui*Vi
 dY <- -Y/tau + (X-Y)/tau_r

 # return the rate of change
 list(c(dX, dY))
 }) # end with(as.list ...
 }

times <- seq(0, 100, by = 0.01)
out <- ode(y = state, times = times, func = Odesolver, parms = parameters)
out.df = as.data.frame(out)
out.m = melt(out.df, id.vars='time')

p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point() +theme_classic()
print(p)

我做对了吗？还有一种方法可以让我绘制1/t改变的值mue吗？这两者都通过第一颂相关联。

Answer 1

ode 系统的转换看起来是合理的，但参数值产生了一些极端的行为。我不知道模型是关于什么的，所以我随意更改了一些参数以获得更合理的输出。我也根据方程组改变了X和。我不明白你的问题是什么意思。您想比较具有不同值的模拟吗？而且，什么是？YV1V2muet

这里是模型的略微修改版本，并提供了如何绘制输出的替代建议：

library("deSolve")
library("reshape")
library("tidyverse")

## parameters of the original poster show extreme behaviour
# parameters <- c(tau = 0.005, tau_r = 0.0025, mui=0, Ve=0.06, Vi=-0.01, s=0.015, mue=10000)

## parameter values arbitrarily changed, 
## so that the output looks somewhat more plausible. 

parameters <- c(tau = 5, tau_r = 50, mui=0, Ve=0.06, Vi=-0.01, s=0.015, mue=0.001)

state <- c(V1 = 0.015, V2 = 0)

## rename to derivative as these are the equations and not the solver 
derivative <-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dV1 <- -(1/tau + mue - mui)*V1 + (V2-V1)/tau_r + mue*Ve - mui*Vi
    dV2 <- -V2/tau + (V1-V2)/tau_r
    
    # return the rate of change
    list(c(dV1, dV2))
  }) # end with(as.list ...
}

times <- seq(0, 100, by = 0.01)

## ode is the ode solver
out <- ode(y = state, times = times, func = derivative, parms = parameters)

## this is deSolve's built-in plot method
plot(out) 

## version of the original poster
out.df   <- as.data.frame(out)
out.m <- melt(out.df, id.vars='time')

p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point() + theme_classic()
print(p)

## alternative with dplyr and tidyr
out %>% 
  as.data.frame() %>%
  pivot_longer(cols = -1) %>%
  ggplot(aes(time, value, color = name)) + geom_point() + theme_classic()