Lesson preamble

Learning objectives

  • Review numerically solving differential equations
  • Calculate and plot solutions with multiple initial conditions
  • Numerically solve two-dimensional models in R
  • Phase portraits of two-dimensional systems
  • Find fixed points for one-dimensional systems analytically (time permitting)

Lesson outline

Total lesson time: 2 hours

  • Review drawing phase portraits and numerically solving equations in one dimension (10 min)
  • Numerical solutions with multiple initial conditions (15 min)
  • Qualitative analysis of two-dimensional models
    • Numerical solutions for two-dimensional models (15 min)
    • Phase portraits of two-dimensional systems by hand (20 min)
    • Using phaseR to create phase portraits in R (40 min)
  • Finding fixed points analytically in one dimension (20 min)


  • install.packages('phaseR')
  • install.packages('deSolve')
  • install.packages('ggplot2') (or tidyverse)
  • install.packages('dplyr') (or tidyverse)
  • install.packages('tidyr') (or tidyverse)

Recap from last week

  • Drawing phase portraits in one dimension:

    • Fixed points: values of \(N\) at which \(\frac{dN}{dt}\), the rate of change of \(N\), is \(0\). To find fixed points, plot \(\frac{dN}{dt}\) vs. \(N\) and find the place(s) where it crosses the \(x\) axis (\(y = 0\)).
    • Stability: if you start at some \(N\) close to the fixed point but not exactly on it, will you go towards (stable) or away (unstable) from the fixed point? The sign of \(\frac{dN}{dt}\) on either side of a fixed point tells you whether \(N\) will increase or decrease in that area. Draw an arrow to the right if \(\frac{dN}{dt}\) is positive, and draw an arrow to the left if \(\frac{dN}{dt}\) is negative.
  • Numerically solving differential equations in R: starting from an initial population size, calculate a sequence of population sizes using the information contained in the differential equation. The result is a trajectory of \(N\) vs. time.

    • Euler’s method: starting from \(N_0\), find the rate of change of \(N\) given by \(\frac{dN}{dt}\) evaluated at \(N_0\). Approximate \(N_{0 + \Delta t}\) using the following expression:

    \[N_{t+\Delta t} \approx N_t + \frac{dN}{dt} \Delta t\] Use a for loop in R to calculate a series of \(N\) at times of interest \(t\).

    • Using R’s ODE-solver ode: define a function that calculates \(\frac{dN}{dt}\) for your model, making sure that it’s in the correct format with arguments t, state, and parameters. Call the function ode and give it the parameters y = state, times = times, func = logistic_fn, parms = parameters

Finding numerical solutions with multiple initial conditions

Let’s make a pretty plot of numerical solutions for the logistic equation from a few different starting points.

##   time       X1       X2       X3       X4 X5       X6       X7       X8
## 1 0.00 10.00000 20.00000 30.00000 40.00000 50 60.00000 70.00000 80.00000
## 2 0.01 10.04006 20.06003 30.05997 40.03994 50 59.94021 69.86063 79.76131
## 3 0.02 10.08024 20.12012 30.11988 40.07976 50 59.88083 69.72250 79.52522
## 4 0.03 10.12054 20.18027 30.17973 40.11946 50 59.82187 69.58560 79.29169
## 5 0.04 10.16096 20.24047 30.23951 40.15904 50 59.76332 69.44992 79.06069
## 6 0.05 10.20150 20.30074 30.29924 40.19850 50 59.70517 69.31544 78.83217

This is exactly what we wanted. We have chosen a set of initial conditions (seq(10, 80, by = 10)) and created a dataframe that has one column of times (time) and a column for each trajectory resulting from each initial condition (X1, X2, …). Now we can use ggplot to plot all the trajectories.

This plot very nicely summarizes the behaviour of the logistic equation: it shows the path that the population size \(N\) will take in time from several initial conditions. Importantly, the carrying capacity \(K\) which we set to be \(50\) is clearly a stable fixed point: all the trajectories go towards \(K\).


Modify the plotting code above so that the legend lists the initial conditions for each line instead of X1, X2, etc.

Qualitative analysis of 2D models

Now we’ll move on to look at how to analyze two-dimensional models. Two-dimensional means that there are now two variables in the system (like \(x\) and \(y\), prey and predators) instead of one (like \(N\)).

Many of the concepts and tools from analyzing one-dimensional models will be used here as well, but we will add some things that are specific to more than one dimension.

To motivate some of these new concepts and the idea of modelling in two dimensions, let’s analyze a model of two species interacting.

\[\frac{dx}{dt} = ax - bxy\] \[\frac{dy}{dt} = cx - dy\]

This is very close to the Lotka-Volterra predator-prey model, but the term for predator growth is \(cx\) instead of \(cxy\). In this model the predators will eat the prey at a constant rate, no matter how many of the predators there are. Make sure you don’t use these exact equations when you do the assignment!

First, we’ll simulate a trajectory for the two species.

It looks like we have oscillations that decay to a fixed point. Notice that the prey population (in red) dips below \(1\) during each oscillation. If this was a real population, the prey would have gone extinct on the first dip below \(1\). This means that the parameters we’ve chosen aren’t very realistic. We could instead choose parameters that would change the fixed points to real-world population sizes.


Modify the plotting code above to include a legend distinguishing ‘predators’ and ‘prey’. Feel free to look up a solution online.

A note on using Euler’s method to solve higher-dimensional systems

