Y2 Computational

Physics Laboratory

Exercise 4

Modelling physical systemsa predator-prey simulation

Aims and objectives:

  • To gain exposure to using computers to model complex, physical systems
  • To make use of ‘expertise’ acquired in other parts of the degree program to make predictions of the behaviour of systems, based in inspection and manipulation of equations; to understand results of the code in the light of the predictions
  • To be comfortable using functions, and passing variables between them

Duration: 3 weeks


Introduction

This exercise is intended to introduce you to the use of computers as invaluable tools for the study of complex systems characterised by interdependent variables. The system we have chosen to scrutinise is the simplest form of the so-called ‘predator-prey’ class of models. As the name suggests, these models have important application in a biological and ecological context. However, they are also of use to economists, and physicists. Examples of physical systems for which predator-prey-like models find useful application are in the complex dynamics of turbulence, and the competition between particles for water in rain clouds.

The system we will study is the Lotka-Volterra model, which is described by two linked differential equations, and was developed, independently, by two mathematicians. Vito Volterra was an Italian mathematician who was asked by his son-in-law (a biologist) to try and explain the varying population numbers of fish in the Adriatic sea, on the basis of available data on growth, death and interactions rates of predator and prey species. Alfred Lotka developed the same set of equations for describing a hypothetical chemical reaction. The most famous example of the use of this simple model in a real biological system has been for explaining numbers of Canadian lynx (predator) and snowshoe hare (prey) [of which more later].

In this exercise you will first make some predictions of the expected behaviour of the Lotka-Volterra system, from inspection and manipulation of the mathematical form of the equations. You will then write a C++ program to solve, numerically, the equations to show the behaviour over time. You will be making use of numerical techniques very similar to those introduced in Exercise 2.

The Lotka-Volterra Model

To introduce the model we will use the language of biology, and talk in terms of species of predator and prey, death, food etc. The system describes the interaction between one predator species and one prey species. Several, key assumptions are implicit in the model (assumptions which would be unrealistic in most real biological scenarios). First, we assume the prey has an unlimited supply of food. In the absence of any form of predator, the numbers of prey, x, would grow exponentially over time, t, i.e.,

.[1]

In the above, the constant, a, provides suitable normalization and reflects, implicitly, natural birth and death rates of the population. Next, let us introduce the effect of the single predator species, y. We assume the rate at which predators ‘find’ prey is proportional to the product of the total numbers of predators and prey; and that some fraction of these encounters leads to death of prey. In mathematical form, we therefore now have:

.[2]

Here, the constant b fixes the death rate.

This is all well and good, but it says nothing about the predator population. In the absence of any source of prey, the number of predators would fall exponentially (at a rate fixed by some constant, c):

.[3]

Let us assume predators rely solely on the prey species for their food. The growth rate of the predator population is assumed to go in proportion to the product of the two population totals. We might choose to normalize this with the same constant, b, which was used to fix the death rate of prey. However, to keep things as general as possible we use a fourth constant, p, so that:

.[4]

Notice the analogy between Equations [2] and [4], and the importance of the signs of the terms. These two equations give us the Lotka-Volterra model:

[5]

Problem 1: Getting to grips with the modelsome predictions

For the first part of the exercise, we begin by making some predictions of the likely behaviour of the model, based on inspection (and some manipulation) of Equations [5].

The problemwhat is required? Your write-upwhat is needed?

  1. Look carefully at the equations. What are the expected equilibrium values for x and y? Include a clear explanation in your notebook. [1 mark]
  1. Combine the equations to eliminate dt; separate the variables and integrate the equations. You should get the following result, which describes a family of curves:

.[6]

Include a neatly laid out derivation in your notebook. [1 mark]

  1. What form do these curves take? How might the values from question (1) fit in? What do the forms of the curves imply for how x and y might vary over time? [1 mark]

You should use your ‘curve sketching’ skills to help you here. It may also help to test some ideas, and check some numbers, with MathCAD or Excel. (Hint: think about the case where all the constants are unity.) Answer the questions, and include clear, concise explanations of your logic in your notebook, together with a rough sketch of the expected forms for the curves. (Very rough, but neat, will do!)

Problem 2: A numerical solution of the Lotka-Volterra model

Equations [5] can be solved numerically, and estimates of the population numbers, x and y, made as a function of time using techniques like Euler’s method, which we met in Exercise 2. Let us recap briefly.

Here, the small increment h used in Euler’s method is a small increment in time (say t). The Euler formulas for the Lotka-Volterra model are then:

[7]

Provided we have initial values for x and y, and chosen values for the constants a, b, c and p, Equations [7] can be used to estimate x and y versus time. A plot of pairs of x and y can then also be made (the ‘trajectories’ of the system; cf. question (2) above).

Exercise 2 also introduced us to some of the shortcomings of Euler’s method. Accuracy can be lost, even for very simple differential equations. Here, you will use a slightly more sophisticated version of Euler’s method.

In the basic method, the slope of the function is determined only at the starting-point of the step. We know this can cause problems if the slope changes significantly over the interval. To help retain accuracy in the solution a small increment, h, is therefore required. A desirable improvement would be to return an estimate of the slope averaged across the entire range of the increment. But this would need an estimate of the slope at the as yet unknown end-point of the interval!

There is a way forward. We can do things in two steps. First, a predictor stepwe use the basic Euler method to get an estimate of the function (say ŷ) at the end of the increment. We then use ŷ in place of the actual value of y (t+t) to get an estimate of the slope at the end of the incrementthis is the second, corrector step. We can then get our prediction for the next value of y using an average value of the estimated slopes from the start and end points of the increment.

This is one example of the family of Runge-Kutta methods. In fact, it is a second-order Runge-Kutta method, i.e., we use an evaluation of the function at two points along each step to get the slope.

