Modeling with ODEs in Matlab – Part 3

Well, I feel like I should apologize for such a long delay between posts. It’s been a crazy summer that has included some vacation time plus an overseas trip to a conference. Regardless, I’m finally back in the swing of things and ready to write up Part 3! To recap: Lesson 1 and Lesson 2 looked at how ODEs are solved numerically and how higher order solutions are more accurate than naive implementations. Today we’ll look at two simulations of living systems (Lotka-Volterra and SIR). Finally, the series will conclude with a post on model fitting and a post about chaotic systems.

Previous posts explained how numerical solutions work and how Matlab will perform the calculations for you automatically. Today, we’ll put that knowledge to good use. We’ll start with a simple Lotka-Volterra predator/prey two-body simulation. There are two populations in question: the predators and the prey. Prey multiply exponentially, similar to our exponential example in the previous lessons. The difference is that prey are also killed off by the predators at a rate directly proportional to both the predator and prey population. Predators are dependent on prey for sustenance and thus grow at a rate dependent on both the predator and prey population. Predators also die off over time due to age.

This gives us two populations each with two separate terms to add to our equation. Because both predator and prey start with ‘p’, we will use the variable x to refer to the prey population and the variable y to refer to the number of predators. Given that both x and y are functions of time, we have the following system of ODEs:

Let’s break things down. The top equation defines how the population of prey, x, changes in relation to the predator and prey populations. The first term, αx, is positive and therefore defines growth. The growth rate is proportional to x, so we know it is positive exponential growth (as seen in previous lessons). The α term is a constant that controls how fast the exponential growth occurs. The second term, -βxy, is negative and thus describes death. The rate is proportional to both x and y. This is because more prey makes it easier for predators to find prey and also because an increase in predators makes it harder for the prey to hide. In contrast, we assume an increase in predators does not affect the birth rate of the prey and thus there is no y in the growth term. Thus, βxy defines the rate at which predators find and kill prey, at an efficiency controlled by the constant β.

The second equation defines the change in the predator population. The first term, δxy is positive and defines growth. Predators rely on prey for sustenance and thus grow proportionally to both their own population and the prey population at a rate controlled by δ. For the second term, because prey can’t kill predators, the predator population only decays proportionally to their own numbers at a rate determined by γ.

It is important to note that both equations have an xy term. This means that both equations are non-linear and thus very difficult to deal with analytically. This is one of the reasons we will turn to numerical integration.

Now, let’s look at the expected dynamics before we actually plug in the equation. We have two populations, so let’s look at the different combinations of high and low populations of each. We want to pick coefficients such that the model is stable. If both populations are low, there needs to be positive growth to restore balance. This growth should come from the prey population being able to multiply without fear of many predators. Thus, the α term must out-pace the β term for small values of x and y. If both values are high, there needs to be death to restore equilibrium. Thus the prey needs to die off due to a large predator population and the β term will dominate. If there is a lot of prey with few predators, the predator population should grow quickly. Conversely, if there are many predators and few prey, the predators population should diminish. If all of this happens, our system should be able to stay in some sort of periodic equilibrium. To this end, I will somewhat arbitrarily chose the constants to be: α = 1, β = .05, δ = .02, and γ = .5.

First, we will code up the equations as a matlab function. The difference between this code and the code from Lesson 2 is that now we are dealing with two variables. Thus, our x parameter is now a vector containing both the prey and predator populations.

% lotka_volterra.m
%
% Imlements the Lotka-Volterra function
%   dx/dt = alpha x - beta xy
%   dy/dt = delta xy - gamma y
%
% Inputs:
%   t - Time variable: not used here because our equation
%       is independent of time, or 'autonomous'.
%   x - Independent variable: this contains both populations (x and y)
% Output:
%   dx - First derivative: the rate of change of the populations
function dx = lotka_volterra(t, x)
  dx = [0; 0];
  alpha = 1; 
  beta = .05; 
  delta = .02;
  gamma = .5;

  dx(1) = alpha * x(1) - beta * x(1) * x(2);
  dx(2) = delta * x(1) * x(2) - gamma * x(2);

Unless you plan on modifying your coefficients through external code it’s better to define them inside the function. Also note that x(1) refers to the prey and x(2) refers to predators.

Next, we need to use odeset to set our options. Note that now both populations 1 and 2 are non-negative.

% Set our preferences for ode45
% The default relative tolerance is 1e-3.
% To set our output to non-negative, we provide an array listing
%   each population that we want to constrain.  Since this example
%   has two populations, we pass the array [1 2]
options = odeset('RelTol', 1e-4, 'NonNegative', [1 2]);

Finally, we use the ode45 function and plot out the results. We will initialize each population to 10 and will integrate from time 0 to 20. The time units here are completely arbitrary since we are using unitless coefficients.

% Use ode45 to solve our ODE
% Place the time points in a vector 't'
% Place the solution in a vector 'x'
[t,x] = ode45('lotka_volterra', [0 20], [10 10], options);
plot(t,x);
legend('prey', 'predators');

