Modeling with ODEs in Matlab – Part 1

Hi everyone!  Today I am posting the first of a planned five part series on using Matlab to simulate systems of ordinary differential equations (ODEs).  This lesson will explore the meaning of a differential equation and look at a few possible ways to solve it.  Lesson Two will look at better ways to evaluate ODEs. Lesson Three will discuss designing and simulating models using systems of ODEs.  Lesson Four will explore ways to fit these models to empirical data.  Finally, we will examine some fun nonlinear ODEs and discuss ways to deal with their complexity in Lesson Five.

First, what is an ordinary differential equation and why do we care?  A differential equation describes a mathematical function that depends on the derivatives of the variables of the equation.  This is a fairly broad and confusing definition, so we’ll simplify things by stating that we are only going to concern ourselves with first-order differential equations of the time domain.  That might not make much sense either, but you can trust me that it simplifies things considerably.  Specifically, because our focus will be on designing and evaluating population models, we will look at systems where the rate of change in population is governed by the size of the population itself.

A simple example of such a formula involves exponential growth.  Here, we assume that the rate of growth is equal to the number of people in the population.  More people, more babies!  The equation can be described by the formula below:

You can see that the equation states that the change of the population over a given time, dx(t), vs the change in time, dt, is equal to the total population at that time: x(t).  The solution to this equation is well known as the exponential function: x(t) = e^t given x(0) = 1 (assuming we start with one person).

Unfortunately, once we delve into nonlinear models, we will not be able to find a closed form solution such as this one.  In anticipation, let’s look at ways to evaluate this simple equation without knowing the solution.  This is known as solving the equation numerically.

The naive approach is to use our initial condition to calculate the rate of change, take a small step forward in time, and use that rate of change to figure out a new population value.  This is known as a first order, or Euler, solution.  We can perform this evaluation with the following Matlab code:

x = zeros(11,1); % Initialize the array, we will evaluate to t=10
x(1) = 1; % The initial population is 1

% For each step, the new population is equal to the 
% old population, plus the rate of change from before 
% (which is equal to the old population)
for i = 2:11 
  x(i) = x(i-1) + x(i-1); 

% Let's see the results!

Because Matlab uses 1 based indexing, our time points 0-10 are stored in indices 1-11. Now, let’s generate the real solution and compare the two results:

% Initialize the new array 
sol = zeros(11,1); 

% Evaluate the solution at each step. 
% Remember that time 0 is stored in index 1. 
% exp is Matlab’s implementation of the exponential function. 
for i = 1:11 
  sol(i) = exp(i-1); 

% Compare the two solutions 
plot([0:10], x, [0:10], sol) legend(‘x’, ‘sol’, 2)

Here is our result:

Whoops, our numerical solution is way too small! Why is that? By using discrete steps, we are ignoring the fact that the differential equation is continuous. At our initial time of t=0, the rate of growth is equal to x(0), or 1. In our first order solution, we use a growth rate of 1 until we reach t=1. That ignores the fact that as our population is steadily increasing, our growth rate is too. As soon as our population grows a little bit, so should the growth rate. Instead, the rate stays the same until we get to our next point of evaluation. So, what can we do? Let’s try updating our rate of growth more frequently (and let’s update our solution for comparison):

% Initialize the array, we will evaluate to t=10 
% This time we will use a step size of .1 
x2 = zeros(101,1); 

% The initial population is 1 
x2(1) = 1; 

% For each step, the new population is equal to the 
% old population, plus the rate of change from before 
% (which is equal to the old population) 
for i = 2:101 
  x2(i) = x2(i-1) + .1*x2(i-1); 

% Initialize the solution array for finer resolution sol = zeros(101,1); 
% Evaluate the solution at each step. 
% We use a finer-grained time step this time to match our 
% new numerical step 
for i = 1:101 
  sol(i) = exp((i-1)/10); 

% Compare the solutions 
plot([0:10], x, [0:.1:10], x2, [0:.1:10], sol) legend(‘x’, ‘x2′, ‘sol’, 2)

We did better than before, but still not well enough. If you haven’t figured it out by now, we can asymptotically approach the solution by taking smaller and smaller step sizes, but we will never be exactly correct. This is actually true of all numerical simulations, but especially so for first order solutions. While it may seem like this small amount of error is acceptable (and it very well may be for certain simulations), I can assure you that this approach will not work once we move on to complex nonlinear systems. Small perturbations in the initial evaluations will propagate into huge errors down the road. This is not acceptable if we want to start building and evaluating models of complex behavior.

So, how do we do better? The answer lies in the Runge-Kutta method of numerical integration and the use of Matlab’s ode45 function, which will be the focus of the next lesson.

11 thoughts on “Modeling with ODEs in Matlab – Part 1

  1. Nice work explaining a difficult subject. However, how could you possibly miss the obvious reality that the growth rate of a population of one is exactly zero? Unless you were making allowances for asexual reproduction (in which case, you ought to have declared that at the outset), I think that most people would find this discussion confusing.

    Just my opinion. However, I did enjoy reading through it.

    • Hi Jim,

      that’s a great point. To keep things simple, we’re ignoring sexual reproduction. If you look even closer, the model has a larger growth rate for a population of 1.1 than it does for a population of 1. This fractional population artifact is an unfortunate side effect of using continuous functions to model a supposedly discrete population. Note that the growth rate would be positive even if the population was 0.001. There is no easy way to fix this problem other than switching to agent based simulations. In the end, we’re usually looking at populations so large that the issue of fractional people isn’t significant.

Leave a Reply

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