WEBVTT
00:00:01.540 --> 00:00:03.910
The following content is
provided under a Creative
00:00:03.910 --> 00:00:05.300
Commons license.
00:00:05.300 --> 00:00:07.510
Your support will help
MIT OpenCourseWare
00:00:07.510 --> 00:00:11.600
continue to offer high quality
educational resources for free.
00:00:11.600 --> 00:00:14.140
To make a donation or to
view additional materials
00:00:14.140 --> 00:00:18.100
from hundreds of MIT courses,
visit MIT OpenCourseWare
00:00:18.100 --> 00:00:19.320
at ocw.mit.edu.
00:00:23.289 --> 00:00:24.830
WILLIAM GREEN JR:
All right, so we're
00:00:24.830 --> 00:00:33.260
going to start a new topic today
about ordinary differential
00:00:33.260 --> 00:00:35.555
equations initial
value problems.
00:00:35.555 --> 00:00:39.020
But before I start
to talk about that,
00:00:39.020 --> 00:00:41.840
I want to remind you next
week's schedule is weird.
00:00:41.840 --> 00:00:45.830
So there's no classes on
Monday, but instead we'll
00:00:45.830 --> 00:00:48.860
have Monday class on Tuesday.
00:00:48.860 --> 00:00:51.950
So I'll see you next
on Tuesday morning.
00:00:51.950 --> 00:00:55.980
And then we're not going to
have a class on Wednesday,
00:00:55.980 --> 00:00:59.490
but we have the quiz
on Wednesday night.
00:00:59.490 --> 00:01:02.250
And then we'll get back
to normal on Friday.
00:01:02.250 --> 00:01:04.709
And I would expect that
possibly the TAs might
00:01:04.709 --> 00:01:08.770
be interested in giving a help
session or review on Wednesday
00:01:08.770 --> 00:01:10.137
the class period.
00:01:10.137 --> 00:01:11.470
AUDIENCE: We'll be here at 3:00.
00:01:11.470 --> 00:01:11.910
WILLIAM GREEN JR:
They'll be here.
00:01:11.910 --> 00:01:13.390
So if you want to
come at class period,
00:01:13.390 --> 00:01:14.764
and just ask
questions and stuff.
00:01:14.764 --> 00:01:18.140
AUDIENCE: And I'll have your
PSAT graded back, if you like.
00:01:18.140 --> 00:01:19.750
WILLIAM GREEN JR:
All right, got that?
00:01:19.750 --> 00:01:22.400
PSAT be ready as well.
00:01:22.400 --> 00:01:26.850
OK, so I'll talk
about today ODE-IVPs.
00:01:26.850 --> 00:01:31.080
AUDIENCE: Yeah, we'll still have
office hours Monday [INAUDIBLE]
00:01:31.080 --> 00:01:33.170
WILLIAM GREEN JR: Can
you guys hear this?
00:01:33.170 --> 00:01:35.210
Yes, office hours
like normal Monday.
00:01:38.047 --> 00:01:40.130
And for those of you who
want to get your homework
00:01:40.130 --> 00:01:41.921
started early, we might
even post homework.
00:01:54.290 --> 00:02:01.632
And I'm going to focus
exclusively on ODE-IVPs
00:02:01.632 --> 00:02:02.590
you can write this way.
00:02:11.600 --> 00:02:16.807
All right, and I
would just comment,
00:02:16.807 --> 00:02:17.890
what do these things mean?
00:02:17.890 --> 00:02:20.060
ODE means ordinary
differential equation.
00:02:20.060 --> 00:02:22.856
It means that the only
differentials in the problem
00:02:22.856 --> 00:02:24.230
are with respect
to one variable,
00:02:24.230 --> 00:02:27.200
so I'm going to call that t.
00:02:27.200 --> 00:02:29.070
It can be z or x or
whatever you want,
00:02:29.070 --> 00:02:31.190
but I'm going to call it t here.
00:02:31.190 --> 00:02:34.220
And so this is not partial
differential equations
00:02:34.220 --> 00:02:38.030
where you have differentials
respect to many variables.
00:02:38.030 --> 00:02:44.460
Initial value problems mean
that all the initial conditions
00:02:44.460 --> 00:02:46.780
are given in at a
single t point where
00:02:46.780 --> 00:02:50.402
t is that the variable that
appears in the differential.
00:02:53.100 --> 00:02:55.100
All right so this is the
form we're going to do.
00:02:55.100 --> 00:02:59.880
This is first order a
way of writing the ODEs.
00:02:59.880 --> 00:03:03.965
If you recall from homework
zero problem number one,
00:03:03.965 --> 00:03:06.590
we showed you there how you can
write higher order differential
00:03:06.590 --> 00:03:07.250
equations--
00:03:10.580 --> 00:03:12.650
convert them into the
form of first order.
00:03:16.100 --> 00:03:21.134
And I guess that's
worthwhile to point down.
00:03:21.134 --> 00:03:23.050
I would also comment
that I going to just talk
00:03:23.050 --> 00:03:28.170
about f of y, however
the ODE solvers in Matlab
00:03:28.170 --> 00:03:34.600
expect f of the t
comma y as their from.
00:03:34.600 --> 00:03:42.820
All right, so let's just explain
all the details about this.
00:03:42.820 --> 00:03:47.900
So suppose I had a differential
equation like this.
00:03:58.887 --> 00:04:02.080
All right, so this is like a
differential equation you'd
00:04:02.080 --> 00:04:06.260
have for a body move moving--
00:04:06.260 --> 00:04:09.520
a particle moving with
some friction force on it
00:04:09.520 --> 00:04:11.760
and some time varying force.
00:04:11.760 --> 00:04:12.378
Yeah?
00:04:12.378 --> 00:04:14.170
AUDIENCE: Are you going to
be posting your notes online?
00:04:14.170 --> 00:04:15.310
WILLIAM GREEN JR: No,
I don't have any notes.
00:04:15.310 --> 00:04:16.268
This is it.
00:04:20.579 --> 00:04:26.180
All right, so this is a
differential equation.
00:04:26.180 --> 00:04:29.420
What you can do is to convert
it into the form I have above,
00:04:29.420 --> 00:04:35.100
is to write a new variable, say
v, that's defined to be dx dt.
00:04:35.100 --> 00:04:38.250
That's the velocity.
00:04:38.250 --> 00:04:41.344
And has a time
dependence, but I told you
00:04:41.344 --> 00:04:43.010
I don't like time
dependence, so I'm not
00:04:43.010 --> 00:04:44.640
going to write it that way.
00:04:44.640 --> 00:04:47.340
So we can also define
another variable.
00:04:47.340 --> 00:04:57.950
So we can define our new y
vector to be x, v, and t.
00:05:01.110 --> 00:05:02.850
So this is y-- the
first component of y,
00:05:02.850 --> 00:05:04.266
this is the second
component of y,
00:05:04.266 --> 00:05:06.210
this is the third
component of y.
00:05:06.210 --> 00:05:10.110
And then I can write
down what f of y is.
00:05:22.000 --> 00:05:27.190
So the equation for dx dt is y2.
00:05:27.190 --> 00:05:32.260
It's the second component
of y, dx dt is equal to v.
00:05:32.260 --> 00:05:38.750
The equation for the dv dt
is given by this equation.
00:05:38.750 --> 00:05:53.442
So that's f of y3
over m minus gamma y2.
00:05:53.442 --> 00:05:57.600
And the equation
dt dt is just one.
00:06:00.390 --> 00:06:03.360
So this is my f of y.
00:06:03.360 --> 00:06:05.780
So this is f1, f2, f3.
00:06:05.780 --> 00:06:06.980
Is that all right?
00:06:10.114 --> 00:06:11.655
So this is just how
you can convert--
00:06:11.655 --> 00:06:13.488
if somebody writes you
an equation this way,
00:06:13.488 --> 00:06:16.520
and you want to get it into the
standard form that I'm going
00:06:16.520 --> 00:06:19.700
to show you here--
00:06:19.700 --> 00:06:21.980
dy dt equals f of y--
00:06:21.980 --> 00:06:22.980
this is how you convert.
00:06:25.710 --> 00:06:29.747
And a skill we haven't
really emphasized so far
00:06:29.747 --> 00:06:31.830
in the course, but I think
a very important skill,
00:06:31.830 --> 00:06:35.310
is to be able to get like a
chemical engineering problem
00:06:35.310 --> 00:06:39.630
and rearrange it so it
becomes a standard form
00:06:39.630 --> 00:06:43.930
that maps up to one of the
solvers you know how to use.
00:06:43.930 --> 00:06:46.260
So you've got a problem has
some differential in it.
00:06:46.260 --> 00:06:48.450
Right away you actually,
can I somehow rewrite it
00:06:48.450 --> 00:06:51.710
so it looks like this standard
form-- or actually for Matlab,
00:06:51.710 --> 00:06:52.480
this form.
00:06:52.480 --> 00:06:58.800
So Matlab will say f of t, y.
00:06:58.800 --> 00:07:02.070
If you can do that, then you can
go ahead and use those Matlab
00:07:02.070 --> 00:07:04.230
solvers to solve it.
00:07:04.230 --> 00:07:08.141
But your variables may not have
a letter y appearing anywhere,
00:07:08.141 --> 00:07:08.640
right?
00:07:08.640 --> 00:07:09.570
You have to figure
out how to write it
00:07:09.570 --> 00:07:10.980
into the standard form.
00:07:10.980 --> 00:07:14.190
And then once you can do
that, you can use the solvers.
00:07:14.190 --> 00:07:16.110
And the same thing
with the non-linear
00:07:16.110 --> 00:07:18.700
equations solvers, with
the optimization software.
00:07:18.700 --> 00:07:20.580
If you can rewrite
it into the form that
00:07:20.580 --> 00:07:23.310
matches the standard form, then
once you're done, you're done.
00:07:23.310 --> 00:07:26.220
Then you can just go and
call the canned programs
00:07:26.220 --> 00:07:27.780
that work well for those forms.
00:07:32.130 --> 00:07:34.160
All right, now what do
we want as the output?
00:07:34.160 --> 00:07:36.900
So this is the problem.
00:07:36.900 --> 00:07:38.840
This is the problem as proposed.
00:07:38.840 --> 00:07:41.150
What do we really
want out of this?
00:07:41.150 --> 00:07:47.190
What we're going to get out
of it is a bunch of points y--
00:07:47.190 --> 00:07:48.260
actually a matrix, right?
00:07:48.260 --> 00:07:50.240
It would be a
vector of y-values.
00:07:50.240 --> 00:07:51.710
Typically y is a vector.
00:07:51.710 --> 00:07:55.220
And we want to know it
a bunch of time points.
00:07:55.220 --> 00:07:56.700
And we would a
lot of time points
00:07:56.700 --> 00:07:59.110
so we're getting pretty plots
of how the components of y
00:07:59.110 --> 00:08:01.404
vary with time.
00:08:01.404 --> 00:08:02.570
That's what we usually want.
00:08:02.570 --> 00:08:07.480
Sometimes we just want the
final value, y time value.
00:08:07.480 --> 00:08:09.770
We integrate to some point
and then stop, and just
00:08:09.770 --> 00:08:10.550
run through the final value.
00:08:10.550 --> 00:08:12.440
But a lot of times, we actually
want to make pretty plots,
00:08:12.440 --> 00:08:15.250
so we really want a lot of
points here and making plots.
00:08:15.250 --> 00:08:16.910
Now how many points do you
need to make a nice plot?
00:08:16.910 --> 00:08:18.701
Maybe, I don't know,
100 points, 50 points,
00:08:18.701 --> 00:08:20.000
you can make a nice plot.
00:08:20.000 --> 00:08:22.520
So that's the kind
numbers you'd like.
00:08:22.520 --> 00:08:26.270
And so the output
of these programs
00:08:26.270 --> 00:08:30.530
is going to be something like
a vector of the time points,
00:08:30.530 --> 00:08:40.380
and then a matrix of the
y-values corresponding
00:08:40.380 --> 00:08:41.750
to each time point.
00:08:41.750 --> 00:08:44.510
And the way it does in
that lab, it's reasonable.
00:08:44.510 --> 00:08:47.870
Each row of y are
the y-values that
00:08:47.870 --> 00:08:52.050
correspond to the same time
point, the first time point.
00:08:52.050 --> 00:08:54.470
So the first row of y
corresponds to the first number
00:08:54.470 --> 00:08:56.000
in the time vector.
00:08:56.000 --> 00:08:57.830
The second run is the second.
00:08:57.830 --> 00:09:01.400
Last one will be at time final.
00:09:01.400 --> 00:09:03.692
And then you can plot them.
00:09:03.692 --> 00:09:04.790
OK so far?
00:09:08.940 --> 00:09:16.940
And ideally, we would like
these y-values that we compute,
00:09:16.940 --> 00:09:19.280
these y's we get
out, ideally one
00:09:19.280 --> 00:09:21.740
that really close to
what the true solution
00:09:21.740 --> 00:09:23.914
of the true differential
equation is.
00:09:23.914 --> 00:09:25.830
But we're doing numerical
stuff, so it's never
00:09:25.830 --> 00:09:27.810
going to be exactly the same.
00:09:27.810 --> 00:09:30.340
And so a big part of the
effort is trying to figure out,
00:09:30.340 --> 00:09:33.155
how can we get the thing
to compute y calculated
00:09:33.155 --> 00:09:34.530
at each time point
to be as close
00:09:34.530 --> 00:09:36.750
as possible to the truth,
of what the true solution
00:09:36.750 --> 00:09:39.390
of the equations are?
00:09:39.390 --> 00:09:41.580
And that turns out
to be difficult,
00:09:41.580 --> 00:09:44.010
so that's a very important
thing to worry about.
00:09:47.237 --> 00:09:48.570
What else can we say about this?
00:09:52.100 --> 00:09:57.280
We specify this from t naught,
and you typically have to input
00:09:57.280 --> 00:10:00.610
your t final, how far you want
to go away from your t naught
00:10:00.610 --> 00:10:04.000
as another input.
00:10:04.000 --> 00:10:07.870
Normally, people write it where
t naught is less than t final,
00:10:07.870 --> 00:10:10.540
and so you're like integrating
from left to right.
00:10:10.540 --> 00:10:12.970
But you could actually write
them, put a negative sign
00:10:12.970 --> 00:10:14.950
equation which corresponds
to a negative t,
00:10:14.950 --> 00:10:17.140
change your variable
to of t prime that
00:10:17.140 --> 00:10:20.404
was equal to t final
minus t if you wanted to,
00:10:20.404 --> 00:10:21.820
and integrate the
other direction.
00:10:21.820 --> 00:10:26.230
And in fact in homework zero
problem one, we did that too.
00:10:26.230 --> 00:10:27.580
So you can integrate frontwards.
00:10:27.580 --> 00:10:28.867
You can integrate backwards.
00:10:28.867 --> 00:10:31.200
If somebody tells you the
initial condition or condition
00:10:31.200 --> 00:10:33.810
in the center of the range,
you integrate forward
00:10:33.810 --> 00:10:36.270
for part of the range, and
backwards the other way,
00:10:36.270 --> 00:10:37.224
so anything like this.
00:10:37.224 --> 00:10:38.640
So what's important
for this to be
00:10:38.640 --> 00:10:40.140
what's called an initial
value problem is just
00:10:40.140 --> 00:10:41.940
that all of the
specifications of y
00:10:41.940 --> 00:10:43.780
are given at the
same value of t.
00:10:43.780 --> 00:10:46.480
It doesn't matter exactly
what value of t it is.
00:10:46.480 --> 00:10:47.980
Later, we'll talk
about what happens
00:10:47.980 --> 00:10:50.005
if you don't know all
the specifications at one
00:10:50.005 --> 00:10:52.470
value of time.
00:10:52.470 --> 00:11:01.250
All right, now how
would you approach this?
00:11:01.250 --> 00:11:04.240
Probably, a bunch of you've done
this already as undergraduates,
00:11:04.240 --> 00:11:05.530
or maybe even high school.
00:11:05.530 --> 00:11:08.690
You can write like
a Taylor expansion,
00:11:08.690 --> 00:11:12.070
y of t plus delta t--
00:11:12.070 --> 00:11:16.550
is y of t plus delta t--
00:11:16.550 --> 00:11:20.261
times dy dt.
00:11:23.560 --> 00:11:25.570
You can write that.
00:11:25.570 --> 00:11:27.040
And then I say,
wow, I know dy dt.
00:11:27.040 --> 00:11:30.045
That's actually f, so this
is actually equal to--
00:11:47.750 --> 00:11:48.910
right?
00:11:48.910 --> 00:11:51.140
And then if I truncate
the Taylor expansion,
00:11:51.140 --> 00:11:52.390
now I have an update formula.
00:11:52.390 --> 00:11:55.600
So if I just forget
this last bit,
00:11:55.600 --> 00:11:58.880
I can put in my current value of
y, like the initial condition.
00:11:58.880 --> 00:12:01.000
The initial condition
here, evaluate f,
00:12:01.000 --> 00:12:03.800
multiply by delta t,
add them together,
00:12:03.800 --> 00:12:06.460
and I'll get a new value of y.
00:12:06.460 --> 00:12:09.300
And then I can repeat
over and over again.
00:12:09.300 --> 00:12:11.120
And in fact, this is
a very famous method.
00:12:11.120 --> 00:12:12.661
It's called the
Forward Euler method.
00:12:12.661 --> 00:12:15.330
Euler was a pretty smart guy, so
it's not completely ridiculous,
00:12:15.330 --> 00:12:17.121
but as I'll show you,
it has some problems.
00:12:21.000 --> 00:12:25.850
So the Forward Euler, another
way to write it is y new
00:12:25.850 --> 00:12:37.026
is equal to y old plus
delta t times f of y old.
00:12:42.510 --> 00:12:44.010
It's called Forward Euler.
00:12:53.090 --> 00:13:01.030
And we can write a little
implementation of this
00:13:01.030 --> 00:13:01.700
pretty easily.
00:13:01.700 --> 00:13:05.707
So you can write y
is equal to y naught,
00:13:05.707 --> 00:13:13.625
t is equal to t
naught for i equals 1
00:13:13.625 --> 00:13:22.540
to the number of
steps, y is equal to y
00:13:22.540 --> 00:13:30.082
plus delta t times f of y, t
is equal to t plus delta t.
00:13:30.082 --> 00:13:31.873
Somewhere up here, I
need to write delta t.
00:13:41.540 --> 00:13:44.765
OK, so there's your code
for doing Forward Euler.
00:13:44.765 --> 00:13:47.700
How many of you have
done this before?
00:13:47.700 --> 00:13:49.630
OK, so I can skip
over this very fast.
00:13:52.630 --> 00:13:55.236
How many did this, and it
didn't work for a problem?
00:13:58.110 --> 00:14:00.170
OK, so you're not doing
hard enough problems.
00:14:00.170 --> 00:14:03.710
Well, we'll correct that.
00:14:03.710 --> 00:14:06.800
All right, so I noticed
that this closely
00:14:06.800 --> 00:14:11.350
resembles the rectangle rule.
00:14:11.350 --> 00:14:16.000
So the rectangle rule
would look a lot the same.
00:14:16.000 --> 00:14:20.830
However, if you're
integrating a function--
00:14:20.830 --> 00:14:29.110
so I have I is equal to
the integral of f of t dt,
00:14:29.110 --> 00:14:31.060
then I can notice
that the derivative dI
00:14:31.060 --> 00:14:36.100
dt is equal to f of t.
00:14:39.900 --> 00:14:42.995
And so if I was going to
do the rectangle rule,
00:14:42.995 --> 00:14:45.578
it would look just the same as
this, except this would be a t.
00:14:48.330 --> 00:14:51.410
Maybe I would pick y0 to be 0.
00:14:51.410 --> 00:14:55.250
And the final value
of y I would get to
00:14:55.250 --> 00:14:59.640
would be my integral,
my i value that I want.
00:14:59.640 --> 00:15:02.490
So how many of you have
done rectangle rule?
00:15:02.490 --> 00:15:04.020
Yes?
00:15:04.020 --> 00:15:06.600
All right, so you
did this already.
00:15:06.600 --> 00:15:09.180
Now you remember in high school,
when they took the rectangle
00:15:09.180 --> 00:15:14.639
rule, they probably mentioned
that it's not very accurate.
00:15:14.639 --> 00:15:16.680
And so then they also
taught you some other ones.
00:15:16.680 --> 00:15:19.496
You guys, how many of you
learned the trapezoid rule?
00:15:19.496 --> 00:15:23.040
How about the midpoint rule?
00:15:23.040 --> 00:15:27.050
How about Simpson's rule?
00:15:27.050 --> 00:15:30.380
All right, so at least a
high school math teachers
00:15:30.380 --> 00:15:36.661
thought that this is
inadequate for doing real work,
00:15:36.661 --> 00:15:38.660
so they tell you all these
other things instead.
00:15:38.660 --> 00:15:42.020
This was like your
first baby thing.
00:15:42.020 --> 00:15:43.760
So let's talk
about why was that.
00:15:59.780 --> 00:16:03.870
So let's do the Taylor
expansion a little bit further.
00:16:03.870 --> 00:16:08.650
So let's do it for the
numeric [INAUDIBLE] case.
00:16:08.650 --> 00:16:17.670
So y t plus delta t is equal
to y of t plus delta t times
00:16:17.670 --> 00:16:26.910
f of t plus 1/2 delta t squared
times the second derivative.
00:16:32.320 --> 00:16:34.510
And so when we made
the approximation
00:16:34.510 --> 00:16:37.000
to do the rectangle rule, we
just threw this term away.
00:16:37.000 --> 00:16:39.400
So that's a leading
term that we threw away.
00:16:39.400 --> 00:16:43.280
So we made an error
order or delta t squared.
00:16:43.280 --> 00:16:46.330
Because normally the second
derivative is not exactly zero.
00:16:46.330 --> 00:16:57.330
So we have an error
in each step that's
00:16:57.330 --> 00:17:00.941
the order of delta t squared.
00:17:00.941 --> 00:17:02.190
But how many steps do we make?
00:17:02.190 --> 00:17:07.225
The number of steps is
equal to t final minus t
00:17:07.225 --> 00:17:08.550
naught over delta t.
00:17:11.270 --> 00:17:16.660
So therefore, we're making
one over delta t order steps.
00:17:16.660 --> 00:17:19.089
So the total error is going
to be approximately the number
00:17:19.089 --> 00:17:21.700
of steps times how much
error you making each step.
00:17:21.700 --> 00:17:24.099
So it's going to be something
the second power in delta t
00:17:24.099 --> 00:17:26.223
divided by something to
the first power of delta t.
00:17:26.223 --> 00:17:33.730
So the total error
in the integral I
00:17:33.730 --> 00:17:39.080
is going to be older of delta
t if you do the rectangle rule.
00:17:39.080 --> 00:17:41.730
And so order delta
t is not very good.
00:17:41.730 --> 00:17:45.870
So you want to increase
your accuracy--
00:17:45.870 --> 00:17:47.550
if you cut your
step size in half,
00:17:47.550 --> 00:17:49.050
that means you
double your CPU time,
00:17:49.050 --> 00:17:51.870
because you have to do twice
as many function evaluations.
00:17:51.870 --> 00:17:55.710
But you only gain-- you reduce
your error by a factor of two.
00:17:55.710 --> 00:17:59.510
So if you want to, say, gain
three significant figures
00:17:59.510 --> 00:18:01.200
of accuracy in your
number, then you
00:18:01.200 --> 00:18:06.770
have to do 1,000 times much
work as you did before.
00:18:06.770 --> 00:18:08.640
So I have some amount
of some precision
00:18:08.640 --> 00:18:09.979
with some number of steps.
00:18:09.979 --> 00:18:12.270
I want to increase, get three
more significant figures.
00:18:12.270 --> 00:18:14.520
I'll have to do 1,000
times as much work.
00:18:14.520 --> 00:18:16.504
That's kind of bad scaling.
00:18:16.504 --> 00:18:18.420
If I want to get six
more significant figures,
00:18:18.420 --> 00:18:20.630
I've got to do a
million times more work.
00:18:20.630 --> 00:18:25.000
So at some point, this going
to be kind of expensive.
00:18:25.000 --> 00:18:31.730
So if you contrast it, you can
do things like trapezoid rule.
00:18:31.730 --> 00:18:35.580
And as you're probably--
00:18:35.580 --> 00:18:38.193
well let's talk about
trapezoid rule for a bit.
00:18:48.030 --> 00:18:52.140
So the idea of this is
you have your function,
00:18:52.140 --> 00:18:56.340
you have some point
t, and t plus delta t.
00:18:56.340 --> 00:18:59.420
You're trying to
integrate between them.
00:18:59.420 --> 00:19:03.384
Let's say t naught, t
naught plus delta t.
00:19:03.384 --> 00:19:04.960
Try to integrate.
00:19:04.960 --> 00:19:07.510
Here's the value of
the function at t.
00:19:07.510 --> 00:19:13.301
Here is the value of the
function that t plus delta t.
00:19:13.301 --> 00:19:15.050
If you do the rectangle
rule, what you say
00:19:15.050 --> 00:19:20.125
is I just draw a line here
and integrate that rectangle.
00:19:20.125 --> 00:19:24.480
But say the real function
really looks like this.
00:19:24.480 --> 00:19:27.560
Then I missing that area
there, so that's the error.
00:19:27.560 --> 00:19:30.280
That's why the error's so big.
00:19:30.280 --> 00:19:37.140
If I do that trapezoid
rule, then the error--
00:19:37.140 --> 00:19:40.570
I actually overestimate
the value here,
00:19:40.570 --> 00:19:42.130
but the integral
of the difference
00:19:42.130 --> 00:19:43.930
between this dotted
line and the solid line
00:19:43.930 --> 00:19:45.990
is less than the integral
between the solid line
00:19:45.990 --> 00:19:48.637
and that dotted line,
so its more accurate.
00:19:48.637 --> 00:19:50.470
So that's why people
like the trapezoid rule
00:19:50.470 --> 00:19:52.886
better than the rectangle rule.
00:19:52.886 --> 00:19:55.270
If you did the
midpoint rule, you'd
00:19:55.270 --> 00:20:00.320
evaluate this function halfway
between these two points
00:20:00.320 --> 00:20:02.815
and then draw a dotted line
here and a dotted line there,
00:20:02.815 --> 00:20:05.020
and integrate that guy.
00:20:05.020 --> 00:20:07.144
And it would
overestimate for part,
00:20:07.144 --> 00:20:09.310
and underestimate for part,
and they would partially
00:20:09.310 --> 00:20:10.940
cancel each other out.
00:20:10.940 --> 00:20:13.589
And so that also would be
more accurate than doing
00:20:13.589 --> 00:20:14.380
the rectangle rule.
00:20:17.745 --> 00:20:18.620
So you did it before.
00:20:23.740 --> 00:20:25.660
And so what we're
going to try to do
00:20:25.660 --> 00:20:27.430
is, instead of
using this formula
00:20:27.430 --> 00:20:31.870
we did before, we're going to
replace this with something
00:20:31.870 --> 00:20:35.300
else that I'm going to call g.
00:20:35.300 --> 00:20:37.040
And g is something
that's supposed
00:20:37.040 --> 00:20:41.510
to be the average value
of f over the step.
00:20:41.510 --> 00:20:45.410
Instead of just taking the
value f at the starting point,
00:20:45.410 --> 00:20:47.239
I want to get sort
of the average.
00:20:47.239 --> 00:20:48.780
Right, because an
integral is related
00:20:48.780 --> 00:20:51.170
to an average over the delta t.
00:20:51.170 --> 00:20:54.500
So g is going to be my way
to estimate the average.
00:20:54.500 --> 00:20:56.429
And there's going to
be a zillion update
00:20:56.429 --> 00:20:58.220
formulas that different
mathematicians have
00:20:58.220 --> 00:20:59.626
proposed over the years.
00:20:59.626 --> 00:21:02.250
They give you different g's that
tell you how to do the update.
00:21:05.240 --> 00:21:07.500
And the key is try to find
a g that's really accurate.
00:21:07.500 --> 00:21:09.420
You want to have one that's
the most accurate you can,
00:21:09.420 --> 00:21:11.060
because then your
errors will be less.
00:21:11.060 --> 00:21:13.910
And if your errors are
less, than your delta t can
00:21:13.910 --> 00:21:17.250
be bigger, which means your
number of steps will be less,
00:21:17.250 --> 00:21:19.700
which means your CPU
time will be less.
00:21:19.700 --> 00:21:21.920
So you might be able
to get more accuracy
00:21:21.920 --> 00:21:24.367
and use less CPU
time both if you
00:21:24.367 --> 00:21:26.575
can find a really good
formula for g, for the update.
00:21:29.830 --> 00:21:32.670
So what's the g for
the trapezoid rule?
00:21:37.300 --> 00:21:46.190
It's f of t plus f of
t plus delta t over 2.
00:21:46.190 --> 00:21:48.690
That was the g you used when
you did the trapezoid rule when
00:21:48.690 --> 00:21:51.440
you were kids.
00:21:51.440 --> 00:21:53.010
What is the g for
the midpoint rule?
00:21:57.280 --> 00:22:02.382
It's f of t plus delta t over 2.
00:22:06.410 --> 00:22:08.640
And so the different g's
are called different update
00:22:08.640 --> 00:22:10.960
formulas.
00:22:10.960 --> 00:22:14.740
And these particular
ones reduce the error
00:22:14.740 --> 00:22:20.770
from delta t squared for
[INAUDIBLE] step to, I think,
00:22:20.770 --> 00:22:23.260
delta t to the cubed power.
00:22:23.260 --> 00:22:27.670
And then when you multiply
by this 1 every delta t,
00:22:27.670 --> 00:22:30.980
it turns out you go from
being order of delta t
00:22:30.980 --> 00:22:35.840
to order of delta t squared,
and that's a really big change.
00:22:35.840 --> 00:22:40.280
So now if I cut the step
size by a factor of 10,
00:22:40.280 --> 00:22:43.904
I gain a factor of
100 in the accuracy.
00:22:43.904 --> 00:22:45.820
If I want to get six
more significant figures,
00:22:45.820 --> 00:22:48.910
I do 1,000 times more work,
not a million times more work.
00:22:48.910 --> 00:22:52.220
So it's a pretty
significant improvement.
00:22:52.220 --> 00:22:54.640
And then people have pushed
this on further and further.
00:22:54.640 --> 00:22:56.320
And so actually very
common integrators
00:22:56.320 --> 00:22:58.270
that you might use
nowadays would go out
00:22:58.270 --> 00:23:02.160
to fourth and fifth
order methods.
00:23:02.160 --> 00:23:04.242
And so they have
complicated update formulas
00:23:04.242 --> 00:23:05.950
that are carefully
designed to cancel out
00:23:05.950 --> 00:23:09.280
all the low order delta t errors
and just leave very high order
00:23:09.280 --> 00:23:10.390
ones.
00:23:10.390 --> 00:23:11.906
And that's kind of the--
00:23:11.906 --> 00:23:13.030
typical ones are like that.
00:23:13.030 --> 00:23:15.446
And then some fancy ones might
go even up to eighth order,
00:23:15.446 --> 00:23:18.337
I've seen, if you're
really pushing it.
00:23:18.337 --> 00:23:20.170
You have really complicated
g formulas then.
00:23:27.260 --> 00:23:31.980
Now, this is all
for the integration,
00:23:31.980 --> 00:23:35.950
but really our problem was
that the ODEs [INAUDIBLE]..
00:23:35.950 --> 00:23:41.780
So we have to do
f of y not f of t.
00:23:41.780 --> 00:23:45.070
Let's see if I
can show you that.
00:23:45.070 --> 00:23:52.140
Yeah, so we did the rectangle
rule, we had f of t.
00:23:52.140 --> 00:23:55.320
But the real problem is
we have to do f of y.
00:23:55.320 --> 00:23:57.690
Now the problem
is, we don't know
00:23:57.690 --> 00:24:00.710
what y is, y is our unknown.
00:24:00.710 --> 00:24:03.880
So we generally
have a problem here.
00:24:03.880 --> 00:24:06.530
In Forward Euler, we can get
away with it because we just
00:24:06.530 --> 00:24:08.930
put the y old value in.
00:24:08.930 --> 00:24:14.450
Like for example, if we
wanted to do midpoint here,
00:24:14.450 --> 00:24:17.690
we'd have to know f at
a different time point
00:24:17.690 --> 00:24:20.260
that we haven't computed yet,
because it's forward in time.
00:24:20.260 --> 00:24:21.635
Same thing with
trapezoidal rule.
00:24:21.635 --> 00:24:24.200
We have to know f at
some future timepoint.
00:24:24.200 --> 00:24:26.540
So it's not so easy,
actually, to go directly
00:24:26.540 --> 00:24:30.860
from these formulas to some
explicit formula like this
00:24:30.860 --> 00:24:33.810
with a g in here, where
the things inside g,
00:24:33.810 --> 00:24:37.810
the y's inside g, are all
y-values we actually know.
00:24:37.810 --> 00:24:39.729
So this is a serious problem.
00:24:45.600 --> 00:24:48.570
But this didn't stop
people from doing it.
00:24:48.570 --> 00:24:54.240
So one idea is you could go
back to the Taylor expansion,
00:24:54.240 --> 00:25:01.450
and now do the Taylor expansion
but in terms of f of y.
00:25:01.450 --> 00:25:03.700
So this is f of y of t.
00:25:06.240 --> 00:25:12.140
And then over here, this
is d squared y dt squared,
00:25:12.140 --> 00:25:20.470
but dy dt is f, so
this is really df dt.
00:25:20.470 --> 00:25:22.970
But f is actually not a function
of t, it's a function of y,
00:25:22.970 --> 00:25:23.540
right?
00:25:23.540 --> 00:25:24.980
But we can do chain rule.
00:25:24.980 --> 00:25:41.440
So df dt is equal to
the sum df dyn dyn dt.
00:25:45.630 --> 00:25:48.680
But dyn dt is f.
00:25:48.680 --> 00:25:50.960
And this is what we call
the Jacobian [INAUDIBLE]..
00:25:50.960 --> 00:25:56.558
So this is really saying that
this is equal to J times f.
00:25:59.550 --> 00:26:00.930
So that would be
one possibility,
00:26:00.930 --> 00:26:02.846
is you could go and write
the Taylor expansion
00:26:02.846 --> 00:26:05.190
out to the next higher
order explicitly,
00:26:05.190 --> 00:26:08.580
putting in J times f
here, and then just
00:26:08.580 --> 00:26:11.920
being in the cubic terms.
00:26:11.920 --> 00:26:13.370
Now, people don't
do this usually.
00:26:13.370 --> 00:26:15.720
Any idea why this is not
a popular thing to do?
00:26:18.996 --> 00:26:19.496
Yeah?
00:26:19.496 --> 00:26:22.059
AUDIENCE: It's computationally
expensive to [INAUDIBLE]..
00:26:22.059 --> 00:26:23.350
WILLIAM GREEN JR: That's right.
00:26:23.350 --> 00:26:23.891
That's right.
00:26:23.891 --> 00:26:25.570
So how many function
evaluations does it
00:26:25.570 --> 00:26:26.730
take to get the Jacobian?
00:26:32.270 --> 00:26:37.390
Right, so each function
is n values, right,
00:26:37.390 --> 00:26:38.930
because it's a vector.
00:26:38.930 --> 00:26:43.530
And you have to do at least,
say, forward differences
00:26:43.530 --> 00:26:46.700
would be n function evaluations.
00:26:46.700 --> 00:26:48.700
So instead of just doing
one function evaluation
00:26:48.700 --> 00:26:50.710
like we were doing with
the Forward Euler, now
00:26:50.710 --> 00:26:52.668
we have to do n plus 1,
or something like that,
00:26:52.668 --> 00:26:54.190
function evaluations.
00:26:54.190 --> 00:26:56.606
And actually, you probably
want to use center differences,
00:26:56.606 --> 00:26:59.260
so you need 2 times n plus 1.
00:26:59.260 --> 00:27:02.421
And if you do them analytically,
unless it's sparse,
00:27:02.421 --> 00:27:04.670
it's still going to be really
expensive to compute it.
00:27:04.670 --> 00:27:06.667
So it's a lot more
effort to do it.
00:27:06.667 --> 00:27:08.750
So you need to think, well,
is it really worth it?
00:27:08.750 --> 00:27:11.410
And so people have found
other formulas that basically
00:27:11.410 --> 00:27:13.510
get rid of the delta
t squared term that
00:27:13.510 --> 00:27:15.710
don't require so much
effort, and so that's
00:27:15.710 --> 00:27:17.224
what people do instead.
00:27:17.224 --> 00:27:18.140
But this is an option.
00:27:18.140 --> 00:27:18.806
You could do it.
00:27:25.580 --> 00:27:27.620
All right, so
instead, what people
00:27:27.620 --> 00:27:32.810
do is they've decided to
generalize the midpoint rule
00:27:32.810 --> 00:27:34.790
and try to figure
out, how can we
00:27:34.790 --> 00:27:39.380
get an estimate of the midpoint
value that's cheap to evaluate?
00:27:39.380 --> 00:27:40.970
And so what they
say is, well, this
00:27:40.970 --> 00:27:44.356
is g midpoint for
numerical integration.
00:27:44.356 --> 00:27:45.980
Let's make an new g
midpoint that works
00:27:45.980 --> 00:27:48.480
for when f is a function of y.
00:27:48.480 --> 00:27:52.700
So let's make it
function of y evaluated
00:27:52.700 --> 00:27:53.930
at t plus delta t over 2.
00:27:59.787 --> 00:28:01.370
So that's what the
formula looks like.
00:28:01.370 --> 00:28:02.870
But then you say,
well, I don't know
00:28:02.870 --> 00:28:05.510
y is at t plus delta t
over 2, because I only
00:28:05.510 --> 00:28:07.360
know it at y of t.
00:28:07.360 --> 00:28:10.140
So then they say, well, let's
do Forward Euler for that.
00:28:10.140 --> 00:28:14.720
So then they say, well, let's
say it's approximately f of--
00:28:18.580 --> 00:28:20.380
what's the Forward
Euler formula?
00:28:20.380 --> 00:28:28.442
It's y of t plus delta t
over 2 times f of y of t.
00:28:31.154 --> 00:28:33.018
How many parenthesis there?
00:28:38.410 --> 00:28:43.680
So this is the formula
that's actually
00:28:43.680 --> 00:28:45.360
used in practice a lot.
00:28:45.360 --> 00:28:48.120
And this one is called
Runge-Kutta two, second order
00:28:48.120 --> 00:28:50.030
Runge-Kutta formula.
00:28:50.030 --> 00:28:52.080
And it's just a--
00:28:52.080 --> 00:28:56.400
it's the midpoint rule where
we estimate the value of y of t
00:28:56.400 --> 00:28:59.074
plus delta t using
Forward Euler.
00:28:59.074 --> 00:29:01.240
And the advantage of this
is it's very cheap, right?
00:29:01.240 --> 00:29:05.152
Just evaluate-- I just have
to evaluate one function here,
00:29:05.152 --> 00:29:06.860
one function here, so
two function calls.
00:29:10.269 --> 00:29:16.860
[INAUDIBLE] and that way I
can get an update formula.
00:29:16.860 --> 00:29:19.810
And this one, turns out
that it has the nice delta t
00:29:19.810 --> 00:29:21.360
scaling that you might like.
00:29:24.640 --> 00:29:27.950
So that's the kind
of thing people do.
00:29:27.950 --> 00:29:30.600
And the Runge-Kutta
guys went crazy,
00:29:30.600 --> 00:29:33.260
and so one of the
most popular solvers
00:29:33.260 --> 00:29:36.920
is called Runge-Kutta four five,
and that's the ODE four five
00:29:36.920 --> 00:29:38.690
program in Matlab.
00:29:38.690 --> 00:29:41.257
And what that does is the
fourth order extension of this,
00:29:41.257 --> 00:29:42.840
and the fifth order
extension of this.
00:29:42.840 --> 00:29:45.050
And it uses them both,
and it compares them,
00:29:45.050 --> 00:29:47.900
and uses the difference
between an estimate the error
00:29:47.900 --> 00:29:48.860
in the calculation.
00:29:48.860 --> 00:29:52.434
And then uses that to decide
what size delta t you need,
00:29:52.434 --> 00:29:54.100
depending on what
tolerances you demand.
00:29:58.016 --> 00:29:59.352
And the formulas for all those--
00:29:59.352 --> 00:30:01.310
[INAUDIBLE] formulas are
given in the textbook.
00:30:06.810 --> 00:30:09.960
Now, this starts to
lead the problem.
00:30:09.960 --> 00:30:14.230
So we weren't able to actually
evaluate the true g midpoint,
00:30:14.230 --> 00:30:15.400
which we'd like to be here.
00:30:15.400 --> 00:30:18.690
Instead, we have to use
some approximation in order
00:30:18.690 --> 00:30:21.480
to extrapolate
forwards to the y.
00:30:21.480 --> 00:30:24.230
So this is hinting
at the whole problem.
00:30:24.230 --> 00:30:26.610
This whole problem is actually
we're just extrapolating.
00:30:26.610 --> 00:30:31.260
We know the true value
of y only at t naught,
00:30:31.260 --> 00:30:34.950
and all the rest we're doing
is extrapolating forwards.
00:30:34.950 --> 00:30:36.450
And we're do it
over and over again,
00:30:36.450 --> 00:30:37.949
so we're going to
do it for n steps.
00:30:37.949 --> 00:30:40.570
Maybe we'll do 1,000 steps.
00:30:40.570 --> 00:30:41.769
So you're extrapolating.
00:30:41.769 --> 00:30:43.310
And then based on
that extrapolation,
00:30:43.310 --> 00:30:44.760
you going to extrapolate again.
00:30:44.760 --> 00:30:47.610
And then based on extrapolation,
you going to extrapolate again.
00:30:47.610 --> 00:30:50.520
And you can see that this is
not really the most robust thing
00:30:50.520 --> 00:30:52.050
to do, right?
00:30:52.050 --> 00:30:54.090
Because if you have
any experience,
00:30:54.090 --> 00:30:57.480
you're not that comfortable
extrapolating at all.
00:30:57.480 --> 00:31:00.350
And then you have to
extrapolate 1,000 times.
00:31:00.350 --> 00:31:03.220
You should be a little bit
worried about what can happen.
00:31:03.220 --> 00:31:05.500
So let's think about how
the error will propagate.
00:31:05.500 --> 00:31:09.180
So we're going to have
some errors while we
00:31:09.180 --> 00:31:10.421
do these extrapolations.
00:31:18.780 --> 00:31:23.520
So we have here a t naught,
t naught, we're happy.
00:31:23.520 --> 00:31:24.810
We know y naught.
00:31:24.810 --> 00:31:26.110
This is exact.
00:31:26.110 --> 00:31:27.570
We're like woohoo.
00:31:27.570 --> 00:31:29.940
We know one value
of y, really good.
00:31:29.940 --> 00:31:34.716
All right, so now the
question is the first step.
00:31:34.716 --> 00:31:36.360
Now we're going
to get to value y1
00:31:36.360 --> 00:31:38.510
that we compute with
one of the formulas,
00:31:38.510 --> 00:31:40.260
and it's going to be
good to some accuracy
00:31:40.260 --> 00:31:43.470
depending on how good or
g is, our update formula.
00:31:43.470 --> 00:31:44.930
And so it's going
to have an error.
00:31:47.750 --> 00:31:50.290
I should comment
that this y naught,
00:31:50.290 --> 00:31:52.980
although we treat it as exact,
it could have an error too.
00:31:52.980 --> 00:31:55.188
But it will probably be more
like a machine precision
00:31:55.188 --> 00:31:57.010
kind of error because
we don't really
00:31:57.010 --> 00:31:58.080
know the initial
conditions that well.
00:31:58.080 --> 00:31:59.290
Or might be an
experimental error
00:31:59.290 --> 00:32:00.940
about how well we know what
the initial conditions are.
00:32:00.940 --> 00:32:02.590
But I'm not going to worry
about that one for now.
00:32:02.590 --> 00:32:04.600
But be aware, this is
true, we don't really
00:32:04.600 --> 00:32:07.700
know initial conditions
perfectly either.
00:32:07.700 --> 00:32:11.075
But certainly, even if we treat
this as being known perfectly,
00:32:11.075 --> 00:32:13.450
we're going to have some error
by the time we extrapolate
00:32:13.450 --> 00:32:15.660
to y1.
00:32:15.660 --> 00:32:17.660
So now what happens at y2?
00:32:17.660 --> 00:32:21.780
So now we've added some
more, so now we're at t2.
00:32:21.780 --> 00:32:28.350
And we see started from
y1 true plus some error.
00:32:28.350 --> 00:32:29.760
This is y-- y naught was true.
00:32:32.600 --> 00:32:35.520
Now we're starting from y1
true with some error in it,
00:32:35.520 --> 00:32:38.290
because of a mistake
we made the first step.
00:32:38.290 --> 00:32:42.240
And now we see how that
step propagates further.
00:32:42.240 --> 00:32:45.170
So let's see what
that's going to say.
00:32:45.170 --> 00:32:53.400
Well say y2 is going to
equal y1 plus delta t times
00:32:53.400 --> 00:32:55.710
some update formula
which would depend on y1.
00:33:00.770 --> 00:33:03.820
But y1 was not true
anymore, so this is actually
00:33:03.820 --> 00:33:05.990
equal to what should have
been the true value of y1
00:33:05.990 --> 00:33:09.176
if we only knew what it
was, but we don't know.
00:33:09.176 --> 00:33:11.480
And we have our error that
we don't know exactly how
00:33:11.480 --> 00:33:14.600
big it is that we added on.
00:33:14.600 --> 00:33:18.980
And then we have another
term from this guy,
00:33:18.980 --> 00:33:22.370
so it's going to
be delta t times
00:33:22.370 --> 00:33:28.520
g of y1 true plus delta 1.
00:33:31.240 --> 00:33:34.240
And then we know our update
formula's not right anyway,
00:33:34.240 --> 00:33:37.621
so there's actually
another error, delta 2,
00:33:37.621 --> 00:33:39.621
just because of the error
in the update formula.
00:33:45.030 --> 00:33:48.089
Now if we had a
great update formula,
00:33:48.089 --> 00:33:49.380
these deltas are kind of small.
00:33:52.240 --> 00:33:57.844
And if somehow, by
luck, our delta 1 was 0,
00:33:57.844 --> 00:33:59.260
and we had a great
update formula,
00:33:59.260 --> 00:34:01.926
then we think that this error in
y2 is going to be pretty small.
00:34:04.440 --> 00:34:05.580
But it's not true.
00:34:05.580 --> 00:34:07.080
So really the errors
are pretty big.
00:34:10.880 --> 00:34:13.790
So let's see if we
can figure out--
00:34:13.790 --> 00:34:17.010
we're going to do a
Taylor expansion here,
00:34:17.010 --> 00:34:27.139
so this'll be g of
y1 true plus d dgy--
00:34:27.139 --> 00:34:29.800
the Jacobian of g
with respect to y--
00:34:29.800 --> 00:34:33.199
times the delta 1.
00:34:33.199 --> 00:34:36.884
That's the first order
Taylor extension.
00:34:36.884 --> 00:34:38.550
All right, so now
let's group the terms.
00:34:48.679 --> 00:34:58.580
So y2 is equal to y1 true
plus the update formula of y1
00:34:58.580 --> 00:35:03.710
true plus the error
in our update formula
00:35:03.710 --> 00:35:08.310
if we had started
at the true value.
00:35:08.310 --> 00:35:09.780
So this is what we--
00:35:09.780 --> 00:35:12.320
if somebody told us what
the true value of y1, this
00:35:12.320 --> 00:35:14.020
is the size error we'd expect.
00:35:14.020 --> 00:35:17.390
It's to be that little
error we calculated before.
00:35:17.390 --> 00:35:20.600
But then on top of it, we
have these other terms.
00:35:20.600 --> 00:35:25.490
We have a delta 1, because
y1 was not the true value.
00:35:25.490 --> 00:35:26.780
Sorry, I lost a delta t here.
00:35:31.490 --> 00:35:46.080
And then we have delta t
times dg dy times delta 1.
00:35:46.080 --> 00:35:47.860
And I'll emphasize
these are all vectors.
00:35:53.283 --> 00:35:56.670
I think that's it if we
only keep to first order.
00:36:00.140 --> 00:36:01.970
So what is this thing?
00:36:01.970 --> 00:36:06.833
This is a matrix, right, dg dy.
00:36:06.833 --> 00:36:09.600
I could write it this way.
00:36:09.600 --> 00:36:13.980
It's the identity matrix plus
delta t times the Jacobian dg
00:36:13.980 --> 00:36:18.500
dy times the vector d1.
00:36:21.990 --> 00:36:24.430
And you can see, at every
iteration when I do this,
00:36:24.430 --> 00:36:27.500
I'm going to get a
similar factor like this.
00:36:27.500 --> 00:36:30.820
And it's going to
again multiply delta 1.
00:36:30.820 --> 00:36:32.385
I keep on multiplying the--
00:36:32.385 --> 00:36:33.760
errors from my
previous step keep
00:36:33.760 --> 00:36:37.090
getting multiplied by
this kind of factor.
00:36:37.090 --> 00:36:38.840
So it's going to be
really important to us
00:36:38.840 --> 00:36:40.540
about what the norm
of this matrix is.
00:36:52.641 --> 00:36:55.099
What's going to happen if the
normal of this matrix is big?
00:36:58.452 --> 00:36:59.410
What's going to happen?
00:36:59.410 --> 00:37:00.534
AUDIENCE: [INAUDIBLE]
00:37:00.534 --> 00:37:02.575
WILLIAM GREEN JR: That's
right, the error's going
00:37:02.575 --> 00:37:04.210
to get multiplied
by some factor,
00:37:04.210 --> 00:37:06.452
say bigger than 1,
every iteration.
00:37:06.452 --> 00:37:08.160
And after we go 1,000
iterations, how big
00:37:08.160 --> 00:37:11.770
is the error going to be?
00:37:11.770 --> 00:37:12.640
Really big, right?
00:37:12.640 --> 00:37:16.060
Even if that thing is quite
close to 1, 1 plus something
00:37:16.060 --> 00:37:18.760
to the 1,000 power
it's gigantic.
00:37:18.760 --> 00:37:21.010
So what was originally
a pretty small error--
00:37:21.010 --> 00:37:22.880
because we chose
a good g method--
00:37:22.880 --> 00:37:26.020
is going to get multiplied by
this huge amplification factor.
00:37:26.020 --> 00:37:27.850
So the key thing
for a method to be
00:37:27.850 --> 00:37:30.190
what is called
numerically stable
00:37:30.190 --> 00:37:35.280
is that the norm of this matrix
has got to be less than 1.
00:37:35.280 --> 00:37:37.700
Or maybe even equal,
equals 1 might be OK.
00:37:37.700 --> 00:37:40.860
If this is much bigger
than 1, you're doomed.
00:37:40.860 --> 00:37:43.390
Your numeric order is
just going to blow up
00:37:43.390 --> 00:37:46.094
from step to step to step.
00:37:46.094 --> 00:37:48.010
And what's bad about it,
is actually sometimes
00:37:48.010 --> 00:37:50.361
you might not even
be aware of this.
00:37:50.361 --> 00:37:52.860
So you'll get some solution,
because it still will calculate
00:37:52.860 --> 00:37:57.070
some numbers for the y's at
each type-- y1, y2, y3, y4--
00:37:57.070 --> 00:38:00.740
it's just running
through the calculation.
00:38:00.740 --> 00:38:04.070
But the numbers may become
increasingly meaningless
00:38:04.070 --> 00:38:07.970
further away from the y true as
time goes on, as the steps go.
00:38:14.275 --> 00:38:20.100
AUDIENCE: [INAUDIBLE]
00:38:20.100 --> 00:38:21.469
WILLIAM GREEN JR: Yes, it can.
00:38:21.469 --> 00:38:24.010
Yeah, so that's what we call a
stable numerical method is one
00:38:24.010 --> 00:38:26.440
that, if you make an
error in an earlier step,
00:38:26.440 --> 00:38:29.980
its impact diminishes
as the steps go on.
00:38:29.980 --> 00:38:31.330
That's what you want.
00:38:31.330 --> 00:38:32.580
You can't always achieve this.
00:38:32.580 --> 00:38:35.490
But if you can,
that's really great.
00:38:35.490 --> 00:38:37.930
Then that's you would call
a really stable, robust kind
00:38:37.930 --> 00:38:39.080
of method.
00:38:39.080 --> 00:38:42.070
Even though I give you
some stupidity early on,
00:38:42.070 --> 00:38:43.570
my stupidity
evaporates away, and it
00:38:43.570 --> 00:38:45.455
goes to the true solution.
00:38:45.455 --> 00:38:46.330
That's what you want.
00:38:48.847 --> 00:38:50.680
OK, so let's think of
cases where it's going
00:38:50.680 --> 00:38:52.660
to be problematic for sure.
00:38:52.660 --> 00:38:56.680
So if dg dy--
00:38:56.680 --> 00:39:00.860
say if all it's
eigenvalues are positive
00:39:00.860 --> 00:39:03.540
and you add it to
identity matrix,
00:39:03.540 --> 00:39:06.890
the thing is still going to
have eigenvalues bigger than 1.
00:39:06.890 --> 00:39:09.350
Can you see that?
00:39:09.350 --> 00:39:12.357
And actually, what do we
say when dg dy is positive,
00:39:12.357 --> 00:39:13.440
has a positive eigenvalue?
00:39:13.440 --> 00:39:14.600
You guys have a word
for this already
00:39:14.600 --> 00:39:16.558
you learn in your
differential equations class.
00:39:16.558 --> 00:39:20.030
What do you call these kind
of differential equations?
00:39:20.030 --> 00:39:20.870
Maybe not dg dy.
00:39:20.870 --> 00:39:22.940
Back up, how about df dy?
00:39:22.940 --> 00:39:26.270
If df dy, if that Jacobian
has a positive eigenvalue,
00:39:26.270 --> 00:39:27.580
what do you say?
00:39:34.230 --> 00:39:35.710
You guys remember this?
00:39:35.710 --> 00:39:39.700
So this is what we call
unstable differential equations.
00:39:39.700 --> 00:39:42.970
So if it has any
positive eigenvalues,
00:39:42.970 --> 00:39:49.900
then two initial conditions that
differ by a little differential
00:39:49.900 --> 00:39:53.650
exponentially separate
as time goes on as long
00:39:53.650 --> 00:39:56.050
as the difference as any
projection in the direction
00:39:56.050 --> 00:39:59.630
of the bad eigenvector.
00:39:59.630 --> 00:40:01.970
So those are called
numerically unstable.
00:40:01.970 --> 00:40:04.980
A lot of those actually
exist in reality.
00:40:04.980 --> 00:40:08.030
So explosions,
that's because you
00:40:08.030 --> 00:40:12.200
had some positive eigenvalue
somewhere for example.
00:40:12.200 --> 00:40:14.090
Amplifiers, like
you buy an amplifier
00:40:14.090 --> 00:40:17.900
to crank up the music, you
get positive feedback when--
00:40:17.900 --> 00:40:20.360
remember Professor Swan
walked to some strange place,
00:40:20.360 --> 00:40:22.967
and for some reason, he got
the feedback from the speakers?
00:40:22.967 --> 00:40:24.800
So he had a positive
eigenvalue going there.
00:40:24.800 --> 00:40:26.716
A little tiny noise and
his microphone, and it
00:40:26.716 --> 00:40:28.700
amplified to a big squeak.
00:40:28.700 --> 00:40:31.790
So there are real systems
like this that amplify.
00:40:31.790 --> 00:40:33.332
And sometimes we
want to model them,
00:40:33.332 --> 00:40:34.540
like explosions for examples.
00:40:34.540 --> 00:40:37.820
That's a pretty important
for chemical engineers.
00:40:37.820 --> 00:40:41.810
But it's going to be really
tough, because if it's
00:40:41.810 --> 00:40:44.990
got positive eigenvalues,
then this is probably going
00:40:44.990 --> 00:40:46.850
to have a norm bigger than 1.
00:40:46.850 --> 00:40:48.320
And then whatever
little errors we
00:40:48.320 --> 00:40:50.070
may get any step in
the whole calculation,
00:40:50.070 --> 00:40:52.130
we're going to amplify
and amplify and amplify.
00:40:52.130 --> 00:40:53.630
And who knows if
we're going to have
00:40:53.630 --> 00:40:57.110
any reality left in our solution
by the time we get to the end.
00:41:02.642 --> 00:41:04.350
All right, so that's
one kind of problem.
00:41:09.520 --> 00:41:13.980
So how about we have a problem
that's intrinsically stable.
00:41:13.980 --> 00:41:16.140
That means that
all the eigenvalues
00:41:16.140 --> 00:41:19.320
of the Jacobian of the
original problem are negative.
00:41:19.320 --> 00:41:22.770
So the thing should
be perfectly stable.
00:41:22.770 --> 00:41:25.254
From whatever initial
condition they started from,
00:41:25.254 --> 00:41:26.920
the whole thing should
sort of come down
00:41:26.920 --> 00:41:28.300
to like an equilibrium point.
00:41:28.300 --> 00:41:29.340
So that's a pretty
common situation
00:41:29.340 --> 00:41:30.506
in chemical engineering too.
00:41:30.506 --> 00:41:34.410
You start from some
state, and it relaxes down
00:41:34.410 --> 00:41:37.040
to like an equilibrium
state or a steady state.
00:41:37.040 --> 00:41:39.490
Well-behaved systems
act like that.
00:41:39.490 --> 00:41:44.180
So what I'm going
to show you is that,
00:41:44.180 --> 00:41:48.470
despite the fact that the
df dy is well-behaved,
00:41:48.470 --> 00:41:52.460
it's possible that this
sum matrix might still
00:41:52.460 --> 00:41:56.960
be poorly behaved and still
have a norm bigger that one.
00:41:56.960 --> 00:41:58.960
And so that means you
started for a problem that
00:41:58.960 --> 00:42:01.920
was pretty well-behaved, and
you broke it by your choice
00:42:01.920 --> 00:42:04.030
in numerical method.
00:42:04.030 --> 00:42:06.550
So this is a very bad one.
00:42:06.550 --> 00:42:09.370
This is the kind that can
make you really embarrassed
00:42:09.370 --> 00:42:11.440
if you're working in a company.
00:42:11.440 --> 00:42:13.000
Your boss gives you a problem.
00:42:13.000 --> 00:42:16.060
He says, please tell me the time
is going to take our reactor
00:42:16.060 --> 00:42:18.070
to relax after a start up.
00:42:18.070 --> 00:42:19.774
the system is
perfectly well-behaved.
00:42:19.774 --> 00:42:21.190
He gives these
beautiful equations
00:42:21.190 --> 00:42:23.121
that your previous
coworker worked out.
00:42:23.121 --> 00:42:23.870
Everything's fine.
00:42:23.870 --> 00:42:25.000
You put it in there.
00:42:25.000 --> 00:42:26.800
You use your stupid
numerical solver,
00:42:26.800 --> 00:42:30.700
and you get a wild oscillation,
and the thing blows up.
00:42:30.700 --> 00:42:33.610
OK, and that's going to be
because of a poor choice of g
00:42:33.610 --> 00:42:39.700
so that it makes this norm of
this matrix bigger than one.
00:42:39.700 --> 00:42:41.360
It's pretty easy
for that to happen.
00:42:41.360 --> 00:42:44.792
So you don't need a very
complicated case to show it.
00:42:52.930 --> 00:42:54.535
So let's just do a scalar case.
00:42:58.150 --> 00:43:00.210
Say dy dt is equal
to negative 2y.
00:43:02.750 --> 00:43:05.370
So this is a well-behaved
differential equation.
00:43:05.370 --> 00:43:08.120
What's the solution?
00:43:08.120 --> 00:43:14.335
So with a y naught y at
t equals 0 is equal to 1.
00:43:17.800 --> 00:43:19.780
What's the solution?
00:43:19.780 --> 00:43:21.770
AUDIENCE: [INAUDIBLE]
00:43:21.770 --> 00:43:25.790
WILLIAM GREEN JR:
Yeah, OK, so it's
00:43:25.790 --> 00:43:27.460
perfectly well-behaved function.
00:43:27.460 --> 00:43:31.460
And the plot when you
plot it, here's one.
00:43:31.460 --> 00:43:33.470
It goes [PLANE DIVING NOISE].
00:43:36.800 --> 00:43:39.890
So now let's solve this with
the Forward Euler method.
00:43:39.890 --> 00:43:46.690
So I'm going to say
that y new is equal to y
00:43:46.690 --> 00:43:48.946
plus delta t times f.
00:43:52.570 --> 00:43:58.360
So in this case, this is f,
right there, negative 2y.
00:43:58.360 --> 00:44:00.200
So for Forward
Euler, this is just
00:44:00.200 --> 00:44:01.926
saying that this is negative 2y.
00:44:05.350 --> 00:44:08.700
All right, so let's
choose a delta t,
00:44:08.700 --> 00:44:12.350
say delta t equals
1 to start with.
00:44:12.350 --> 00:44:16.480
So we're going to compute.
00:44:16.480 --> 00:44:17.540
We start at t equals 0.
00:44:17.540 --> 00:44:21.540
We're computing out to t
equals 1. so we're out here.
00:44:21.540 --> 00:44:24.400
So we put 1 in here.
00:44:24.400 --> 00:44:28.590
Our initial condition
was y was 1, so that's 1.
00:44:28.590 --> 00:44:33.840
So as 1 times negative
2 times 1 is negative 2.
00:44:33.840 --> 00:44:35.410
This is a 1.
00:44:35.410 --> 00:44:40.770
1 minus 2 is negative 1.
00:44:40.770 --> 00:44:44.898
So at the first step, it
goes from here down to there.
00:44:49.379 --> 00:44:50.420
This is not looking good.
00:44:54.200 --> 00:44:55.750
This physical
variable is supposed
00:44:55.750 --> 00:44:57.770
to be gradually
relaxing towards 0.
00:44:57.770 --> 00:45:00.580
Going negative is
not what you want.
00:45:00.580 --> 00:45:02.990
And then I think we could go
the next step, if you want.
00:45:02.990 --> 00:45:07.432
So we can put in
suppose y is negative 1.
00:45:07.432 --> 00:45:09.640
And we're still doing
a delta t of one step.
00:45:09.640 --> 00:45:14.610
So negative 1 times
negative 2 is positive 2.
00:45:14.610 --> 00:45:18.731
Plus negative 1, so this goes
back to the starting point.
00:45:18.731 --> 00:45:24.484
We'll get the sawtooth,
is what you get.
00:45:24.484 --> 00:45:26.900
And it turns out, actually,
if you make delta t any larger
00:45:26.900 --> 00:45:30.630
than 1, then this
sawtooth amplifies,
00:45:30.630 --> 00:45:34.320
so you have a oscillating
solution that explodes.
00:45:34.320 --> 00:45:37.170
Which is a little bit different
than the physical solution,
00:45:37.170 --> 00:45:38.100
which gently relaxes.
00:45:42.230 --> 00:45:46.220
And so this is just
a scalar equation
00:45:46.220 --> 00:45:48.260
that you guys could have
done in high school.
00:45:48.260 --> 00:45:49.710
If you high school math
teacher had been brave,
00:45:49.710 --> 00:45:51.160
he would have showed it you.
00:45:51.160 --> 00:45:54.682
But he don't want to ruin your
faith in the Forward Euler
00:45:54.682 --> 00:45:57.140
method, so he just showed you
a simple one that worked out.
00:46:00.320 --> 00:46:01.990
Now you can make
this work better
00:46:01.990 --> 00:46:04.650
by making the delta t smaller.
00:46:04.650 --> 00:46:07.950
If you made the delta
t small enough--
00:46:07.950 --> 00:46:11.550
in this case, I think it's
delta t is less than half--
00:46:11.550 --> 00:46:15.250
then this actually behaves
somewhat reasonably.
00:46:15.250 --> 00:46:21.420
And so what's going on there
is that the identity matrix
00:46:21.420 --> 00:46:23.660
is positive.
00:46:23.660 --> 00:46:27.880
My Jacobian is negative.
00:46:27.880 --> 00:46:30.260
And if I make delta
t small enough,
00:46:30.260 --> 00:46:32.755
this big positive number
is bigger than this smaller
00:46:32.755 --> 00:46:34.880
negative number, and the
whole thing stays positive
00:46:34.880 --> 00:46:38.910
but less than 1, and so
then the normal is good.
00:46:38.910 --> 00:46:41.110
But if I make delta
t too large, I
00:46:41.110 --> 00:46:43.615
make this second
term much larger.
00:46:43.615 --> 00:46:46.390
If it gets much larger than 1,
it overwhelms the first term.
00:46:46.390 --> 00:46:51.480
The thing's negative, but I
what do the norm it's positive.
00:46:51.480 --> 00:46:53.360
And so it has an expansion.
00:46:53.360 --> 00:46:55.790
Then the negative is why
we get these oscillations.
00:46:55.790 --> 00:46:58.290
So if you ever do a numerical
differential equation solution
00:46:58.290 --> 00:47:00.330
and you start seeing
this kind of stuff,
00:47:00.330 --> 00:47:06.090
it's almost certainly related
to this numerical instability
00:47:06.090 --> 00:47:07.590
of the method.
00:47:07.590 --> 00:47:09.420
and not really a physical thing.
00:47:09.420 --> 00:47:10.140
I mean, there are
physical things
00:47:10.140 --> 00:47:12.473
that do this, but not too
often in chemical engineering,
00:47:12.473 --> 00:47:14.355
fortunately.
00:47:14.355 --> 00:47:18.230
All right, questions?
00:47:18.230 --> 00:47:19.040
No questions, OK.
00:47:19.040 --> 00:47:20.956
So in the textbook, they
have a big discussion
00:47:20.956 --> 00:47:25.190
about how to evaluate
different g methods in order
00:47:25.190 --> 00:47:27.800
to figure out if they're
guaranteed to be stable,
00:47:27.800 --> 00:47:30.410
or under what conditions
they'll be stable.
00:47:30.410 --> 00:47:34.184
And so that's worth
a look to see.
00:47:34.184 --> 00:47:35.600
And they have
detailed derivations
00:47:35.600 --> 00:47:36.543
for several of them.
00:47:40.000 --> 00:47:42.444
Questions?
00:47:42.444 --> 00:47:44.360
All right, so for the
homework, homework four,
00:47:44.360 --> 00:47:47.870
you're going to have to
write your own ODE solver.
00:47:47.870 --> 00:47:52.390
So you'll do update
formulas and compare them.
00:47:52.390 --> 00:47:56.180
When we come back
on Tuesday, I'll
00:47:56.180 --> 00:47:57.830
show you how this
is usually dealt
00:47:57.830 --> 00:48:04.530
with by switching to implicit
solvers and talk about that.
00:48:04.530 --> 00:48:06.980
All right, have a good weekend.
00:48:06.980 --> 00:48:09.470
Enjoy the three days off.