Which gives us a final result pictured here:

Notice the nice cyclical behavior. At first, prey multiply rapidly and predators decline a bit. As the prey population soars, predators find more food and begin to grow as well. Eventually the predators eat too many prey and lose numbers to starvation. This lets the prey rebound and the cycle begins again.

Let’s look at one more example. An SIR model is an epidemiological example of an infection invading a population. In the beginning most people are healthy and the infection spreads slowly. As more people get sick, the infection begins to grow noticeably. Eventually, the population begins to recover and gains an immunity to the infection. Thus, the susceptible population disappears and is replaced by a resistant population. The three populations modeled in the standard SIR equation are susceptible, infected, and recovered.

The first thing to note is that there are only two terms. The first, βSI describes the rate of infection. This term scales positively with both the susceptible population and the infected population. The more people susceptible, the easier it is to find a new host. The more infected people, the more there are to spread the disease. Also note that the exact same quantity that is removed from the susceptible population is added to the infected population. This represents a direct state transition from one population to the other. Similarly, as infected people gradually recover at a rate δI, they are removed from the infected population and added to the recovered population. Finally, note that if you add all three equations together you get a net effect of (dS+dI+dR)/dt = 0. This shows that the total population is stable: there is no addition or removal of people, only changes in their state.

Let’s code this up! First, the function definition (I will set the constants to be β = .003 and δ = 1:

% sir.m
%
% Imlements a SIR infection model
%   dS/dt = -beta SI
%   dI/dt = beta SI - delta I
%   dR/dt = delta I
%
% Inputs:
%   t - Time variable: not used here because our equation
%       is independent of time, or 'autonomous'.
%   x - Independent variable: this contains three
%       populations (S, I, and R)
% Output:
%   dx - First derivative: the rate of change of the populations
function dx = sir(t, x)
  dx = [0; 0; 0];
  beta = .003; 
  delta = 1;

  dx(1) = -beta * x(1) * x(2);
  dx(2) = beta * x(1) * x(2) - delta * x(2);
  dx(3) = delta * x(2);

We will use the same options, but add one extra non-negative population (as we now have three).


options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3]);

Finally, we will model an infection of a population of 1000 healthy individuals and 1 sick person. The time interval will go from 0 to a unitless 10.

% Use ode45 to solve our ODE
% Place the time points in a vector 't'
% Place the solution in a vector 'x'
[t,x] = ode45('sir', [0 10], [1000 1 0], options);
plot(t,x);
legend('S', 'I', 'R');

And here is our result:

Notice how the infection initially picks up steam and spreads quickly through the population but soon runs out of gas as the susceptible population quickly drops. The simulation will continue indefinitely as the infected population will slowly decay and the recovered population will asymptotically approach 1001.

I will end with one final point about running continuous simulations of discrete populations. ODE models are a very nice way to simulate complicated population dynamics. Unfortunately, they can only be at best a good approximation. Populations are not continuous. Through ODE modeling, it is possible to represent positive populations of sizes less than one. Our SIR model would work just as well if we started it with an infected population of one tenth a person. Our Lotka-Volterra system would still be cyclical even if our predator population dropped down to one one-hundredth of an individual. Clearly this is not realistic. The lesson here is that you must always take the time to evaluate your model from a rational perspective. Don’t fall into the trap of only looking at the quantitative dynamics. Before you use ode45, take the time to think about what kind of behavior you would expect to see. Think about stable equilibrium points: population values that result in one or more of the equations evaluating to zero, or no change. Think about what should happen at different relative population levels. When you do finally run the model, you’ll be able to look at it with a much better understanding of what is going on behind the scenes.

