Modeling with ODEs in Matlab – Part 5A

We are going to wrap up this tutorial series with a fun exploration of complex systems. Complex systems behave in unpredictable ways. This often makes it difficult to design and use models to examine their behavior. In this lesson we will look at some hallmarks of complex systems and examine a canonical example. Finally, in the next installment we will look at how differential equation models of complex systems can be difficult to examine using numerical solutions.

If you have not read the previous lessons, I encourage you to start there (at least read through Part 2). Lesson 1 discussed the meaning of an Ordinary Differential Equation and looked at some simple methods for solving these equations. Lesson 2 looked at the Runge-Kutta approach to solving ODEs and showed us how to use Matlab’s built in function to do so. Lesson 3 looked at two sample multi-population ODE models and explored how to interpret ODE models of living systems. Lesson 4 examined non-linear optimization approaches to fitting ODE modeled systems to empirical data. Lesson 5 continues soon!

A complex system can be made up of many different components, but in my experience there are a few key traits that most complex systems all share.

  1. First and foremost, all complex systems must be non-linear. Linear systems are easily solved using simple analytic techniques. As we have seen from previous lessons, non-linear systems cannot be solved analytically and must be approximated numerically.
  2. Due to their non-linearity, most complex systems are made up of simple independent components. These components react with each other in unpredictable ways.
  3. Most complex systems exhibit broken symmetry on multiple scales. Broken symmetry refers to systems that have critical thresholds past which the system’s behavior transitions sharply.
  4. Although not required, often complexity may arise due to feedback loops among components. You can imagine the scenario where a microphone is placed in front of a speaker and a simple phrase is spoken. This feedback loop creates a very complex and unpredictable reaction.

Complex systems created using these traits will exhibit certain key behaviors.

  1. The most important property of complex systems is their sensitivity to initial conditions. This means that if you change the starting position of the system by a very small amount, the change in the output will eventually become intractably large. This is the key reason why we can’t predict the outcomes of complex systems beyond a certain point.
  2. A complex system will often exhibit emergent behavior that has order seemingly beyond the scope of its components.
  3. Unlike purely chaotic systems, complex systems will often be bounded. This means that while we cannot predict their behavior precisely, we will be able to make statements regarding the scope of the system. This will be examined later when we look at the canonical Lorenz strange attractor.

Before we look at some examples of complex systems, let’s take a quick look at a system that is not complex. To do so, let’s go back to our first example, the exponential growth equation:

Exponential Growth

Let’s see how it lines up with our two lists above. First, this is a linear system. We know this because the only term is the x(t). It is not multiplied by another variable and it is not taken to a power that is not one. That alone is enough to classify it as non-complex, but let’s keep going. Second, this equation is made up of a simple component, but because there is only one component there is nothing for it to interact with in non-linear ways. Third, there is no broken symmetry or sharp thresholds in this equation. Finally, because this is an autonomous system (no t parameter in the equation other than the default x(t)), there can be no feedback.

I will use the exponential.m file from Part 2:

% exponential.m
% Imlements the function dx/dt = x
% Inputs:
%   t - Time variable: not used here because our equation
%       is independent of time, or 'autonomous'.
%   x - Independent variable: this is our population
% Output:
%   dx - First derivative: the rate of change of the
%        population, given the population
function dx = exponential(t, x)
  % The exponential ODE function is simple!
  dx = x;

I will run ode45 twice, once with an initial value of 1.0, and once with an initial value of 2.0.

>> [t, y] = ode45('exponential', [0 10], 1.0);
>> [t2, y2] = ode45('exponential', [0 10], 2.0);
>> plot(t, y, t2, y2)
>> legend('X(0) = 1', 'X(0) = 2', 2)

and here is the result:


So, does this system exhibit any complex behavior? No! First, the system is not sensitive to initial conditions! What do I mean by this? Let’s look at the case where we start with x(0) = 1 and evaluate x(10). In this case, x(10) will be approximately 22,026. Ok, so what? Well, what happens if we start with x(0) being double what we originally set it to: x(0) = 2. In this case, x(10) equals approximately 44,052… exactly double that of the original case! In fact, for any initial value, the output will simply be scaled by the same factor! This shows a linear dependence on our initial condition (what a coincidence, considering our system is linear). Thus, this system is not sensitive to initial conditions. Yes, changes in the initial conditions do change the output, but they do so predictably and the changes are bounded by a known factor! Finally, this system does not have any emergent behavior beyond its basic structure, nor is the system bounded (it grows exponentially).

So, what does a complex system look like? Let’s take the most simple example of a complex system. This equation is known as the logistic map. It has many fascinating details that I will not go into here (I do have suggested follow-up readings at the bottom of this post). The equation is discrete in time and can be described with the following recurrence relation:


