Solving Initial Value Problems in R Using deSolve

How to solve initial value problems (ODEs) using the deSolve package for R.

This post is a repeat of my earlier post on solving initial value problems, but now using the deSolve package (documentation).

library(deSolve)
library(ggplot2)

options(repr.plot.width=12, repr.plot.height=8)

Solving Linear Growth

The first, and simplest, system that we’ll solve is linear growth. This system is given by

For simplicity, we’ll choose $a=1$. Of course, we can easily solve this exactly - the solution is $y(t) = y_0 e^{at}$. Now let’s see how we can solve it numerically.

For ordinary differently equations, the deSolve package provides the ode function. This function takes in the initial state (y=), the times (times=), any parameters to be passed into the system itself (parms=), and a function func representing the system. The function func has signature func(t, state, parameters), where

• t is the current time of the system
• state is the current state of the system (i.e. the value of all variables)
• parameters represents any parameters to the system The function should return a list of time derivatives of the variables.

parameters <- c(a=1) # here the growth rate is 1
state <- c(y=1) # initial state is 1
times <- seq(0, 5, by=0.1)
# Function which takes the current state and returns a list of derivatives
linearGrowth <- function(t, state, parameters) {
    with(as.list(c(state, parameters)),{
        dydt <- a * y
        list(dydt)
    })
}

While most of this is relatively straightforward, there is one wrinkle - the with function. What this does is make the variables and parameters stored in state and parameters respectively available using just their names, rather than having to access them using parameters[['a']] or whatnot. Although confusing at first, the R documentation does a good job of explaining how it works.

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

# convert to a dataframe for graphing
df <- data.frame(time=out[,"time"], y=out[,'y'])
head(df)

ggplot(df, aes(time, y)) + geom_point() + geom_line()

From this, it seems like our numerical solution to this system is pretty good. But how good is it? Since we can solve this system exactly, we can actually compute the maximum deviation that the numerical solution produced:

df['exact'] <- exp(df['time'])
max(abs(df['y'] - df['exact']))

0.000236225065123108

So our numerical approximation is pretty good!

Two Variables - Lotka-Volterra

Now let’s make this a little more complicated by moving to the two-variable Lotka-Volterra system. This system is given by

Where

• $x$ is the population of the prey species
• $y$ is the population of the predator species

And the parameter

• $b$ represents the growth rate of the prey species in the absence of predation
• $p$ has something to do with how often / likely it is that the predator species eats the prey species
• $r$ has something to do with how much food each member of the prey species provides
• $d$ is the death rate of the predator species

Luckily, our approach to solving this hardly has to change at all.

lotkaVolterra <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        dxdt <- x * (b - p * y)
        dydt <- y * (r * x - d)
        list(c(dxdt, dydt))
    })
}
parameters <- c(b=1, p=1, r=1, d=1)
state <- c(x=0.5, y=0.5)
times <- seq(0, 15, by=0.01)

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

# again, converting to a dataframe for graphing
df <- data.frame(time=out[,'time'], x=out[,'x'], y=out[,'y'])
head(df)

In order to see what is going on with our solution, we’ll first plot both the predator and prey populations against the time.

ggplot(df, aes(x=time)) +
    geom_line(aes(y=x, colour="x"), linewidth=2) +
    geom_line(aes(y=y, colour="y"), linewidth=2) +
    guides(colour=guide_legend(title="Animal")) +
    scale_colour_manual(labels=c("Prey", "Predator"), values=c("red", "blue")) +
    ylab("Population") +
    ggtitle("Lotka-Volterra Predatory-Prey Dynamics") +
    theme(plot.title=element_text(hjust=0.5))

From this, is seems like the system demonstrates some periodic behaviour. Just by eyeballing this graph, it seems like the period is roughly seven years. To see this more clearly, let’s now plot the prey population against the predator one for about one cycle.

ggplot(df[df['time'] < 7,], aes(x, y, colour=time)) +
	scale_colour_gradient(low='blue', high='red', aesthetics='colour') +
	geom_path(linewidth=4) +
	labs(title="Lotka-Volterra Phase Space Diagram") +
	theme(plot.title = element_text(hjust=0.5)) +
	xlab("Prey") +
	ylab("Predator")

This make it far more clear that the system is periodic.

3 Variables - Lorenz System

Let’s take this up another dimension by seeing the chaotic behaviour of the Lorenz system. This set of differential equations is given by

with $\sigma=10$, $\beta=\frac{8}{3}$, and $\rho = 28$ (as originally used by Lorenz).

parameters <- c(sigma=10, beta=8/3, rho=28)
state <- c(x=1, y=1, z=1)

lorenz <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        dxdt <- sigma * (y - x)
        dydt <- x * (rho - z)
        dzdt <- x * y - beta * z
        list(c(dxdt, dydt, dzdt))
    })
}

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

# convert to dataframe
df <- data.frame(
    time=out[,'time'],
    x=out[,'x'],
    y=out[,'y'],
    z=out[,'z']
)
head(df)

# 'melt' the dataframe for graphing
melted_df <- cbind(df['time'], stack(df[c('x', 'y', 'z')]))
head(melted_df)

The first way that we’ll examine this is through plots of each individual variable x, y, and z against the time.

ggplot(melted_df, aes(x=time, y=values)) +
	geom_line() +
	facet_wrap(~ind)

It’s a little hard to tell anything from these! To better see the behaviour, we should plot this in 3D. To do so, we’ll use plotly for R to see th results in three dimensions.

library(plotly)
plot_ly(df, x=~x, y=~y, z=~z, type='scatter3d', mode='lines')

Conclusion

There we have it - several examples of solving initial value problem using the deSolve library! The API for this library is quite nice, and as you can see, the set ups for each different system are remarkably similar to each other. The documentation for this library is also very nicely laid out with clear examples and good explanations of not only what they are doing, but why they’re doing it.

Sources / Further Reading

• deSolve Package (CRAN)
• plotly
• R Documentation - with