Note: if you are using Euler’s method to solve a two-dimensional system, make sure you calculate the rates of change using the variables from the previous time step. For example, in the predator-prey model, both \(\Delta x\) and \(\Delta y\) depend on \(x\) and \(y\), and it’s important that you calculate \(\Delta x\) and \(\Delta y\) before updating the values of \(x\) and \(y\).

Here’s an example.

This is the right way to do it.

This is the wrong way.

Phase portraits in two dimensions


What were the axes of the one-dimensional phase portraits we drew? What would the axes of a phase portrait for the predator-prey system above be?

Continuing with the predator-prey example: instead of plotting our numerical trajectory as \(X\) and \(Y\) vs. \(t\), we’ll plot \(Y\) vs. \(X\): predator population size vs. prey population size.

A phase portrait in two dimensions has two axes: one for each of the variables in the system. Here, the two variables are \(X\) and \(Y\), prey and predators. We want to include the same information in this 2D phase portrait as in the 1D portrait: at minimum, we want to mark the fixed points and draw a vector field: arrows indicating in which direction the rate of change is for different values of \(X\) and \(Y\).

We can also add other information, and in the plot above we haven’t yet drawn arrows or fixed points, but we have drawn a trajectory - the dots trace out what happens to \(X\) and \(Y\) over time starting from \(X = 2\), \(Y = 2\). Time is no longer an axis in the plot, but it’s still there, hidden in the path that the populations take.

Comparing these two ways of plotting a trajectory, we can see that oscillations in time look like circular motion in the phase portrait. The GIF below illustrates this.

Predator-prey simulation

Predator-prey simulation

Drawing qualitative phase portraits by hand

Before we look at how to add more to phase portraits in R, let’s practice the process by hand to get familiar with it. Let’s use the following two-dimensional system as an example. The two state variables are \(x\) (position) and \(v\) (velocity), and this is a model of a mass attached to a spring. \(dx/dt\) is the rate of change of position, and \(dv/dt\) is the rate of change of velocity.

\[\frac{dx}{dt} = v\]

\[\frac{dv}{dt} = -\omega^2 x\]

We want to draw a vector field in two dimensions to show how the system will evolve through time. Our axes are \(x\) and \(v\). It doesn’t matter which one is the vertical axis and which one is horizontal, but let’s put \(x\) on the x-axis and \(v\) on the y-axis. These variables aren’t population sizes, so they can realistically be negative and non-integers.

Next, we choose some pairs of \(x\) and \(v\) and evaluate both \(dx/dt\) and \(dv/dt\), drawing an arrow at \((x,v)\) in the direction of the derivative. Let’s also set \(\omega = 1\) for convenience. Let’s start with \((0,0)\): both derivatives are 0, so there’s no arrow. For \((0,1)\), \(dx/dt = 1\) and \(dv/dt = 0\), so we draw an arrow pointing horizontally in the \(+x\) direction. Similarly for \((1,0)\). For \((1,1)\), we draw an arrow that goes \(1\) unit in the \(x\) direction and \(-1\) unit in the \(y\) direction. As we keep filling in this grid, we start to get a picture of how the system can behave for different starting points of \(x\) and \(v\).

This is analogous to plotting \(dN/dt\) vs. \(N\) and drawing arrows to the right where \(dN/dt\) is positive and arrows to the left when it’s negative, except that now we have to consider that the rate of change is affecting both variables at the same time.


Define a function called mass_spring that uses the format required by ode to calculate and return \(dx/dt\) and \(dv/dt\). We will use this function later to plot a phase portrait in R. Remember that the arguments for your function must be called t, y, and parameters.

Using phaseR to draw phase portraits in R

phaseR is a package specifically for drawing phase portraits. It can automatically calculate and plot things like trajectories and vector fields.

Let’s plot a vector field and some trajectories for the predator-prey model.

## [1] "Note: colour has been reset as required"


Use phaseR and the function flowField to make a phase portrait for the mass-spring system of differential equations. You will need to use the function you wrote in the previous challenge. Is there a fixed point, and if so, where do you think it is?

Finding fixed points analytically (time permitting)

Up until now we’ve been plotting \(\frac{dN}{dt}\) vs. \(N\) to find the fixed points of a system. But we can also solve \(\frac{dN}{dt}\) \(=0\) analytically to find the fixed points. One advantage of finding the fixed points analytically is that you can see how the fixed points depend on the model parameters, whereas before we had to choose parameters in order to plot \(\frac{dN}{dt}\) and find the fixed points qualitatively.

\[\frac{dN}{dt} = rN(1-\frac{N}{K})\]

First, we set \(\frac{dN}{dt}\) \(=0\).

\[0 = rN(1-\frac{N}{K})\]

Next, we find all the values of \(N\) that satisfy this equation. If either \(rN = 0\) or \(1-\frac{N}{K} = 0\), the equation is satisfied. This means we have two fixed points, one at \(N = 0\) (since \(r \neq 0\)) and the other at \(N = K\). These are the same fixed points we found graphically last week. Notice that one of the fixed points depends on \(K\), but that neither of the fixed points depend on \(r\).


Find the fixed point(s) of the following differential equation:

\[\frac{dx}{dt} = x^2 - 1\]

  • First find the fixed points analytically on paper.
  • Plot \(\frac{dx}{dt}\) vs. \(x\) for \(x\) ranging from \(-2\) to \(2\) in \(0.1\) increments.
  • Sketch a phase portrait. Are the fixed point(s) stable or unstable?

This work is licensed under a Creative Commons Attribution 4.0 International License. See the licensing page for more details about copyright information.