This simple equation exhibits some very specific traits! First, it is non-linear due to the multiplication of the X(n) and 1-X(n) terms (this creates a X(n)2 term). This equation is made up of simple components. The equation does not have a clear example of broken symmetry, but it does exhibit feedback as it is a recurrence relation (the value of X(n) depends on the value of X(n-1)).

Ok, let’s take a look! We want to create a function in Matlab that takes in an initial value for X(0) and outputs all the values for X(0) to X(n). Thus, the function will take in three arguments: the initial value, the value for r and the number of values to calculate. It will output a 1D array containing the values of X from 0 to n.

% logistic.m
%   Returns the array of values for the logistic function:
%                X(n+1) = r*X(n)[1-X(n)]
%     given an initial value and a value for r
% Inputs:
%   init - value of X(0)
%   r    - A value for the constant r
%   len  - Number of extra values to calculate
% Output:
%   x    - The values for X from 0 to len
function x = logistic(init, r, len)

% I like to initialize my arrays
x = zeros(1, len+1); 
x(1) = init;

% Loop through the equation one value at a time
for i=2:len+1
    x(i) = r*x(i-1)*(1-x(i-1));

Let’s take a look. For this lesson I am going to ‘randomly’ choose a value for r of 3.75. In case you’re interested, the logistic equation exhibits complex behavior for values of r between 3.57 and 4. I will choose an initial value of 0.5 and calculate the first 100 values past the starting point:

>> x = logistic(.5, 3.75, 100);

I will plot the result as follows:

>> figure('Position',[100 100 960 360])
>> plot([0:100], x)

and viola! (Click for a larger image.)


Immediately we see a few interesting things. The behavior of this system is chaotic. We see where it ended at step 100. Could you predict where it will be at step 110 with much certainty? Likely not! But, given its complex behavior, it is bounded between 0 and 1. What about its sensitivity to initial conditions? This is where it gets fun!

What if instead of starting at 0.5, we start at 0.51? Let’s compare the two:

>> x2 = logistic(.51, 3.75, 100);
>> figure('Position',[100 100 960 360])
>> plot([0:100], x, 'b', [0:100], x2, 'r')


If you look closely, the two systems do mirror each other for about the first 10 steps. Around step 13 or 14, the two start to diverge, after which they never reunite! Thus, a small difference in where we started creates very different behavior, although the behavior is still bounded within the same region.

Well, how much do we gain if we get even closer to where we started? Let’s try a few more starting points and see how they compare!

>> x3 = logistic(.501, 3.75, 100);
>> x4 = logistic(.5001, 3.75, 100);
>> x5 = logistic(.50001, 3.75, 100);
>> x6 = logistic(.500001, 3.75, 100);
>> x7 = logistic(.5000001, 3.75, 100);
>> figure('Position',[100 100 960 960])
>> subplot(6, 1, 1)
>> plot([0:100], x, 'b', [0:100], x2, 'r')
>> title('0.51');
>> subplot(6, 1, 2)
>> plot([0:100], x, 'b', [0:100], x3, 'g')
>> title('0.501');
>> subplot(6, 1, 3)
>> plot([0:100], x, 'b', [0:100], x4, 'm')
>> title('0.5001');
>> subplot(6, 1, 4)
>> plot([0:100], x, 'b', [0:100], x5, 'c')
>> title('0.50001');
>> subplot(6, 1, 5)
>> plot([0:100], x, 'b', [0:100], x6, 'k')
>> title('0.500001');
>> subplot(6, 1, 6)
>> plot([0:100], x, 'b', [0:100], x7, 'r')
>> title('0.5000001');


Each graph starts 10 times closer to the actual (blue) system as the previous one. Do we see a 10x increase in accuracy each time? No! If you look closely, each time we increase our accuracy by a factor of 10, we only delay the divergence of the two models by an extra 10 or so time steps. Thus, we are getting a linear gain for an exponential increase in accuracy. Thus, improving our initial accuracy by a factor of 10 billion will only give us a 10x improvement in our predictive capabilities! A very very small change in our initial condition will result in completely different behavior past a certain point. This is what I mean when I state that complex systems have a sensitive dependence on their initial conditions.

If you ever wanted to know why it’s so hard to predict the weather beyond a few days, this is a great example of the difficulties of modeling complex systems. It takes almost infinite accuracy to be able to predict behavior. Furthermore, Huge gains in accuracy only result in minor gains in predictive capabilities! This is Chaos Theory at its finest!

We will continue and finish this series in Part 5B where I examine a canonical complex system of ordinary differential equations. See you soon!

PS – If you are interested in learning more about complex systems, I highly recommend these two books:

Complexity: A Guided Tour
The Computational Beauty of Nature

3 thoughts on “Modeling with ODEs in Matlab – Part 5A

  1. I merely such as helpful information an individual give in the content articles. I most certainly will take note of your current weblog and check out yet again in this article consistently. I’m just marginally confident I’m going to find out quite a few completely new stuff correct the following! All the best for one more!

Leave a Reply

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