35 thoughts on “Modeling with ODEs in Matlab – Part 3

  1. hi,
    I get the same error as some of the people commenting here; “Undefined function or variable” while copying the exact example as you presented over for the SIR model. How can this be solved?

  2. Pingback: K2 Seo Tricks And Tips – Best Trick SEOBest Trick SEO

  3. Hi, a very helpful tutorial indeed but I have an issue with my sir model’s input arguments. I copied the code exactly as you have showed it but each time I try to run it a yellow box pop out saying that I require more input arguments while highlighting the variable t in the sir function dx=sir(t,x). Can someone help me asap!!!

  4. With regard to the Lotka-Volterra program, what was the rationale for assigning both initial values of y as 10? That this occurs at one of two times per cycle where the curves cross and the populations are equal is clear from the solution, but how can it be known ahead of time? Without this knowledge, the program will not be useful for other sets of coefficients.

    • Actually almost any pair of numbers will suffice, as long as they represent valid population levels at the same time. But they will reflect different proportions of predator and prey in the population. Or so it seems with a bit of experimenting.

  5. Hi, thanks for the help 🙂
    I have been trying to run the lotka volterra function but it returns the error: Undefined function or variable ‘t’.
    How do i define t outside the function?
    Why is my code not running?

  6. hi ,
    I have an SEIRS model i.e,ds/dt=B-dS-(kSI/1+alpha1I+alpha2I^2)+vR
    dE/dt=(kSI/1+alpha1 I+alpha2 I^2)-(epselon+d)E
    dI/dt-epselonE-(gamma+d)I
    dR/dt=gamma I -(v+d) R
    For B=15,k=0.398,alpha1=0.7,alpha2=0,d=0.04,v=0.0033,epselon=1,gamma=0.143
    I want help about this how can i plot a graph b/w this perameter . Though i can do the rest but little help i needed

  7. Where do the values of x come from in the function if the output is the derivative and the x is found by integrating?

    I’m not sure how the function can solve for dx when x within the function is unknown and partially what you are solving for.

  8. Pingback: Modeling with ODEs in Matlab – Part 5B | Business Intelligence Info

  9. For the Votka-Voltera model, I am trying to find the Fisher Information using matlab coding to find a multiple integral but unfortunately it does not work, could you advice me please?

  10. Thank you a lot sir this has been so helpful. Can you post some more codes especially those with sensitivity analysis

  11. It’s hard to find your blog in google. I found it on 14 spot, you should build quality
    backlinks , it will help you to rank to google top 10.
    I know how to help you, just type in google – k2 seo tips and tricks

  12. It’s hard to find your blog in google. I found it on 14 spot, you should build quality
    backlinks , it will help you to rank to google top 10.
    I know how to help you, just type in google – k2 seo tips and tricks

  13. Dear Sir,

    What can i do if i want to find the answers for x for various values of beta in the same program?

    Can i use a for loop for that?

    Should i use the loop in the first m file or the m file using the ode function?

  14. I have an SEIR model dy(1)=b-(mu+delta*a)*y(1)-lambda*(1-a)*y(3)*y(1)/(1+alpha*y(3))-lambda*(1-delta)*a*y(3)*y(1)/(1+alpha*y(3));
    dy(2)=lambda*(1-a)*y(3)*y(1)/(1+alpha*y(3))+lambda*(1-delta)*a*y(3)*y(1)/(1+alpha*y(3))-(mu+epsilon)*y(2);
    dy(3)=epsilon*y(2)-(alpha+beta+mu)*y(3);
    dy(4)=delta*a*y(1)+beta*y(3)-mu*y(4);
    i want vary to lambda and Epsilon.
    For different values of epsilon=0.1,0.2,0.3,….
    lambda=0.1,0.2,0.3…
    I want help about this how can i implement. Though i can do the rest but little help mis needed

  15. Hi – I have a question about my dynamics problem.
    It is of the form:

    dy1/dt = F(t)-k1*y1+k2*y2

    Where F(t) is a discrete signal not continuous. How can I implement this using the ode-solver?

  16. wow! what a great blog! thank you so much, your explanations are so clear and easy to follow! I’m understanding a whole lot more now than what I was before~

    Do you have any suggestions on how to keep the solutions of an ODE system within set values? For example, I’m modelling and oscillating system measuring concentration, and being in percent I want my results to be within o and 1. Obviously the NonNegative option helps with half, but how would I go about keeping the solution values under 1?
    Thanks again for a great page!

  17. Hi all,

    Ran into a problem. Trying to used ode45 in a for loop. Basically, the following gives me problems while the same version without k coefficients runs fine. Please hint if you can. Thx.
    for k-1:…..
    [t(k),rv(:,k)] = ode45(@Orbit2bodyODE,[t(k):dt:t(k+1)],rv(:,k),options);
    ….
    rv is an array of state 6×1 vectors

    function drv = Orbit2bodyODE(t,rv)

    mu_grav = 398600.4415*(1E3)^3; % geocentric gravitational constant, m^3/s^2

    drv = zeros(6,1);

    drv(1) = rv(4);
    drv(2) = rv(5);
    drv(3) = rv(6);

    drv(4) = -mu_grav*rv(1) / norm(rv(1:3))^3;
    drv(5) = -mu_grav*rv(2) / norm(rv(1:3))^3;
    drv(6) = -mu_grav*rv(3) / norm(rv(1:3))^3;

    drv = [drv(1); drv(2); drv(3); drv(4); drv(5); drv(6)];

    • Hi Andy,

      It’s hard to tell exactly what’s going on, but at the very least you’re going to have an issue assigning the solution of the ode to rv(:,k). The output in this scenario will be a 2D array consisting of values over time for each of your six parameters. Assigning this to rv(:,k) essentially defines rv to be a 3D array. Unfortunately, Matlab is a bit quirky when it comes to 3D arrays. Sorry I can’t help more.

      PS – The last line in your ODE function is redundant.

  18. hi, i have seen your tutorial, but i have a problem with modifying constants through external code. i can’t get the input to be accepted to ode solver, can anyone help with?

    • Hi Mark,

      Take a look at the next section (4a). The trick is to use global variables so that parameters inside utility functions can be set outside their local scope.

  19. hi, i have a problem using ode solver that i have to use if condition for some function inputs but unable to use so can any body provide some help.

Leave a Reply

Your email address will not be published. Required fields are marked *