In this exercise we are going to use a more accurate fourth-order approach, where the mean slope is instead determined from four points along the step. Derivations of the fourth-order formulae (which make use of the Taylor expansion) can be found in numerical textbooks. The formulae you will need are as follows:

Where:

To clarify: f (x, y) is the right-hand side of the first Equation in [7], assessed at time t. The more complicated f (x + 0.5k1, y + 0.5m1) means x(t) should be substituted by x(t)+0.5k1 before the right-hand side is computed; and, likewise, y(t) should be substituted by y(t)+0.5m1 prior to calculation. And so on through the progression.

The problemwhat is required?

  1. Write a C++ program to generate estimates of the predator and prey populations, x and y, over time, using the fourth-order Runge-Kutta method.
  1. Code should be written to allow the user to alter the values of the constants, and the step increment t. It should output the results to file [i.e., t, x(t), y(t)] to facilitate graphical output of the results, using MathCAD or Excel.
  1. Start by investigating the case where the constants a = b = c = p = 1. Plot, on the same graph, estimates of x(t) and y(t) versus time, covering t = 0 up to about 30 seconds. [You will need to find a suitable value for t for the computations.] Also produce a plot of y(t) versus x(t) for this time period. Is the phase relation between the curves roughly what you might expect? How do these results compare with your predictions from Problem 1?

Your write-upwhat is needed?

  1. You should include a printout of your code, and a description of it. Detail should be sufficient to allow the marker to easily follow, and understand, the logic and function of the code. The code itself should be well commented. [2 marks]
  1. You should include clearly labelled plots of the population numbers versus time, and versus each other, for the case where a = b = c = p = 1. [0.5 marks]
  1. You should comment on the relation of the two curves: is this what you might naively expect (think in ‘predator’ and ‘prey’ terms). Also discuss the results in the light of your predictions from Problem 1. How do your predictions of the equilibrium numbers fit in? Does your prediction for the shapes of the curves match what is seen? [1.5 marks]

Getting started: Passing variables between functions

To keep things tidy, you ought to be thinking in terms of dividing your code into functions to do the job. Your main program could contain a loop, over the increments in time; inside the loop functions would then be called to calculate the Runge-Kutta estimate for that particular time step. Think carefully about what is required (you may even want to have a function call another function).

If you have forgotten some of the procedures and syntax pertaining to functions, you will need to look back at the appropriate pages of the Y1 C++ course.

‘Introduction to User-defined functions in C++’

There are different types of functions. Those that have a few variables passed to them, and which then use these values to calculate a value, which is passed back by the function:

‘Functions that return values’

Then there are functions which have parameter values passed to them, but which do not return a function value:

‘Functions with parameters and no return value’

You might ask: how might this second set of functions be useful? They might simply do some final calculations in the routine, or output the data etc. However, they can also modify the parameter values passed to them, and return the modified values to the main program, so that other routines can work on the modified values. This sort of function will be of use to you here. The population numbers, x and y are variables that change over time. Initial values might be passed to a Runge-Kutta subroutine, modified to get the new values (at the next time step), and then returned to the main program for use as initial values for the next loop.

You must, however, be careful. Unless you tell the program that the function can alter the actual parameter values, it will only work on copies of them, and the values will not be altered in the rest of the program. You therefore need to specify the parameters to be altered are of the ‘call-by-reference’ type.

‘Call-by-reference parameters’

An example of such usage can be found in:

‘Passing parameters between functions’ [From the module pages]

Problem 3: Lynx and hare

Over the course of the second half of the nineteenth, and the first part of the twentieth centuries, the Hudson Bay Company kept records of numbers of pelts of Canadian lynx (predator) and snowshoe hare (prey). These fur data, some of which are reproduced in the table below, may be regarded as a reasonable proxy for the actual populations (after multiplication by a suitable correction factor).

The problemwhat is required?

  1. On the same plot, make a graph of the numbers of predator and prey pelts as a function of time.
  1. See if you can alter the constants in your computer model to match the pelt records. You should take initial values from the table, i.e., x = 30 and y = 4.

Year / Hares (103) / Lynx (103) / Year / Hares (103) / Lynx (103)
1900 / 30 / 4 / 1911 / 40.3 / 8
1901 / 47.2 / 6.1 / 1912 / 57 / 12.3
1902 / 70.2 / 9.8 / 1913 / 76.6 / 19.5
1903 / 77.4 / 35.2 / 1914 / 52.3 / 45.7
1904 / 36.3 / 59.4 / 1915 / 19.5 / 51.1
1905 / 20.6 / 41.7 / 1916 / 11.2 / 29.7
1906 / 18.1 / 19 / 1917 / 7.6 / 15.8
1907 / 21.4 / 13 / 1918 / 14.6 / 9.7
1908 / 22 / 8.3 / 1919 / 16.2 / 10.1
1909 / 25.4 / 9.1 / 1920 / 24.7 / 8.6
1910 / 27.1 / 7.4 /  /  / 

A couple of hints:

  • An average, over time, of the pelt records will give an estimate of each equilibrium value. Remember how these relate to the constants (see Problem 1).
  • Approximate estimates for a, and c, can be made from the steepest rising and falling parts of the lynx and hare curves. Assume exponential growth (or decay) dominates; simplified equations then describe the evolution of the hare or lynx over time. Manipulate these into a form where you can use data on the ‘steepest’ part to estimate the constants.

Your write-upwhat is needed?

  1. You should include a plot showing the observed time variation of the fur pelts, and your best prediction from your computer model. State the values of the constants, which best matched the records. You should include a brief discussion of how you were led to the quoted values for the constants (cf. ‘hints’ above). [2 marks]