Welcome to my office.
I'll now give you a look into
how I usually work when I'm interested in
following the behavior of a dynamical system
and I need to rely on
numerical integration or numerical methods to do so.
I will start by introducing
a slightly more general Lotka-Volterra model
than what we've seen so far.
The reason being that the two-species Lotka-Volterra model,
we've been able to study its behavior,
look at its fixed points, flow diagrams
with simple pen and paper, using null clients
and solving for our derivatives equal to zero.
In the case of three-species Lotka-Volterra,
things get a little crazier.
The algebra gets more complex, it gets harder to visualize.
Four-species, even more so.
And so numerical methods are going to be very useful.
We're gonna use a slightly different form.
So this is what's often called a generalized Lotka-Volterra.
It has those logistic-looking terms
like X I times one minus X I,
but it's also defined over a set N of species,
so we would have X one, x two, up to X N for X N species,
and all interactions would be encoded in this A matrix.
So the interaction A I J would represent
the impact of species J on species I.
So for example, in the case of a three-species model,
we would have nine parameters for interaction.
It's important to realize that
because of the negative term in the first equation,
so the fact that the X I times X J term Is a negative term,
it means that by default a positive parameter,
a positive element in the matrix A
would result in a decrease,
so a negative impact on the derivative of X I.
And so positive parameter means something like predation,
So A one two means the species two is hurting species one
if A one two is positive,
and a negative parameter means the opposite.
We're gonna rely on numerical integration.
If I go back, you'll see that in this matrix
throughout this video, I'll keep everything fixed,
but only vary this A three one element
which I'm parameterizing with alpha.
And so you can think of the middle row, species two,
as being a species that predates on the other two
because both species one and two have a positive
because of the negative element here.
It's hard to wrap your arm around that.
Positive impact on species two,
meaning it's sort of a predator.
And by tuning species three,
we're deciding whether it's above.
By tuning alpha, I mean,
we're deciding whether species three
is above or below species one in the trophic level.
That's one way to think about this.
But we don't have to bother too much
about the reality behind this model.
We're really interested in
following its behavior through time,
and the only way to do it well
is to rely on numerical integration.
So from a previous video, you might remember
that the goal of numerical integration
is to capture how X of T is gonna change
and if I make small steps in time,
so if I go to S T plus H forward in time.
And from Taylor series, we know that
value X at time T plus H is gonna depend on the derivative.
Here I'm using primes to denote
the derivative of the solution X of T,
but it's also gonna depend on the curvature X prime prime,
on the third derivative, X prime prime prime,
and so on and so forth.
But all we have is the current value at that time T, X of T,
and a function that gives me my time derivative.
So remember that whether it is,
D X over D T, X that or X prime, that doesn't matter.
That derivative is equal to some function
of X of T and maybe even of time.
So we know how to calculate X prime
because we have this differential equation,
we have that function that gives us X prime.
So Euler's method is going to say
I'm just gonna take what I have,
I'm gonna assume that I'm at point X of T,
I'm following a slope for a length of time H
and I'm gonna know where I am after that time step H,
and I'm gonna ignore the higher order derivatives.
So the local error, as we've said before, is of order two
because we're ignoring this H squared term
which is the most significant
since we're assuming H to be small.
So the local error is gonna be an H square.
The global error over, you know,
a capital T times series that we're interested in
is gonna take capital T over H steps
and therefore capital T over H time H square error,
which gives us the capital T times H global error.
The global error is always gonna be
one order smaller than the local error
since the smaller H is, the more steps we have to do.
Heun's method we've also discussed,
as you'll see in the assignment,
it combines multiple derivatives
to evaluation of the derivative in a smart way
such that it can approximate
or actually capture, I should say,
the second derivative
without having an explicit equation for it.
So that's a really interesting method.
So we're gonna use some of those.
I'll rely on some more advanced integrator
to give us an idea of a ground truth.
There is no, like, ground truth here,
no perfect solution for X of T,
but we can rely on more advanced algorithms
just to compare the error done by those simpler approaches.
So we're gonna look at
this generalized Lotka-Volterra model in greater detail.
To do so, I will take you to my usual workspace.
So you should now be seeing my codes in the top left.
So this is in Octave,
which is an open-source free alternative to Matlab.
If you want to use this code,
I'll be sharing it on the Blackboard site.
There are some special functions,
like calls to this O D E four five integrator,
which you might have to play with a little bit,
but it should almost work out of the box
for your Matlab software.
In the bottom left, what you have is
this command line version of Octave
where I'll be able to simply call my different functions.
So very briefly here, how these codes are organized is
you have a main code, which is "compare integrators"
because that's what I want to do.
I'm gonna be taking this time series, capital T,
which is gonna go from time zero to time max T
which is an argument of my function call,
and I'm gonna do so in time steps of H.
So this, you know, in Python you might use NumPy
and this would be something like
a NumPy arrange function call
to create my different time steps.
And I also have initial conditions X zero,
so that's the starting point for my time series.
And with that I can call O D E four five
which is an advanced Runge Kutta
with adaptive time step integrator
that we're gonna assume is
close enough to the the real trajectory,
and we're gonna compare
what the orders method does against it.
So this line here is creating my time series of T
and Y, the state of my system over time.
And I'm doing the same thing
with my own Euler integrator
using time step H to compare them.
We won't need this for now.
And we'll be comparing,
actually, let me remove all of this for now.
Remove this, keep it as simple as possible.
Okay, so let's look at what these other functions look like.
So O D E four five, I have to pass it my set of equation.
It's always good to define your equation
in a separate function.
The reason for that is if I wanted to use a different model,
I wouldn't have to change
my comparing integrator function here.
I would just have to give O D E four five
a different function to integrate.
So if I go in my predator-prey system,
I've hard-coded my matrix A
because I said the only parameter that I'll be changing
is this alpha, or here lower case A,
and this is in matricial form
calculating the three equations for D X one over D T,
D X two over D T, and D X three over D T
that you saw in the previous slide.
It's doing the exact same form in matricial operations.
The reason for that being that
whenever you use Matlab, remember that the "Mat" in Matlab
is for matrix not mathematics
and therefore it is optimized for matricial operations.
So this makes it a little faster than
if I have three lines for the three different equations.
The other function that I have is for my Euler integrator,
so this Euler pred-prey function,
and if we go in there, we'll see that
I have the same hard-coding of my A matrix
with my interaction parameter A,
and then I wanna go from time T equals zero
to time T F in steps of H.
So again, this is similar to a N P R range.
So this four loop is just doing as many time steps of H
that it needs until it reaches time T F.
At every time step, it's gonna look at the current state.
So this is just asking,
"What is the current state of the system?"
It's gonna calculate
in the same way I was doing in the previous code,
calculate the derivatives at this point,
and then simply add to the current state Y,
add this first term of our Taylor series.
So follow this derivative D X times D T for time H, right?
So D X times D T, or D X D T, I should say,
is simply the speed of our system.
So if you're going somewhere at 10 miles per hour,
and you're doing it for a time step of 10 hours,
you get 10 times 10 miles per hour,
a hundred miles done in 10 hours.
So this is the amount of change done over that time steps,
and we simply keep track of all time
with this little trick here.
Going back to that compare integrators code,
you'll see another good programming practice here, which is
always include a sample call.
So when I first opened this file,
there were some modifications
that were done in discussion groups.
Sometimes you'll modify things and forget about them.
It might take a year until you go back to your code.
It's always useful to be able
and copy a function call to get it running.
So here, if I run this command,
we expect to run our model with alpha equals 0.25.
I have my initial values here set randomly
or set to 0.2, 0.3, 0.4,
but those were not chosen in any meaningful way.
I expect two figures to appear.
Figure one will be a plot comparing the time series
from the O D E four five, our really advanced integrator,
to, in black, the time series obtained
with our Euler integrator.
So red will be the advanced integrator and black the Euler.
And in the second plot, I expect a three D plot
of the trajectory obtained with the Euler integrator.
So here's our first figure.
So again, in red you have the really advanced integrator
and in black, which we obtained with Euler,
so notice that they don't match perfectly, right?
And we know that Euler does better when H is very small.
So if H was five times as big, for example,
the difference would be even bigger.
And if H was instead five times smaller,
it takes a little longer to run
but now my two curves are almost indistinguishable.
So going back to that original one,
and then we get the three D plot of that Euler curve.
So we start in the bottom corner here,
and as our system evolves,
it eventually settles around a steady state,
and we know this by looking at this curve.
Now, I can play with different values of alpha,
so if I make alpha bigger, we'll notice that
an overall time, capital T, equal to 50
is not enough to settle to a steady state value,
and we see the same thing here.
We see an interesting behavior, the start of a swirl here.
So we can integrate for longer, to 200,
to see what that actually looks like.
And we see that we're spiraling toward a fixed point,
eventually reaching it.
And our integrator is still doing really well.
But again, if our time steps were too big,
we might get a little, or a big difference in this case,
in time to convergence to that fixed point.
So one important thing to keep in mind is that
it's really important to run your integrator
with maybe a larger time step at first
just to make sure that things are running.
But if you want to validate a trajectory,
then you really have to decrease your time step
and make sure that you get the same time series.
And once your time series doesn't change much,
with your times that you know
you've converged onto something meaningful.
Let's bring alpha up a little more,
and now we get some really nice oscillations going.
We can integrate further to see
if it's just the beginning of a spiral
like we had with the previous value.
And in this case we see that we're really capturing
sort of what we had in the numerical simulation
of the Lotka-Volterra model.
We have resilient, robust oscillation and a periodic system
where we're now orbiting, our system is orbiting
in this three D space, three different species.
I'm only plotting species one on this top figure right now.
I can increase my alpha value a little bit more.
If I do so, what I notice is that
I'm slightly changing the length and duration of my orbit.
So in this case here,
the period is now roughly twice as much,
and we see that we have two different peaks
and we see the same thing in the trajectories
where we're now tracing a slightly more complicated orbit.
So as we increase alpha,
we're getting richer and richer behavior.
If I increase it to 1.5,
well, now, even though I'm still using a times step of 0.1,
I start seeing a really big difference between the behavior
predicted by my good integrator in red
and by Euler in black.
So let me decrease my time steps,
because I know something's up.
That's starting to look better, but not quite.
Smaller time step.
And remember, every time I decrease my time step,
I have to do more integration steps so it runs slower.
So that's why I'm not overdoing it.
All right, I'm starting to do pretty well
and I can go see what the future looks like for my system.
This is gonna take three times as long
because I'm now gonna integrate to 150 instead of 50.
And once that's done running we'll be able to see
whether or not our time step of 0.005 was really enough.
Eventually we see that we're not predicting
the state of the system accurately at all.
We have a huge error.
So even though we were doing great at first,
eventually our two trajectories have sort of diverged.
And we could keep playing with our time step
to try and see if there is a time step small enough
to get the same trajectories over a long period of time,
but I can tell you that's not gonna be the case.
The reason for that is that the system has now entered
what we call a chaotic regime
where it is sensitive to initial conditions.
So there is a buildup of error
in that first part of the integration
and eventually the two trajectories,
even though they're very similar, diverge.
To visualize that, I'll introduce
a second set of initial conditions
in my code in the top left.
So I'm gonna use random numbers
to just randomly perturb my initial condition
between zero and 0.001.
So very small perturbation on my initial conditions.
I'm gonna integrate from these conditions
using the Euler algorithm.
So I'm gonna have a series of trajectories,
a trajectory used by Euler,
which I'm now gonna call Y E two here
for a second set of Euler.
Instead of comparing the Euler integrator
to the good integrator, the O D E four five approach,
we're gonna compare our old set of initial condition in red
and the new randomly perturbed set of initial condition
in black.
We're gonna keep our three D figure here for figure two,
but I'm gonna introduce a third figure
which is gonna look at the difference
between the black and the red curve in the first plot.
To edit some of this out,
you've probably just had a weird time jump.
I was getting an error because in my code I forgotten
to transpose my second solution.
I just need, because of my comment here,
I wanted to make sure I was using
vector columns for all of them.
Okay, so we can run and now compare our solutions.
You see one already on the screen.
So now our two Euler solutions are very similar.
Again, I'm gonna go a little further in time,
and eventually they start diverging, right?
So even though they start very closely
and you might expect them to stay close forever,
because there is this sensitivity to initial conditions
in the chaotic system,
we expect eventually the two time series to diverge.
And this is what I'm looking at
with my third Figure here
is the distance between these two time series
as time goes on.
And I'm using a semi-log scale,
so a log scale in the Y axis,
and as you can see, it goes up roughly exponentially.
I say roughly because you're in a limited space
and they're moving around,
so you do expect some up and down
just due to the constrained nature of the space.
But you have this general exponential increase in distance
between the two time series.
The exponent, this exponential rate
of separation after initial condition
is what's called the Lyapunov exponent
and it is a measure of how chaotic the system might be,
how sensitive to initial condition the system might be.
So I can now run my system for much longer,
say, up time 400 instead of 200.
I expect the two time series
to eventually diverge completely.
I've used a larger time step
to be a little more accurate in my integration.
And I'm gonna go even further.
Eventually, as you'll see,
the separation between the two methods,
we see this right now in figure three, saturates
because there is just so much separation they can take
because they're constrained to a final space
just because of the way this model is specified.
And as we integrate more and more,
we keep exploring more of the space.
If you're looking at this bottom figure here, my figure two,
I keep exploring a larger fraction of space.
So this nice shape here is what's called a strange attractor
because no matter where I start in this system,
I will end up on this surface or volume
and once I'm on it, my chaotic system
will keep exploring it, exploring every point.
This has come up in discussion before.
It's what we call ergodicity, and chaos is exactly this.
It's a deterministic system.
It's not about being unpredictable or random.
It's a deterministic system
that is sensitive to initial condition.
It has this exponential separation of trajectories
and those trajectories are aperiodic and ergodic,
meaning they don't repeat themselves
and they eventually explore all of space.
So if I kept integrating for longer and longer,
I would fill up this entire surface completely.
So that's it for now, that's our first taste of chaos.
The textbook has an entire chapter of it.
There is a great book by Steven Strogatz
on non-linear dynamics if you wanna study this topic more.
We also have a course on it here at UVM.
Otherwise, I'll see you in discussion
where we can go a little deeply into
the theory and implications of chaos.
And I'll see you in the next video to wrap up
our first module on dynamical systems.
Thank you.