如何求解和绘制 R 中的微分方程?
how to solve and plot a differential equation in R?
我想求解并绘制指数增长的微分方程,但我不太了解如何使用 deSolve 库。我的等式是 N = N_0 * e^(rt) 我试过的代码是
library(deSolve)
## Time
t <- seq(0, 5, 1)
## Initial population
N0 <- 2
## Parameter values
r = 1
fn <- function(t, N0, r) with(r, list(N0 * exp(r*t)))
## Solving and ploting
out <- ode(N0, t, fn, params)
plot(out, lwd=2, main="exp")
但我希望的输出不是我想要的。我要获取的图表如下:
希望你能帮助我。谢谢
模型函数fn
应该包含导数,然后积分由求解器完成。一阶增长当然可以解析求解,但对于更复杂的模型,这并不总是可行的。
library(deSolve)
## == derivative ==
fn <- function(t, N, r) {
# dN/dt = r * N
list(r * N)
}
r <- 1 # Parameter value
N <- 0:100 # sequence of N
t <- 0 # dummy as the derivative is not time dependent
plot(N, fn(t, N, r)[[1]], type="l")
## == integration ==
t <- seq(0, 5, .1) # time
N0 <- 2 # initial state
## numerical solver
out <- ode(N0, t, fn, r)
plot(out, lwd=2, main="exp")
## for comparison: analytical integration
lines(t, N0*exp(r*t), lwd=2, lty="dotted", col="red")
或者您可以尝试 curve
函数。
op <- par(mfrow=c(1, 2), mar=c(5, 5, 4, 3))
curve(r*x, from=0, to=100, xlab="N", ylab=bquote(dot(N)), main=bquote(dot(N)==N))
curve(N0 * exp(r*x), from=0, to=5, xlab="Time t", ylab="N(t)", main="Exponential growth")
par(op)
我想求解并绘制指数增长的微分方程,但我不太了解如何使用 deSolve 库。我的等式是 N = N_0 * e^(rt) 我试过的代码是
library(deSolve)
## Time
t <- seq(0, 5, 1)
## Initial population
N0 <- 2
## Parameter values
r = 1
fn <- function(t, N0, r) with(r, list(N0 * exp(r*t)))
## Solving and ploting
out <- ode(N0, t, fn, params)
plot(out, lwd=2, main="exp")
但我希望的输出不是我想要的。我要获取的图表如下:
希望你能帮助我。谢谢
模型函数fn
应该包含导数,然后积分由求解器完成。一阶增长当然可以解析求解,但对于更复杂的模型,这并不总是可行的。
library(deSolve)
## == derivative ==
fn <- function(t, N, r) {
# dN/dt = r * N
list(r * N)
}
r <- 1 # Parameter value
N <- 0:100 # sequence of N
t <- 0 # dummy as the derivative is not time dependent
plot(N, fn(t, N, r)[[1]], type="l")
## == integration ==
t <- seq(0, 5, .1) # time
N0 <- 2 # initial state
## numerical solver
out <- ode(N0, t, fn, r)
plot(out, lwd=2, main="exp")
## for comparison: analytical integration
lines(t, N0*exp(r*t), lwd=2, lty="dotted", col="red")
或者您可以尝试 curve
函数。
op <- par(mfrow=c(1, 2), mar=c(5, 5, 4, 3))
curve(r*x, from=0, to=100, xlab="N", ylab=bquote(dot(N)), main=bquote(dot(N)==N))
curve(N0 * exp(r*x), from=0, to=5, xlab="Time t", ylab="N(t)", main="Exponential growth")
par(op)