WEBVTT
00:00:04.550 --> 00:00:06.900
CLIVE MOLER: The
Lotka-Volterra Altera predator
00:00:06.900 --> 00:00:12.820
prey equations are the
granddaddy of all models
00:00:12.820 --> 00:00:16.950
involvement competition
between species.
00:00:16.950 --> 00:00:19.480
They are the
foundation of fields
00:00:19.480 --> 00:00:22.860
like mathematical ecology.
00:00:22.860 --> 00:00:29.440
Think of the two species
as rabbits and foxes
00:00:29.440 --> 00:00:36.850
or moose and wolves or
little fish in big fish.
00:00:36.850 --> 00:00:42.290
Y1 represents the prey,
who would live peacefully
00:00:42.290 --> 00:00:46.020
by themselves if there
were no predators.
00:00:46.020 --> 00:00:49.920
I've chosen the units
of time and population
00:00:49.920 --> 00:00:51.840
so that the
coefficients in front
00:00:51.840 --> 00:00:57.130
of the leading
linear terms are one.
00:00:57.130 --> 00:01:04.370
So y1 prime equals y1
represents exponential growth
00:01:04.370 --> 00:01:07.625
of the prey in the
absence of any predators.
00:01:10.610 --> 00:01:15.460
The predators need the
prey, live on the prey.
00:01:15.460 --> 00:01:19.500
So in the absence of any
prey, this minus sign
00:01:19.500 --> 00:01:21.090
is all important.
00:01:21.090 --> 00:01:27.930
So y2 prime equals minus y2
represents exponential decay.
00:01:27.930 --> 00:01:31.900
And the predators
die off exponentially
00:01:31.900 --> 00:01:36.970
in the absence of any prey.
00:01:36.970 --> 00:01:39.680
But then here are
the non-linear terms.
00:01:39.680 --> 00:01:45.300
These are like logistics terms,
except with the interaction
00:01:45.300 --> 00:01:47.580
between the two species.
00:01:47.580 --> 00:01:52.880
The growth of Y1 is limited
by the presence of y2.
00:01:52.880 --> 00:02:04.390
So y1 will grow until this term
becomes 0 one y2 read reaches
00:02:04.390 --> 00:02:06.440
mu2.
00:02:06.440 --> 00:02:13.510
On the other hand,
the decay of y2
00:02:13.510 --> 00:02:18.000
becomes 0 when y1 reaches mu1.
00:02:21.430 --> 00:02:25.940
To complete this specification,
we need the initial conditions.
00:02:25.940 --> 00:02:30.220
So we have two values
eta 1 and eta 2,
00:02:30.220 --> 00:02:35.150
which are the initial
values of y1 and y2.
00:02:35.150 --> 00:02:40.760
So these four parameters,
two mus and two etas,
00:02:40.760 --> 00:02:45.885
are the four parameters we have
in our predator prey model.
00:02:49.340 --> 00:02:52.820
Don't worry about the fact that
these are continuous variables
00:02:52.820 --> 00:02:57.110
and that we can have non-integer
numbers of individuals.
00:02:57.110 --> 00:03:02.580
We can have half of a rabbit
or a tenth of a moose.
00:03:02.580 --> 00:03:05.200
These are, after
all, models that
00:03:05.200 --> 00:03:09.280
are idealized versions of
what's happening in nature.
00:03:14.160 --> 00:03:19.740
The critical points are when
the derivatives become 0.
00:03:19.740 --> 00:03:21.920
There's a critical
point at the origin.
00:03:21.920 --> 00:03:26.140
But the interesting one is
when these terms become 0.
00:03:29.920 --> 00:03:32.800
So that's the point
where y1 and y2
00:03:32.800 --> 00:03:35.350
are equal to the mu1 and mu2.
00:03:38.510 --> 00:03:41.380
We have to look at the Jacobian.
00:03:41.380 --> 00:03:44.420
And here's the
Jacobian in general.
00:03:44.420 --> 00:03:49.650
And at the critical point, the
Jacobian-- here's the Jacobian.
00:03:49.650 --> 00:03:55.810
The eigenvalues of the
Jacobian are plus or minus I.
00:03:55.810 --> 00:04:02.730
And so this critical
point is a stable center
00:04:02.730 --> 00:04:05.035
with a period 2pi.
00:04:11.440 --> 00:04:13.790
So these are
non-linear equations.
00:04:13.790 --> 00:04:16.360
We can't express the
solution in terms
00:04:16.360 --> 00:04:19.640
of simple analytic functions.
00:04:19.640 --> 00:04:22.990
We have to compute
them numerically.
00:04:22.990 --> 00:04:26.370
But we do know this
about their behavior.
00:04:26.370 --> 00:04:30.820
If the initial conditions are
close to the critical point,
00:04:30.820 --> 00:04:33.410
the solution is periodic.
00:04:33.410 --> 00:04:36.740
The period is close to 2pi.
00:04:36.740 --> 00:04:41.310
And the orbit is
close to an ellipse.
00:04:41.310 --> 00:04:45.220
On the other hand, if the
initial conditions are far
00:04:45.220 --> 00:04:48.060
from the critical
point, then it turns out
00:04:48.060 --> 00:04:51.210
the solution is still periodic.
00:04:51.210 --> 00:04:56.660
But the period is greater
than 2pi and the orbit
00:04:56.660 --> 00:05:00.060
is far from an ellipse.
00:05:04.680 --> 00:05:08.990
OK, let's bring up MATLAB
and compute a solution.
00:05:08.990 --> 00:05:10.450
We need the parameters.
00:05:10.450 --> 00:05:15.410
Here's a mu and an eta.
00:05:15.410 --> 00:05:18.830
And I'll set the
differential equation.
00:05:18.830 --> 00:05:26.410
And then choose ODE 45 and
we'll integrate from 0 to 25.
00:05:26.410 --> 00:05:28.340
And here's the solution.
00:05:28.340 --> 00:05:30.080
It's periodic.
00:05:30.080 --> 00:05:32.390
A predator and a prey.
00:05:32.390 --> 00:05:36.550
And it looks like
the period, it's
00:05:36.550 --> 00:05:41.120
returning back to the initial
condition is 100 and 400.
00:05:41.120 --> 00:05:45.420
And it's returning back
to that along about here.
00:05:45.420 --> 00:05:47.420
And we've integrated
over something
00:05:47.420 --> 00:05:49.960
more than three periods.
00:05:49.960 --> 00:05:56.900
I happen to know that
the period is 6.5.
00:05:56.900 --> 00:06:02.860
And so I can-- I'm
going to set-- I want
00:06:02.860 --> 00:06:05.570
to compute to higher accuracy.
00:06:05.570 --> 00:06:08.940
10 to the minus 6.
00:06:08.940 --> 00:06:13.650
And let's capture the
solution and integrate
00:06:13.650 --> 00:06:14.840
over three periods.
00:06:23.810 --> 00:06:29.460
It generated 337 points.
00:06:29.460 --> 00:06:37.730
Let's plot it with
the finer dots.
00:06:37.730 --> 00:06:42.610
And we can see here, we've
gone over three periods
00:06:42.610 --> 00:06:49.310
and return back to our
initial values of 100 and 300.
00:06:49.310 --> 00:06:56.260
And now I'm going
to use something
00:06:56.260 --> 00:07:00.290
that'll show off the
periodicity of function
00:07:00.290 --> 00:07:01.940
in MATLAB called Comet.
00:07:15.350 --> 00:07:15.850
One orbit.
00:07:18.760 --> 00:07:21.520
Two orbits.
00:07:21.520 --> 00:07:22.265
Three orbits.
00:07:29.210 --> 00:07:32.150
Determining the period
of a periodic solution
00:07:32.150 --> 00:07:36.530
is often the important
part of a calculation.
00:07:36.530 --> 00:07:42.190
In the MATLAB ODE suite, this
is done with an event handler.
00:07:42.190 --> 00:07:51.790
So I'm going to use ODE's set to
provide an event handler called
00:07:51.790 --> 00:07:53.260
Pit Stop.
00:07:53.260 --> 00:07:57.960
And here's the code in Pit Stop.
00:07:57.960 --> 00:08:06.630
This code is called every
time step in the integration.
00:08:06.630 --> 00:08:10.310
And I'm not going to
go into detail here,
00:08:10.310 --> 00:08:15.160
but it computes a
function g that we're
00:08:15.160 --> 00:08:17.530
looking to see when this is 0.
00:08:17.530 --> 00:08:22.660
And when g returns to 0,
the integration is stopped.
00:08:22.660 --> 00:08:28.380
I'm going to-- I have a
display function in here.
00:08:28.380 --> 00:08:29.990
Ordinarily, this
wouldn't be here.
00:08:29.990 --> 00:08:31.760
But I want to look
at this to see
00:08:31.760 --> 00:08:39.030
how the calculation is done.
00:08:39.030 --> 00:08:41.760
It says, terminate when
y returns to the point
00:08:41.760 --> 00:08:45.770
where its angle mu is the
same as the angle between eta
00:08:45.770 --> 00:08:46.990
and mu.
00:08:46.990 --> 00:08:52.440
This is more reliable
than just finding
00:08:52.440 --> 00:08:59.500
when the solution returns
back to its initial condition.
00:08:59.500 --> 00:09:03.820
So let's run ODE 45 and
tell it to integrate over
00:09:03.820 --> 00:09:06.922
an infinite time step,
over an infinite interval.
00:09:06.922 --> 00:09:08.630
That's not going to
happen, because we're
00:09:08.630 --> 00:09:13.710
going to stop as soon as g is 0
but with the option vector set
00:09:13.710 --> 00:09:16.260
with this event handler.
00:09:16.260 --> 00:09:17.830
So here it is.
00:09:17.830 --> 00:09:23.350
And here's the output
produced by Pit Stop
00:09:23.350 --> 00:09:25.630
as we integrate along.
00:09:25.630 --> 00:09:30.130
Here's the values
that it's found.
00:09:30.130 --> 00:09:36.090
And here's where g returns to 0.
00:09:36.090 --> 00:09:41.940
I've produced a t vector
with 117 values in it.
00:09:41.940 --> 00:09:49.310
And the final value
of t is this value
00:09:49.310 --> 00:09:54.530
that I said was the
period of this solution.
00:09:54.530 --> 00:09:57.720
So that's how the
period was determined
00:09:57.720 --> 00:10:02.805
by this event handling
feature in the ODE solvers.
00:10:06.110 --> 00:10:08.410
Great.
00:10:08.410 --> 00:10:12.260
I have a program called
Predator Prey that's
00:10:12.260 --> 00:10:18.820
in the collection of programs
that comes with NCM, Numerical
00:10:18.820 --> 00:10:21.160
Computing with MATLAB.
00:10:21.160 --> 00:10:25.375
My book that's available
on the MathWorks website.
00:10:29.610 --> 00:10:38.830
Predator prey offers this
graphic user interface
00:10:38.830 --> 00:10:41.810
to demonstrate what we've been
talking about the predator prey
00:10:41.810 --> 00:10:44.660
equations.
00:10:44.660 --> 00:10:48.950
The top display shows the
phase plane plot, the plot
00:10:48.950 --> 00:10:52.690
of prey versus predator.
00:10:52.690 --> 00:10:58.840
And the bottom display shows
the time series plot, the plot
00:10:58.840 --> 00:11:00.870
of the two populations.
00:11:00.870 --> 00:11:02.920
And initially, it's
set at the conditions
00:11:02.920 --> 00:11:04.880
we've been talking about.
00:11:04.880 --> 00:11:12.540
Here's 400 rabbits and 100
foxes around the critical point
00:11:12.540 --> 00:11:17.370
of 300 rabbits and 200 foxes.
00:11:17.370 --> 00:11:22.110
And here's the period
of 6.5 some odd
00:11:22.110 --> 00:11:24.250
that we've been working with.
00:11:24.250 --> 00:11:27.510
Now, it says drag either dot.
00:11:27.510 --> 00:11:34.530
And you can change
the initial conditions
00:11:34.530 --> 00:11:37.380
or the critical point.
00:11:37.380 --> 00:11:40.280
If I bring the initial
conditions close
00:11:40.280 --> 00:11:46.200
to the critical point,
then the phase plane plot
00:11:46.200 --> 00:11:48.480
becomes close to ellipse.
00:11:48.480 --> 00:11:53.760
And the period
becomes close to 2pi.
00:11:53.760 --> 00:12:00.590
This is 6.28, which
is twice 3.14.
00:12:00.590 --> 00:12:08.390
But if I change it so
that the two are far away
00:12:08.390 --> 00:12:13.310
from each other, then
the phase plane plot
00:12:13.310 --> 00:12:21.250
becomes quite different
from an ellipse.
00:12:21.250 --> 00:12:22.550
It's always periodic.
00:12:22.550 --> 00:12:28.400
That's an amazing thing about
these nonlinear equations.
00:12:28.400 --> 00:12:32.090
They always have a
periodic solution.
00:12:32.090 --> 00:12:38.490
But now as you can see, the
period is greater than 2pi.
00:12:38.490 --> 00:12:47.340
And the solutions are
very far from sinusoids.
00:12:47.340 --> 00:12:51.860
So that's the pred
prey app that's
00:12:51.860 --> 00:12:59.520
available with the programs
in from the mathworks website
00:12:59.520 --> 00:13:04.705
under NCM for Numerical
Computing with MATLAB.
00:13:10.910 --> 00:13:13.780
Whoops, I stand corrected.
00:13:13.780 --> 00:13:20.250
If you Google Moler predprey,
it tries to talk me out of it,
00:13:20.250 --> 00:13:23.570
but then it shows it's in
my second book, Experiments
00:13:23.570 --> 00:13:24.500
with MATLAB.
00:13:24.500 --> 00:13:26.000
There are two books.
00:13:26.000 --> 00:13:30.410
Numerical Computing with MATLAB
and Experiments with MATLAB.
00:13:30.410 --> 00:13:35.190
And pred prey is in the second
one, Experiments with MATLAB.
00:13:35.190 --> 00:13:37.640
You can go to
website and you can
00:13:37.640 --> 00:13:42.890
download all of the
programs from EXM
00:13:42.890 --> 00:13:50.830
or you can go down here
and here's pred prey.
00:13:50.830 --> 00:13:59.001
So it's in Experiments with
MATLAB on the MathWorks
00:13:59.001 --> 00:13:59.500
website.
00:14:05.390 --> 00:14:10.980
Just Google Moler predator
prey and you'll find it.