WEBVTT
00:00:03.670 --> 00:00:05.250
PROFESSOR: Many
mathematical models
00:00:05.250 --> 00:00:07.750
involve high order derivatives.
00:00:07.750 --> 00:00:11.760
But the MATLAB ODE
solvers only work
00:00:11.760 --> 00:00:16.149
with systems of first
order ordinary differential
00:00:16.149 --> 00:00:17.720
equations.
00:00:17.720 --> 00:00:21.230
So we have to rewrite the
models to just involve
00:00:21.230 --> 00:00:23.380
first order derivatives.
00:00:23.380 --> 00:00:28.690
Let's see how to do that with a
very simple model, the harmonic
00:00:28.690 --> 00:00:30.290
oscillator.
00:00:30.290 --> 00:00:33.220
x double prime plus x equals 0.
00:00:33.220 --> 00:00:36.270
This involves a second
order derivative.
00:00:36.270 --> 00:00:39.040
So to write it as a
first order system,
00:00:39.040 --> 00:00:42.580
we introduced the vector y.
00:00:42.580 --> 00:00:45.740
This is a vector with
two components, x,
00:00:45.740 --> 00:00:47.660
and the derivative of x.
00:00:47.660 --> 00:00:52.830
We're just changing
notation to let y
00:00:52.830 --> 00:00:56.360
have two components,
x and x prime.
00:00:56.360 --> 00:01:00.340
Then the derivative
of y is the vector
00:01:00.340 --> 00:01:03.400
with x and x double prime.
00:01:03.400 --> 00:01:06.990
So the differential
equation now becomes
00:01:06.990 --> 00:01:11.750
y2 prime plus y1 equals zero.
00:01:11.750 --> 00:01:14.910
Do you see how we've just
rewritten this differential
00:01:14.910 --> 00:01:16.180
equation.
00:01:16.180 --> 00:01:26.750
so y2 prime is playing
of x double prime?
00:01:26.750 --> 00:01:29.515
Once you've done that,
everything else is easy.
00:01:32.290 --> 00:01:41.790
The vector system is now
y1, y2 prime is y2 minus y1.
00:01:41.790 --> 00:01:45.670
The first components
says y1 prime is y2.
00:01:45.670 --> 00:01:48.850
That's just saying
that the derivative
00:01:48.850 --> 00:01:51.580
of the first component
is the second.
00:01:51.580 --> 00:01:54.070
Here's the differential
equation itself.
00:01:54.070 --> 00:02:00.300
Y2 prime is minus y1 is the
actual harmonic oscillator
00:02:00.300 --> 00:02:01.370
differential equation.
00:02:04.370 --> 00:02:08.330
When we write this as an
autonomous function for MATLAB,
00:02:08.330 --> 00:02:11.130
here's the autonomous function.
00:02:11.130 --> 00:02:15.040
f is an autonomous
function of t and y,
00:02:15.040 --> 00:02:17.800
that doesn't depend upon t.
00:02:17.800 --> 00:02:22.340
First it's a vector
now, a column vector.
00:02:22.340 --> 00:02:25.600
The first component of f is y2.
00:02:25.600 --> 00:02:30.390
And the second component is -y1.
00:02:30.390 --> 00:02:34.450
The first component here is
just a matter of notation.
00:02:34.450 --> 00:02:38.717
All the content is in the second
component, which expresses
00:02:38.717 --> 00:02:39.800
the differential equation.
00:02:43.810 --> 00:02:46.120
Now for some
initial conditions--
00:02:46.120 --> 00:02:50.600
suppose the initial conditions
are that x of 0 is 0,
00:02:50.600 --> 00:02:53.990
and x prime of 0 is 1.
00:02:53.990 --> 00:02:57.940
In terms of the vector
y, that's y1 of 0,
00:02:57.940 --> 00:03:01.340
the first component of y is 0.
00:03:01.340 --> 00:03:04.230
And the second component is 1.
00:03:04.230 --> 00:03:10.140
Or in vector terms, the
initial vector is 0, 1.
00:03:10.140 --> 00:03:15.105
That implies they solution
is sine t and cosine t.
00:03:18.130 --> 00:03:22.350
When we write the initial
condition in the MATLAB,
00:03:22.350 --> 00:03:27.480
it's the column vector 0, 1.
00:03:27.480 --> 00:03:29.870
Let's bring up the
MATLAB command window.
00:03:29.870 --> 00:03:32.050
Here's the
differential equation.
00:03:32.050 --> 00:03:34.700
y1 prime is y2.
00:03:34.700 --> 00:03:37.010
And y2 prime is -y1.
00:03:37.010 --> 00:03:40.120
Here's the harmonic oscillator.
00:03:40.120 --> 00:03:42.350
We're going to
integrate from 0 to 2pi,
00:03:42.350 --> 00:03:45.230
because they're trig functions.
00:03:45.230 --> 00:03:51.140
And I'm going to ask for output
in steps of 2 pi over 36,
00:03:51.140 --> 00:03:54.440
which corresponds
to every 10 degrees
00:03:54.440 --> 00:03:58.370
like the runways at an airport.
00:03:58.370 --> 00:04:02.410
I'm going to need an
initial condition.
00:04:02.410 --> 00:04:06.680
y0 not is 0.
00:04:06.680 --> 00:04:15.170
I need a column vector, 0,
1, for the two components.
00:04:15.170 --> 00:04:18.300
I'm going to use
ODE45, and if I call it
00:04:18.300 --> 00:04:25.876
with no output arguments, ODE45
of the differential equation f,
00:04:25.876 --> 00:04:32.970
t span the time interval,
and y0 the initial condition.
00:04:32.970 --> 00:04:37.100
If I call ODE45 with
no output arguments,
00:04:37.100 --> 00:04:40.470
it just plots the
solution automatically.
00:04:40.470 --> 00:04:45.970
And here we get a graph
of cosine t starting at 1,
00:04:45.970 --> 00:04:50.390
and sine t starting at 0.
00:04:50.390 --> 00:04:53.910
Now if I go back to
the command window,
00:04:53.910 --> 00:05:01.400
and ask to capture
the output in t and y,
00:05:01.400 --> 00:05:08.390
I then get vectors of output.
00:05:08.390 --> 00:05:15.040
37 steps, vector t,
and two components
00:05:15.040 --> 00:05:20.790
y, the two columns
containing sine and cosine.
00:05:20.790 --> 00:05:24.130
Now I can plot them
in the phase plane.
00:05:24.130 --> 00:05:28.610
Plot ya against y2.
00:05:28.610 --> 00:05:37.390
If I regularize the axes, I
get a nice plot of a circle
00:05:37.390 --> 00:05:43.970
with small circles every
10 degrees, as I said,
00:05:43.970 --> 00:05:46.035
like the runways at an airport.
00:05:48.570 --> 00:05:52.660
The Van der Pol oscillator
was introduced in 1927
00:05:52.660 --> 00:06:00.080
by Dutch electrical engineer, to
model oscillations in a circuit
00:06:00.080 --> 00:06:02.730
involving vacuum tubes.
00:06:02.730 --> 00:06:07.620
It has a nonlinear damping term.
00:06:07.620 --> 00:06:11.200
It's since been used
to model phenomena
00:06:11.200 --> 00:06:13.470
in all kinds of
fields, including
00:06:13.470 --> 00:06:16.740
geology and neurology.
00:06:16.740 --> 00:06:20.980
It exhibits chaotic behavior.
00:06:20.980 --> 00:06:24.790
We're interested in it
for numerical analysis
00:06:24.790 --> 00:06:28.860
because, as the
parameter mu increases,
00:06:28.860 --> 00:06:34.500
the problem becomes
increasingly stiff.
00:06:34.500 --> 00:06:37.100
To write it as a first
order system for use
00:06:37.100 --> 00:06:42.840
with the MATLAB ODE solvers,
we introduce the vector y,
00:06:42.840 --> 00:06:46.890
containing x and x prime.
00:06:46.890 --> 00:06:52.460
So y prime is x prime
and x double prime.
00:06:52.460 --> 00:06:54.360
And then the
differential equation
00:06:54.360 --> 00:07:02.850
is written so that the first
component of y prime is y2.
00:07:02.850 --> 00:07:05.960
And then the
differential equation
00:07:05.960 --> 00:07:10.540
is written in the
second component of y.
00:07:10.540 --> 00:07:15.730
Here's the nonlinear
damping term minus y1.
00:07:15.730 --> 00:07:22.350
When mu is 0, this becomes
the harmonic oscillator.
00:07:22.350 --> 00:07:26.655
And here it is as the
anonymous function.
00:07:32.080 --> 00:07:36.160
Let's run some experiments with
the Van der Pol oscillator.
00:07:36.160 --> 00:07:39.850
First of all, I have to
define the value of mu.
00:07:39.850 --> 00:07:41.860
And I'll pick a
modest value of mu.
00:07:41.860 --> 00:07:43.770
Mu equals 100.
00:07:43.770 --> 00:07:49.860
And now with mu set, I can
define the anonymous function.
00:07:49.860 --> 00:07:52.100
It involves a value of mu.
00:07:52.100 --> 00:07:56.690
And here is the Van
der Pol equation
00:07:56.690 --> 00:07:59.820
in the second component of f.
00:07:59.820 --> 00:08:04.350
I'm going to gather statistics
about how hard the ODE
00:08:04.350 --> 00:08:06.270
solver works.
00:08:06.270 --> 00:08:09.210
And for that, I'm
going to use ODE set,
00:08:09.210 --> 00:08:11.940
and tell it I want
to turn on stats.
00:08:16.760 --> 00:08:18.265
I need an initial condition.
00:08:23.970 --> 00:08:28.030
Now I'm going to run
ODE45 on this problem.
00:08:28.030 --> 00:08:32.740
And I'm specifying just
a starting value of t,
00:08:32.740 --> 00:08:35.990
and a final value of t.
00:08:35.990 --> 00:08:41.490
ODE45 is going to pick
its own time steps.
00:08:41.490 --> 00:08:47.280
And I know with t
final equals 460,
00:08:47.280 --> 00:08:51.570
it's going to integrate
over it about two
00:08:51.570 --> 00:08:52.735
periods of the solution.
00:08:56.450 --> 00:09:01.890
Now we can watch
it go for a while.
00:09:01.890 --> 00:09:04.380
It's taking lots of steps.
00:09:04.380 --> 00:09:07.030
And it's beginning
to slow down as it
00:09:07.030 --> 00:09:08.270
takes more and more steps.
00:09:16.440 --> 00:09:21.060
Now this begins to get painfully
slow as it runs into stiffness.
00:09:25.660 --> 00:09:27.970
I'll turn off the
camera for a while here,
00:09:27.970 --> 00:09:30.010
so you don't have to
watch all these steps.
00:09:35.130 --> 00:09:38.140
We're trying to
get up here to 460.
00:09:38.140 --> 00:09:42.880
And I'll turn it back on
as we get close to the end.
00:10:07.620 --> 00:10:10.980
OK, well, the camera's been
off about three minutes.
00:10:10.980 --> 00:10:13.920
And you can see how
far we've gotten.
00:10:13.920 --> 00:10:16.050
We're nowhere near the end.
00:10:16.050 --> 00:10:20.450
So I'm going to press
the stop button here.
00:10:20.450 --> 00:10:29.660
And we'll go back to
the command window.
00:10:29.660 --> 00:10:35.440
And oh, so we didn't
get to the end here.
00:10:38.320 --> 00:10:42.166
Let me back off on
the time interval
00:10:42.166 --> 00:10:47.309
and try this value here.
00:10:47.309 --> 00:10:48.100
See how that works.
00:10:53.320 --> 00:10:56.650
So this is going to
still take lots of steps.
00:10:56.650 --> 00:11:04.810
But we'll be able to-- This
will go over about one period.
00:11:04.810 --> 00:11:07.400
We'll actually get
to the end here.
00:11:24.020 --> 00:11:26.750
I'll leave the camera
on until we finish.
00:11:44.480 --> 00:11:47.780
OK so that took a
little under a minute.
00:11:47.780 --> 00:11:52.390
And it took nearly 15,000 steps.
00:11:52.390 --> 00:11:54.355
So let's run it
with a stiff solver.
00:12:04.150 --> 00:12:06.550
There.
00:12:06.550 --> 00:12:11.670
So it took half a second,
and only 500 steps.
00:12:11.670 --> 00:12:18.800
So there's a modest
example of stiffness here.
00:12:18.800 --> 00:12:25.420
So let's examine the
Van der Pol equation
00:12:25.420 --> 00:12:28.520
using the stiff solver.
00:12:28.520 --> 00:12:30.930
Let's capture the
output and plot it.
00:12:30.930 --> 00:12:33.890
Because that plot
wasn't very interesting.
00:12:33.890 --> 00:12:36.780
I want to plot it a
couple of different ways.
00:12:36.780 --> 00:12:42.600
And again, I want to go
back up to the-- capture
00:12:42.600 --> 00:12:43.620
a couple periods.
00:12:51.130 --> 00:12:53.450
Let's plot one of the
current components
00:12:53.450 --> 00:12:54.695
as a function of time.
00:12:57.380 --> 00:12:58.130
And here it is.
00:12:58.130 --> 00:13:00.410
Here's the Van der Pol equation.
00:13:00.410 --> 00:13:03.840
And you can see it has
an initial transient,
00:13:03.840 --> 00:13:10.770
and then it settles into
this periodic oscillation
00:13:10.770 --> 00:13:20.400
with these really
steep spikes here.
00:13:20.400 --> 00:13:26.960
And even this stiff
solver is working hard
00:13:26.960 --> 00:13:28.610
at these rapid changes.
00:13:28.610 --> 00:13:30.960
We've got a fair
number of points
00:13:30.960 --> 00:13:35.720
in here, as it is it
chooses the step size.
00:13:35.720 --> 00:13:42.550
And now, let's go back
to the command line
00:13:42.550 --> 00:13:47.502
and do a phase plane plot.
00:13:51.230 --> 00:13:55.650
So here's the phase plane
plot of this oscillator
00:13:55.650 --> 00:13:56.750
with damping.
00:13:56.750 --> 00:14:02.840
And it's nowhere near a circle,
which it would be if mu was 0.
00:14:02.840 --> 00:14:09.390
And this is the characteristic
behavior of the Van der Pol
00:14:09.390 --> 00:14:10.940
oscillator.