how to solve and plot a differential equation in R?

user_11014533 0
Viewed 492

I want to solve and graph a differential equation for exponential growth but I can't quite understand how to use the deSolve library. My equation is N = N_0 * e^(rt) and the code that I tried is

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")

but the output that I hope is not what I want. The graphs that I want to obtain are the following:

enter image description here

I hope you can help me. Thank you

differential-equations plot r
2 Answers

The model function fn should contain the derivative, the integration is then done by the solver. First order growth can of course be solved analytically, but this is not always possible for more complex models.

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")
user_36775760

Alternatively you could try the curve function.

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)

enter image description here

user_65740380
Related