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:18.980
at ocw.mit.edu.
00:00:23.865 --> 00:00:25.740
PROFESSOR: All right,
I think we're gathered.
00:00:25.740 --> 00:00:28.270
So let's get going.
00:00:28.270 --> 00:00:30.510
We were going to
DAEs today, but we
00:00:30.510 --> 00:00:33.812
decided to postpone for one
day and do a little catch up
00:00:33.812 --> 00:00:35.520
on different things
that you should know.
00:00:38.850 --> 00:00:43.110
And I thought I'd
first back up a minute
00:00:43.110 --> 00:00:47.050
and recall what
you've learned so far.
00:00:47.050 --> 00:00:50.970
So we started out talking
a lot about linear algebra
00:00:50.970 --> 00:00:54.300
and, particularly, the
problems of this type
00:00:54.300 --> 00:00:56.640
where you're given a matrix
and you're given a vector,
00:00:56.640 --> 00:00:58.306
and you want to figure
out what the x is
00:00:58.306 --> 00:01:00.330
that makes this equation true.
00:01:00.330 --> 00:01:02.010
And you've become
extremely familiar
00:01:02.010 --> 00:01:06.340
with the convenient
Matlab backslash function.
00:01:06.340 --> 00:01:08.330
You've probably used
this a lot of times.
00:01:08.330 --> 00:01:11.430
And if you remember, this is one
of the only problems where we
00:01:11.430 --> 00:01:14.310
often have a unique solution.
00:01:14.310 --> 00:01:16.320
We know a solution exists.
00:01:16.320 --> 00:01:18.630
And we can do a
finite number of steps
00:01:18.630 --> 00:01:20.310
in order to get the solution.
00:01:20.310 --> 00:01:22.870
And actually, the method
we have is really good.
00:01:22.870 --> 00:01:25.272
So it almost always
gives us the solution
00:01:25.272 --> 00:01:27.580
and to machine accuracy.
00:01:27.580 --> 00:01:30.030
So this is a brilliant method.
00:01:30.030 --> 00:01:32.640
And so everything is based
on this in our whole field,
00:01:32.640 --> 00:01:35.140
actually.
00:01:35.140 --> 00:01:36.750
So we learned this first.
00:01:36.750 --> 00:01:43.200
We learned it takes
n cubed operations
00:01:43.200 --> 00:01:46.012
in order to get
the solution, where
00:01:46.012 --> 00:01:47.220
n is the size of the vectors.
00:01:51.140 --> 00:01:59.470
And when we do the solution,
the x we get for backslash,
00:01:59.470 --> 00:02:01.690
invariably, it
either is a solution,
00:02:01.690 --> 00:02:03.537
solves this equation
to machine accuracy,
00:02:03.537 --> 00:02:05.245
in the sense that you
multiply m times x,
00:02:05.245 --> 00:02:08.900
you get b to as many
digits as you can read.
00:02:08.900 --> 00:02:12.220
However, we talked about how
m could be poorly conditioned
00:02:12.220 --> 00:02:13.870
so that there might
be multiple x's
00:02:13.870 --> 00:02:16.880
that would solve this
equation to machine accuracy.
00:02:16.880 --> 00:02:19.800
And so that's the problem
with ill conditioning.
00:02:19.800 --> 00:02:21.770
So it's not that it
doesn't solve the equation.
00:02:21.770 --> 00:02:24.077
It's that we're not really
sure that's the right x.
00:02:24.077 --> 00:02:25.910
You could get a bunch
of different x's that,
00:02:25.910 --> 00:02:29.890
basically, solve it as
well as we can tell.
00:02:29.890 --> 00:02:31.986
All right, so do you
guys remember this?
00:02:31.986 --> 00:02:34.450
All right, and so
that also implies
00:02:34.450 --> 00:02:36.550
that if you make a
very small change in b,
00:02:36.550 --> 00:02:39.496
you might get a
gigantic change in x.
00:02:39.496 --> 00:02:41.995
Because even for one b, you
have a quite a range of x's that
00:02:41.995 --> 00:02:45.930
could work, all right.
00:02:45.930 --> 00:02:48.570
Then we try to solve
problems like this.
00:02:51.880 --> 00:02:55.470
Right, systems of equations,
not just linear ones,
00:02:55.470 --> 00:02:57.400
but nonlinear ones as well.
00:02:57.400 --> 00:03:00.890
And what we decided to do
as one of our best methods
00:03:00.890 --> 00:03:06.115
is Newton-Raphson,
where we keep on doing j
00:03:06.115 --> 00:03:10.680
backslash f
[INAUDIBLE] negative f.
00:03:10.680 --> 00:03:15.200
Right, so all we did in
this method, a lot of times,
00:03:15.200 --> 00:03:19.410
is we just kept on calling this
method over and over again.
00:03:19.410 --> 00:03:21.520
All right, so this
took order of n cubes.
00:03:21.520 --> 00:03:24.730
This thing is going
to take order of n,
00:03:24.730 --> 00:03:27.589
I don't know what you call
it, steps or something.
00:03:27.589 --> 00:03:29.130
Maybe if you guys
have a word for it.
00:03:29.130 --> 00:03:31.560
What do you call the iterations?
00:03:31.560 --> 00:03:33.770
Iterations.
00:03:33.770 --> 00:03:38.746
So it takes the number of
iterations times n cubed.
00:03:38.746 --> 00:03:42.110
Now, if you really can
make Newton-Raphson work,
00:03:42.110 --> 00:03:43.910
and you have a really
great initial guess,
00:03:43.910 --> 00:03:45.440
the number of iterations
is really tiny.
00:03:45.440 --> 00:03:47.606
Because Newton-Raphson is
such a fantastically great
00:03:47.606 --> 00:03:48.440
convergence.
00:03:48.440 --> 00:03:51.410
However, we're usually
not such good guessers.
00:03:51.410 --> 00:03:53.720
And so it might take
a lot of iterations.
00:03:53.720 --> 00:03:56.120
And in fact, we use
the dogleg method
00:03:56.120 --> 00:03:59.060
to solve this as we've
discussed about where you mix up
00:03:59.060 --> 00:03:59.840
different methods.
00:03:59.840 --> 00:04:02.445
And you have a trust region
and all this baloney.
00:04:02.445 --> 00:04:04.070
And so it might
actually be quite a lot
00:04:04.070 --> 00:04:05.194
of iterations to get there.
00:04:09.170 --> 00:04:12.140
And then we went and tried
to solve this problem.
00:04:19.459 --> 00:04:22.145
And what we decided to do
for both of these was almost
00:04:22.145 --> 00:04:23.550
do the same thing.
00:04:23.550 --> 00:04:28.810
So I said, well, we can
really solve this one
00:04:28.810 --> 00:04:30.190
by instead solving this.
00:04:41.510 --> 00:04:44.820
Right, that's really
what we ended up doing.
00:04:44.820 --> 00:04:47.310
That's what happens
inside fsolve is
00:04:47.310 --> 00:04:49.215
it converts the
problem of finding
00:04:49.215 --> 00:04:53.490
f of x equals 0 into
a minimization problem
00:04:53.490 --> 00:04:56.040
to try to make the norm
of f as small as possible.
00:04:56.040 --> 00:04:59.320
Because the answer will
be when the norm is zero.
00:04:59.320 --> 00:05:01.660
And if f of x equals
0 has a solution,
00:05:01.660 --> 00:05:04.690
then the global minimum
here is a solution,
00:05:04.690 --> 00:05:06.850
as long as your initial
guess is good enough
00:05:06.850 --> 00:05:11.320
that you're in the valley that
leads to the true minimum.
00:05:11.320 --> 00:05:12.590
Then this will work.
00:05:12.590 --> 00:05:14.620
And that's actually
what fsolve does.
00:05:14.620 --> 00:05:18.490
So these two guys get
kind of coupled together,
00:05:18.490 --> 00:05:23.805
the minimization problem and
the root-finding problem,
00:05:23.805 --> 00:05:25.180
or the problem of
solving systems
00:05:25.180 --> 00:05:26.830
of non-linear equations.
00:05:26.830 --> 00:05:29.740
And in both of them, we're
going to call the first method
00:05:29.740 --> 00:05:31.400
a million times.
00:05:31.400 --> 00:05:36.940
OK, so that's the way
we're working so far.
00:05:36.940 --> 00:05:39.220
And actually, when
you do fmincon,
00:05:39.220 --> 00:05:41.890
which was the next
level, then you actually
00:05:41.890 --> 00:05:43.390
end up looking
for saddle points,
00:05:43.390 --> 00:05:46.540
which, again, we look
for as minimizing
00:05:46.540 --> 00:05:51.340
the norm of a gradient
squared, which gets it back
00:05:51.340 --> 00:05:53.254
into this kind of form.
00:05:53.254 --> 00:05:54.670
And so they're all
the same thing.
00:05:57.470 --> 00:06:01.270
And you can, again,
work out estimates
00:06:01.270 --> 00:06:03.160
of how much effort it takes.
00:06:03.160 --> 00:06:06.430
And it's something times
n cubed, Typically,
00:06:06.430 --> 00:06:09.500
if most your work is this.
00:06:09.500 --> 00:06:12.880
Now, I should make a comment
that, in practice, that's
00:06:12.880 --> 00:06:14.060
not always true.
00:06:14.060 --> 00:06:17.320
So in practice, sometimes,
although the highest order term
00:06:17.320 --> 00:06:20.820
is n cubed, the linear
algebra programs that do this
00:06:20.820 --> 00:06:24.370
are so great that the
prefactor is tiny.
00:06:24.370 --> 00:06:26.910
So in order to get
the limit where
00:06:26.910 --> 00:06:29.280
the n cubed term was actually
dominate the CPU time,
00:06:29.280 --> 00:06:31.230
you had to have a gigantic n.
00:06:31.230 --> 00:06:34.170
And in fact, the limit,
it's often the effort
00:06:34.170 --> 00:06:35.770
to construct the matrix j.
00:06:35.770 --> 00:06:37.920
And so it really
goes as n squared.
00:06:37.920 --> 00:06:40.589
So although this nominally
is like n cubed, in practice,
00:06:40.589 --> 00:06:41.880
a lot of times, it's n squared.
00:06:44.610 --> 00:06:46.680
And I say it for
these other methods.
00:06:46.680 --> 00:06:48.540
And then a lot of
them get into how good
00:06:48.540 --> 00:06:49.620
the initial guesses are.
00:06:49.620 --> 00:06:53.467
Because the number of
iterations can get crazy.
00:06:53.467 --> 00:06:55.050
And in fact, maybe,
it never converges
00:06:55.050 --> 00:06:56.820
if your initial guess is bad.
00:06:56.820 --> 00:06:58.904
So that can be a
really critical thing.
00:06:58.904 --> 00:06:59.820
[INAUDIBLE] is better.
00:06:59.820 --> 00:07:01.380
Because you can't
change n I mean,
00:07:01.380 --> 00:07:04.469
how many unknowns you've got is
how many unknowns you've got.
00:07:04.469 --> 00:07:06.510
But you can change the
number iterations by, say,
00:07:06.510 --> 00:07:08.410
more clever initial guessing.
00:07:08.410 --> 00:07:11.550
And that was the whole concept
about continuation and homotopy
00:07:11.550 --> 00:07:13.110
and all this kind
of stuff was trying
00:07:13.110 --> 00:07:15.210
to help reduce this number.
00:07:18.440 --> 00:07:23.490
Last time, we talked about
implicitly solving ODEs.
00:07:23.490 --> 00:07:26.740
And so we had an update formula.
00:07:26.740 --> 00:07:28.320
So we were solving an ODE.
00:07:34.890 --> 00:07:37.122
The initial value
problem, y of--
00:07:41.010 --> 00:07:44.580
Right, and we said
that, first, we went
00:07:44.580 --> 00:07:46.440
through the explicit methods.
00:07:46.440 --> 00:07:48.687
And they were perfectly
straightforward to do.
00:07:48.687 --> 00:07:51.270
And you don't have to solve any
systems of nonlinear equations
00:07:51.270 --> 00:07:52.020
to do them.
00:07:52.020 --> 00:07:53.500
And so that was great.
00:07:53.500 --> 00:07:55.375
But then we pointed out
that, a lot of times,
00:07:55.375 --> 00:07:56.980
they're numerically unstable.
00:07:56.980 --> 00:08:01.050
And so then you have to
instead use implicit methods.
00:08:01.050 --> 00:08:03.960
But the implicit methods
generally work this way.
00:08:03.960 --> 00:08:05.210
They have your new value of y.
00:08:05.210 --> 00:08:09.660
What you're trying to compute
is your current value of y
00:08:09.660 --> 00:08:13.930
plus delta t sums g.
00:08:13.930 --> 00:08:18.305
It depends on the size
of delta t and y new.
00:08:18.305 --> 00:08:19.518
Maybe y old also.
00:08:24.000 --> 00:08:26.870
And these g's are things like
the Crank-Nicolson, which
00:08:26.870 --> 00:08:28.590
you're going to work
on in the homework,
00:08:28.590 --> 00:08:31.364
and we show it in class too.
00:08:31.364 --> 00:08:33.030
And so there's different
update formulas
00:08:33.030 --> 00:08:36.059
of people suggesting
how to do this.
00:08:36.059 --> 00:08:39.870
The trick of this is that the
y new is inside the function.
00:08:39.870 --> 00:08:45.470
So therefore, this is
really another problem
00:08:45.470 --> 00:08:48.160
like we had before.
00:08:48.160 --> 00:08:49.820
You can really
rewrite this that way
00:08:49.820 --> 00:08:51.860
by just putting this
on the one side.
00:08:51.860 --> 00:08:53.455
Put them all on
one side of zero.
00:08:53.455 --> 00:08:55.455
And now, I'm trying to
figure out what y new is.
00:08:55.455 --> 00:08:58.350
And I'm going to solve
that using this method.
00:08:58.350 --> 00:09:01.800
And I'm going to solve
that by using this method.
00:09:01.800 --> 00:09:05.960
All right, so now, how is
the scaling going to go?
00:09:05.960 --> 00:09:12.140
At each time step, I have
to solve one of these guys.
00:09:12.140 --> 00:09:15.199
How much effort did that take?
00:09:15.199 --> 00:09:15.990
We did that, right.
00:09:15.990 --> 00:09:19.020
Here it is, n [? after ?]
times n squared.
00:09:19.020 --> 00:09:21.050
And then how many times
do I have to do it?
00:09:23.894 --> 00:09:25.728
AUDIENCE: It's h times
7 times [INAUDIBLE]..
00:09:25.728 --> 00:09:27.852
PROFESSOR: Right it depends
on how many times steps
00:09:27.852 --> 00:09:29.320
I have [INAUDIBLE] my delta t.
00:09:29.320 --> 00:09:34.976
So it's something like t
final minus t 0 over delta t.
00:09:34.976 --> 00:09:36.350
That's the number
of times steps.
00:09:36.350 --> 00:09:38.230
But I have constant delta t.
00:09:38.230 --> 00:09:40.960
Times the number of
iterations that it
00:09:40.960 --> 00:09:44.920
takes to do the
nonlinear solve times n
00:09:44.920 --> 00:09:50.940
squared, which is a cost of
forming the Jacobian matrix.
00:09:50.940 --> 00:09:53.084
Now, this can be pretty large.
00:09:53.084 --> 00:09:55.000
So as we talked about,
this might be as big 10
00:09:55.000 --> 00:09:59.914
to the eighth in a bad case.
00:09:59.914 --> 00:10:02.330
Certainly be less than 10 to
the eighth because otherwise,
00:10:02.330 --> 00:10:03.621
you'll have numerical problems.
00:10:03.621 --> 00:10:05.600
Maybe, it's really
going to be 10
00:10:05.600 --> 00:10:08.910
to the sixth for a real problem.
00:10:08.910 --> 00:10:10.430
So you have a
million time steps.
00:10:10.430 --> 00:10:14.740
At each times step, you have
to solve the number of iter--
00:10:14.740 --> 00:10:16.310
the solve costs this much.
00:10:16.310 --> 00:10:17.780
Number of iterations
might be what?
00:10:17.780 --> 00:10:19.280
How many iterations
do you guys take
00:10:19.280 --> 00:10:21.634
to solve nonlinear equations?
00:10:21.634 --> 00:10:22.800
You solve a lot of them now.
00:10:22.800 --> 00:10:26.245
How many iterations
does it take?
00:10:26.245 --> 00:10:27.134
AUDIENCE: 10.
00:10:27.134 --> 00:10:28.800
PROFESSOR: 10, does
that 10 sound right?
00:10:32.084 --> 00:10:33.250
AUDIENCE: Yeah, [INAUDIBLE].
00:10:33.250 --> 00:10:35.260
PROFESSOR: You've got
an initial guess 10.
00:10:35.260 --> 00:10:37.177
If you have a medium
initial guess, maybe 100.
00:10:37.177 --> 00:10:39.635
It it's more than 100, you're
not going to converge, right?
00:10:39.635 --> 00:10:40.300
AUDIENCE: Yeah.
00:10:40.300 --> 00:10:41.133
PROFESSOR: Yeah, OK.
00:10:41.133 --> 00:10:43.025
So this is order of 10 maybe.
00:10:43.025 --> 00:10:45.700
So the million 10
and then n squared
00:10:45.700 --> 00:10:47.710
just depends on how
big your problem is.
00:10:47.710 --> 00:10:52.060
So in my group, the
problems we typically do,
00:10:52.060 --> 00:10:54.070
n is approximately 200.
00:10:54.070 --> 00:10:57.610
So this is like 200 squared.
00:10:57.610 --> 00:11:00.345
And then you can see that the
number of iterations, and this
00:11:00.345 --> 00:11:01.845
is maybe a little
bit overestimated.
00:11:01.845 --> 00:11:05.330
Maybe I can get away with 10 to
the fifth, something like that.
00:11:05.330 --> 00:11:08.350
But you can see it
starts to get pretty big.
00:11:08.350 --> 00:11:11.980
But for writing it, as you
guys are writing software,
00:11:11.980 --> 00:11:13.090
it's not that big a deal.
00:11:13.090 --> 00:11:17.180
So I can assign it as part 1 of
a three-part homework problem.
00:11:17.180 --> 00:11:18.270
It's the right ODE solver.
00:11:18.270 --> 00:11:20.346
It solves implicit equations.
00:11:20.346 --> 00:11:21.220
And you can write it.
00:11:21.220 --> 00:11:24.330
And the reason you can is
because you've already written
00:11:24.330 --> 00:11:25.330
methods that solve this.
00:11:25.330 --> 00:11:26.620
Actually, Matlab did it for you.
00:11:26.620 --> 00:11:27.703
But you write it yourself.
00:11:27.703 --> 00:11:29.380
I assigned to the
[? 10:10 ?] students.
00:11:29.380 --> 00:11:32.890
So you write the
[? backsub ?] solution
00:11:32.890 --> 00:11:34.780
to solve the equations.
00:11:34.780 --> 00:11:38.340
And we we've already
wrestled this.
00:11:38.340 --> 00:11:40.724
And I think you guys wrote
one of these yourself.
00:11:40.724 --> 00:11:42.140
And also, we can
use fsolve, which
00:11:42.140 --> 00:11:44.139
is just another little
bit better implementation
00:11:44.139 --> 00:11:45.040
of the same thing.
00:11:45.040 --> 00:11:49.000
And then now, you
write your next level.
00:11:49.000 --> 00:11:52.180
So you can see you're doing
a composite process, where
00:11:52.180 --> 00:11:54.234
you're doing something
at the top level.
00:11:54.234 --> 00:11:55.900
It's calling something
at a lower level.
00:11:55.900 --> 00:11:57.850
That's calling something
at a lower level.
00:11:57.850 --> 00:11:59.200
And then what you really
have to be concerned
00:11:59.200 --> 00:12:00.520
about if you're
worried about CPU time
00:12:00.520 --> 00:12:03.061
or whether the thing is going
to solve before the homework is
00:12:03.061 --> 00:12:07.062
due, is how big does this get.
00:12:07.062 --> 00:12:08.770
How many operations
are you really asking
00:12:08.770 --> 00:12:09.780
the computer to do?
00:12:09.780 --> 00:12:11.815
Now, remember the
computer is fast.
00:12:11.815 --> 00:12:13.690
Right, so it can do a
gigahertz or something.
00:12:13.690 --> 00:12:16.390
So you can do like 10 to the
eighth somethings per second.
00:12:19.320 --> 00:12:21.050
And what should
I say about this?
00:12:21.050 --> 00:12:27.290
I have not included,
but I should that it's--
00:12:27.290 --> 00:12:30.415
well, I'm assuming here
that the function evaluation
00:12:30.415 --> 00:12:34.080
is sort of order n of how
many variables they have.
00:12:34.080 --> 00:12:35.980
OK, sometimes,
function evaluations
00:12:35.980 --> 00:12:38.177
are much more
difficult than that,
00:12:38.177 --> 00:12:40.510
primarily because, sometimes,
they're actually including
00:12:40.510 --> 00:12:42.125
further nested stuff inside.
00:12:42.125 --> 00:12:44.500
So your f is actually coming
from some really complicated
00:12:44.500 --> 00:12:46.090
calculation itself.
00:12:46.090 --> 00:12:48.190
And then you have
another big giant factor
00:12:48.190 --> 00:12:50.949
of how much it costs to do f.
00:12:50.949 --> 00:12:53.240
But on this right here, I'm
just assuming it goes as n.
00:12:57.900 --> 00:12:59.680
So actually, let's
just multiply this out.
00:12:59.680 --> 00:13:02.330
So this is what, 10 to the
fourth, 10 to the 10th.
00:13:02.330 --> 00:13:04.320
So this is 10 to the 10th.
00:13:04.320 --> 00:13:08.550
10 to the 10th kind of flops.
00:13:08.550 --> 00:13:12.630
But your computer does 4 times
10 to the ninth per second.
00:13:12.630 --> 00:13:14.231
Is that right?
00:13:14.231 --> 00:13:16.240
Is that [? what they ?]
told us, four gigahertz?
00:13:16.240 --> 00:13:19.570
Yeah, so is this really not
bad with just a few seconds
00:13:19.570 --> 00:13:20.500
on your computer.
00:13:20.500 --> 00:13:22.570
So that's why you can do
these homework problems,
00:13:22.570 --> 00:13:23.960
and you still get them done.
00:13:23.960 --> 00:13:25.330
And if you have even
[? porous ?] limitations,
00:13:25.330 --> 00:13:26.440
it takes 1,000 times longer.
00:13:26.440 --> 00:13:28.856
Still, 1,000 [INAUDIBLE] take
20 minutes on your computer.
00:13:31.340 --> 00:13:34.640
So you're OK, but you can see
we're starting to get up there.
00:13:34.640 --> 00:13:37.580
And it's not that hard to
really make it a big problem
00:13:37.580 --> 00:13:39.320
if you do something wrong.
00:13:44.010 --> 00:13:51.740
And just a notation thing, a
lot of these things we're doing
00:13:51.740 --> 00:13:58.220
are writing software that
takes as input not just numbers
00:13:58.220 --> 00:13:59.889
or vectors.
00:13:59.889 --> 00:14:01.680
So those are things
that we call functions.
00:14:01.680 --> 00:14:03.890
Right, things that take a
number of vectors and input
00:14:03.890 --> 00:14:07.235
and give you, say, a function,
a number of vectors and output.
00:14:07.235 --> 00:14:09.860
But oftentimes, we actually want
to take the names of functions
00:14:09.860 --> 00:14:11.150
as inputs.
00:14:11.150 --> 00:14:13.880
So for example, the
root-finding problem,
00:14:13.880 --> 00:14:16.330
it takes as input what is f.
00:14:16.330 --> 00:14:19.570
Right, you have to tell it the
name of a function for fsolve
00:14:19.570 --> 00:14:20.890
to be able to call it.
00:14:20.890 --> 00:14:23.830
So fsolve is an object
that takes as input
00:14:23.830 --> 00:14:26.290
the name of a function.
00:14:26.290 --> 00:14:35.600
And things that do that
are called functionals
00:14:35.600 --> 00:14:43.870
So for example, this one is x
is equal to this root-finding
00:14:43.870 --> 00:14:48.295
functional, fsolve.
00:14:48.295 --> 00:14:50.380
And I'm going to write
with square brackets.
00:14:50.380 --> 00:14:51.742
It depends on f.
00:14:51.742 --> 00:14:52.450
That's it, right?
00:14:52.450 --> 00:14:55.630
[INAUDIBLE] f.
00:14:55.630 --> 00:14:57.550
So you give it the function f.
00:14:57.550 --> 00:14:59.060
And if fsolve's
good, it should be
00:14:59.060 --> 00:15:02.080
able to tell you what the x is.
00:15:02.080 --> 00:15:04.450
Now, in reality, we help it out.
00:15:04.450 --> 00:15:07.970
So we give it x 0 [? 2 ?]
just give it a better chance.
00:15:11.034 --> 00:15:13.700
I write the square brackets just
indicate that one of the things
00:15:13.700 --> 00:15:18.110
inside is a function,
not just numbers.
00:15:18.110 --> 00:15:20.010
We do this a lot.
00:15:20.010 --> 00:15:24.050
So you guys have done
this since you were kids.
00:15:24.050 --> 00:15:29.880
So derivatives Yeah, ddx of f.
00:15:33.720 --> 00:15:38.224
So we have, I don't know,
for example, grad of f.
00:15:41.160 --> 00:15:44.510
That's a functional, the grad.
00:15:44.510 --> 00:15:47.740
The symbol means the
operator, the functional,
00:15:47.740 --> 00:15:49.090
it operates on functions.
00:15:49.090 --> 00:15:50.980
You could write a
general gradient function
00:15:50.980 --> 00:15:54.472
that takes any scalar
valued function and computes
00:15:54.472 --> 00:15:56.180
the gradients with
respect to the inputs.
00:15:58.654 --> 00:16:00.820
This week, we give you one
that computes a Jacobian.
00:16:00.820 --> 00:16:02.230
Give it any f of x.
00:16:02.230 --> 00:16:05.710
It gives you the Jacobian
of f with respect to x.
00:16:05.710 --> 00:16:07.270
So that's a functional.
00:16:07.270 --> 00:16:09.001
So anyway, we do a
lot with functionals.
00:16:09.001 --> 00:16:10.500
You've been using
them all the time.
00:16:10.500 --> 00:16:12.730
I'm just trying get you
to think about abstractly.
00:16:12.730 --> 00:16:15.577
And then there's a whole
math about functionals.
00:16:15.577 --> 00:16:17.910
And there's a whole thing
called calculus of variations.
00:16:17.910 --> 00:16:20.080
It's all about, if the
objects you're worrying about
00:16:20.080 --> 00:16:22.300
are functionals, not
functions, then you
00:16:22.300 --> 00:16:24.739
have another whole bunch
of theorems and stuff
00:16:24.739 --> 00:16:26.530
you can do with that,
which you'll probably
00:16:26.530 --> 00:16:28.570
end up doing before you
get out of this place.
00:16:28.570 --> 00:16:30.195
But I'm not going to
do that right now.
00:16:32.740 --> 00:16:35.500
I would comment that one
kind of functional that's
00:16:35.500 --> 00:16:41.800
particularly important to a
lot of people in your research
00:16:41.800 --> 00:16:49.120
is the fact that the
energy of any molecule
00:16:49.120 --> 00:16:51.640
is a functional of
the electron density.
00:16:56.127 --> 00:16:57.710
And this is the basis
of what's called
00:16:57.710 --> 00:16:58.793
density functional theory.
00:17:09.160 --> 00:17:11.888
And so this is why a lot of
people, you and many of you
00:17:11.888 --> 00:17:14.304
included, are going to end up
using computer programs that
00:17:14.304 --> 00:17:17.752
are trying to find what the
electron density shape is,
00:17:17.752 --> 00:17:19.460
things like molecular
orbitals and stuff.
00:17:19.460 --> 00:17:21.490
You've probably
heard these before.
00:17:21.490 --> 00:17:23.560
And then there's the
functional that gives
00:17:23.560 --> 00:17:24.215
the energy of the molecule.
00:17:24.215 --> 00:17:26.006
And what you care about,
mostly, is energy,
00:17:26.006 --> 00:17:28.510
because that connects
to all your thermo.
00:17:28.510 --> 00:17:31.402
You're doing [? 1040, ?]
it's all about the energies.
00:17:31.402 --> 00:17:32.860
And once you know
all the energies,
00:17:32.860 --> 00:17:34.443
you can actually
compute the entropies
00:17:34.443 --> 00:17:35.730
and all kinds of stuff too.
00:17:35.730 --> 00:17:40.972
And so this is one application
that's really important.
00:17:40.972 --> 00:17:42.430
It was a little
surprising this was
00:17:42.430 --> 00:17:44.900
true and able to be
proven to be true.
00:17:44.900 --> 00:17:47.400
And so the guy who did that
got the Nobel Prize-- his name
00:17:47.400 --> 00:17:48.850
was Cohen-- a few years ago.
00:17:51.540 --> 00:17:53.365
Actually, in liquid
phase stuff, you
00:17:53.365 --> 00:17:55.180
work with Professor
Blankschtein.
00:17:55.180 --> 00:17:58.105
They also have a version of
this for properties of solution.
00:17:58.105 --> 00:17:59.980
And then where it's not
the electron density,
00:17:59.980 --> 00:18:03.100
but it's actually the density
of the material in the solution
00:18:03.100 --> 00:18:03.600
phase.
00:18:03.600 --> 00:18:04.891
Maybe Professor [? Swan ?] too.
00:18:04.891 --> 00:18:05.620
You guys do that?
00:18:05.620 --> 00:18:07.090
Yeah, so if you work with
Professor [? Swan, ?] you
00:18:07.090 --> 00:18:08.390
might learn about that too.
00:18:08.390 --> 00:18:10.110
So that's one where
the actual functional
00:18:10.110 --> 00:18:14.486
is like a big focus
of the entire project.
00:18:14.486 --> 00:18:16.110
And the key is that
this is a function.
00:18:16.110 --> 00:18:18.100
The density is a
function of the position.
00:18:18.100 --> 00:18:20.650
And then the functional
converts that into a number.
00:18:26.025 --> 00:18:27.310
That's all page one here.
00:18:29.850 --> 00:18:33.280
Now, I'm going to sort of
change topics, but it's related.
00:18:33.280 --> 00:18:37.360
It's about the connections
between numerical integration,
00:18:37.360 --> 00:18:39.969
interpolation, and
basis functions.
00:18:39.969 --> 00:18:42.010
These things are all
intimately tied up together.
00:18:42.010 --> 00:18:43.540
So I'll try to show you.
00:18:43.540 --> 00:18:46.240
So we try to solve a
numerical integral.
00:19:00.320 --> 00:19:03.170
And then let's look at
what we're really doing.
00:19:03.170 --> 00:19:04.970
When we're doing
this, we typically
00:19:04.970 --> 00:19:07.190
only do a few
function evaluations.
00:19:07.190 --> 00:19:12.260
So we have a few points where
we know the numbers f of tn.
00:19:12.260 --> 00:19:14.570
We've picked a few t's.
00:19:14.570 --> 00:19:16.190
We evaluated f.
00:19:16.190 --> 00:19:18.909
And from those few
points, that small number
00:19:18.909 --> 00:19:21.200
of function evaluations, we
want to get a good estimate
00:19:21.200 --> 00:19:22.712
of this whole integral.
00:19:22.712 --> 00:19:24.170
Now, it's a little
goofy, actually,
00:19:24.170 --> 00:19:25.970
if you think about what
we're trying to do.
00:19:25.970 --> 00:19:28.386
So suppose you're trying to
do an integral from negative 1
00:19:28.386 --> 00:19:29.134
to 1.
00:19:29.134 --> 00:19:31.550
We have some function that has
a point here, a point here,
00:19:31.550 --> 00:19:34.070
a point here.
00:19:34.070 --> 00:19:37.010
And just for concreteness,
let's say this is 1 over 26.
00:19:37.010 --> 00:19:38.000
And this is 1.
00:19:38.000 --> 00:19:39.724
And this is 1 over 26.
00:19:39.724 --> 00:19:41.390
And I have a particular
function in mind
00:19:41.390 --> 00:19:43.700
that has those three points.
00:19:43.700 --> 00:19:45.620
Now, the first thing
to keep in mind
00:19:45.620 --> 00:19:48.211
is there's an infinite
number of functions that
00:19:48.211 --> 00:19:49.460
go through those three points.
00:19:51.990 --> 00:19:55.850
Now, what we've
been doing is we've
00:19:55.850 --> 00:20:01.840
been fitting some
function to these points
00:20:01.840 --> 00:20:05.365
and then integrating that
function that we fitted.
00:20:05.365 --> 00:20:07.740
So when you're in high school,
you did the rectangle rule
00:20:07.740 --> 00:20:15.720
and what the function was
used was this and this.
00:20:15.720 --> 00:20:17.820
And you knew the
integral of this function
00:20:17.820 --> 00:20:20.220
and the integral
of this function.
00:20:20.220 --> 00:20:21.760
And you added them up.
00:20:21.760 --> 00:20:24.940
And you decided from
the rectangle rule
00:20:24.940 --> 00:20:28.019
that the integral was
equal to 1 and 1/26.
00:20:33.052 --> 00:20:34.760
But then some other
people in high school
00:20:34.760 --> 00:20:35.880
are a little smarter.
00:20:35.880 --> 00:20:38.262
And they said it's not
so smart [INAUDIBLE]..
00:20:38.262 --> 00:20:39.970
This thing was
symmetrical to begin with.
00:20:39.970 --> 00:20:41.595
And now, I'm fitting
it with a function
00:20:41.595 --> 00:20:42.890
that's wildly unsymmetrical.
00:20:42.890 --> 00:20:45.729
So let's instead not
evaluate those points.
00:20:45.729 --> 00:20:47.270
Let's instead evaluate
the midpoints.
00:20:47.270 --> 00:20:49.186
And maybe we find out
that the midpoint values
00:20:49.186 --> 00:20:51.060
are here and here.
00:20:51.060 --> 00:20:58.569
And so I would assume that the
function is just a flat line.
00:20:58.569 --> 00:21:00.610
That looks not so good
either, because it doesn't
00:21:00.610 --> 00:21:01.190
get through these points.
00:21:01.190 --> 00:21:02.690
But anyway, for the
two points I knew,
00:21:02.690 --> 00:21:03.770
if I only used those
two points, that
00:21:03.770 --> 00:21:05.240
would be the midpoint method.
00:21:05.240 --> 00:21:07.800
And I would get the
I midpoint method.
00:21:07.800 --> 00:21:08.870
And I'd get something.
00:21:08.870 --> 00:21:11.390
And I did it for this particular
function I have in mind.
00:21:11.390 --> 00:21:12.590
And it was 8 over 29.
00:21:15.149 --> 00:21:16.690
And then some other
people say, well,
00:21:16.690 --> 00:21:18.030
let's use the trapezoid rule.
00:21:18.030 --> 00:21:19.440
So let's approximate
the function
00:21:19.440 --> 00:21:25.310
is this over the interval,
which looks kind of sensible.
00:21:25.310 --> 00:21:26.810
And if you do that,
it turns out you
00:21:26.810 --> 00:21:29.101
get the same numbers as you
do with the rectangle rule.
00:21:29.101 --> 00:21:30.447
So I trapezoid.
00:21:34.830 --> 00:21:38.915
Now, let me just as a warning,
if you just do the rectangle
00:21:38.915 --> 00:21:40.290
rule and then the
trapezoid rule,
00:21:40.290 --> 00:21:41.915
and you got the same
number both times,
00:21:41.915 --> 00:21:44.082
you might be tempted to
conclude you have converged.
00:21:44.082 --> 00:21:45.873
And you have the true
solution, because you
00:21:45.873 --> 00:21:47.730
did two methods of
really different orders,
00:21:47.730 --> 00:21:49.524
and you got exactly
the same number.
00:21:49.524 --> 00:21:50.940
And you're super
happy, and you're
00:21:50.940 --> 00:21:53.622
ready to publish your paper.
00:21:53.622 --> 00:21:55.830
But of course, the function
might not look like this.
00:21:55.830 --> 00:22:00.510
So the particular function that
I always had in mind was f of t
00:22:00.510 --> 00:22:05.520
is equal to 1 over
1 plus 25 t squared.
00:22:05.520 --> 00:22:07.460
OK, that was just the
function I had in mind.
00:22:07.460 --> 00:22:16.616
And the way that function
looks like is like that.
00:22:16.616 --> 00:22:18.240
And the real integral
of that function,
00:22:18.240 --> 00:22:20.160
which you guys if you remember
how do it from high school,
00:22:20.160 --> 00:22:21.950
you can actually
evaluate that integral.
00:22:21.950 --> 00:22:26.520
It's nothing to do with
either 8/29 or 1 in 126.
00:22:26.520 --> 00:22:28.020
It's something
completely different.
00:22:28.020 --> 00:22:30.530
Now, those of you who went to
really advanced high schools,
00:22:30.530 --> 00:22:32.135
you also learned
Simpson's rule too.
00:22:32.135 --> 00:22:35.150
And Simpson's rule says
let's fit it to a parabola.
00:22:35.150 --> 00:22:38.930
And so that would give
you something like this.
00:22:38.930 --> 00:22:41.870
And if you do Simpson's
rule, I ended up getting this
00:22:41.870 --> 00:22:43.330
if I did the arithmetic right.
00:22:46.390 --> 00:22:51.420
1 and 1/3 plus 2/78.
00:22:51.420 --> 00:22:52.520
Who knew?
00:22:52.520 --> 00:22:55.070
OK, which also has
little relation
00:22:55.070 --> 00:22:57.750
to the true value
of the integral.
00:22:57.750 --> 00:23:00.000
Now, if you're smart, and
you did trapezoid and then I
00:23:00.000 --> 00:23:02.790
it to Simpson's, then you
say, oh my goodness, these
00:23:02.790 --> 00:23:03.990
are not really that similar.
00:23:03.990 --> 00:23:06.448
And so I'm not converged, and
I'd better do something else.
00:23:08.324 --> 00:23:11.800
So I just want you
to be aware when
00:23:11.800 --> 00:23:13.900
we do these rules, what
we're actually doing
00:23:13.900 --> 00:23:19.360
is we're fitting our points
to some continuous function
00:23:19.360 --> 00:23:21.040
that we know how to integrate.
00:23:21.040 --> 00:23:22.623
And then we do the
analytical integral
00:23:22.623 --> 00:23:25.570
of that function, which
is a fine thing to do.
00:23:25.570 --> 00:23:28.814
But just be aware that
it's like extrapolation.
00:23:28.814 --> 00:23:30.730
It's like we're imagining
what the function is
00:23:30.730 --> 00:23:34.610
doing between these points, when
we actually haven't tested it.
00:23:34.610 --> 00:23:35.305
Now, we could.
00:23:35.305 --> 00:23:36.596
We could just call more points.
00:23:36.596 --> 00:23:37.388
We could calculate.
00:23:37.388 --> 00:23:39.179
But that would cost us
more function evals.
00:23:39.179 --> 00:23:41.127
And we're trying to
save our computer time.
00:23:41.127 --> 00:23:42.710
So we want to finish
off before the TG
00:23:42.710 --> 00:23:43.900
so we can go buy some beer.
00:23:43.900 --> 00:23:51.370
And so we have to take some
choices about what to do.
00:23:51.370 --> 00:23:53.914
And we can't do an infinite
number of points for sure.
00:23:53.914 --> 00:23:56.080
All right, so in all these
methods, what we're doing
00:23:56.080 --> 00:23:59.110
is we're approximating
our true f
00:23:59.110 --> 00:24:04.300
with some approximate f, which
is equal to a sum of some basis
00:24:04.300 --> 00:24:12.657
functions like this.
00:24:12.657 --> 00:24:13.990
So we pick some basis functions.
00:24:13.990 --> 00:24:15.820
We know how to integrate.
00:24:15.820 --> 00:24:18.730
We fit, somehow determine
what these c's are.
00:24:18.730 --> 00:24:20.890
And then we decide
that the integral
00:24:20.890 --> 00:24:27.260
is equal to the
integral of f tilde, not
00:24:27.260 --> 00:24:29.360
the actual original f.
00:24:29.360 --> 00:24:31.690
And then that's actually
going to be the sum
00:24:31.690 --> 00:24:33.320
of ck times the integral.
00:24:37.047 --> 00:24:38.630
And we carefully
chose basis functions
00:24:38.630 --> 00:24:42.260
that we don't how to do
this integral; analytically.
00:24:42.260 --> 00:24:44.030
So we just plug
that in, and we just
00:24:44.030 --> 00:24:48.040
have to figure out that the
c's are to get this to work.
00:24:48.040 --> 00:24:49.750
OK, so that's what
you really did
00:24:49.750 --> 00:24:51.290
when you're doing this trapezoid
rule and Simpson's rule,
00:24:51.290 --> 00:24:52.630
and all that kind of stuff.
00:24:52.630 --> 00:24:54.671
And it's just different
choices of what the basis
00:24:54.671 --> 00:24:56.200
functions were that you used.
00:24:59.495 --> 00:25:01.120
Now, there's a lot
of choices about how
00:25:01.120 --> 00:25:04.720
to do that kind of approximation
to make f tildes that
00:25:04.720 --> 00:25:09.320
are sort of like f's with some
limited amount of information.
00:25:09.320 --> 00:25:12.140
And it generally has to do
with how many k's we have.
00:25:12.140 --> 00:25:15.845
So if k is less than n, then
it's really an approximation.
00:25:15.845 --> 00:25:17.380
Because there's no
way that we would
00:25:17.380 --> 00:25:20.690
expect that with just
a few parameters,
00:25:20.690 --> 00:25:22.780
we could get a lot f of tn.
00:25:22.780 --> 00:25:30.870
So we would expect that f tilde
of tn is not equal to f of tn,
00:25:30.870 --> 00:25:32.710
at least not at all the points.
00:25:32.710 --> 00:25:34.030
Because if n is--
00:25:34.030 --> 00:25:36.670
if we have a lot of points and
we only have a few parameters,
00:25:36.670 --> 00:25:40.900
we don't expect to be
able to fit it perfectly.
00:25:40.900 --> 00:25:43.410
So we know this
is what we expect.
00:25:43.410 --> 00:25:45.577
It might happen
if, by some luck,
00:25:45.577 --> 00:25:47.410
we chose a basis function
which was actually
00:25:47.410 --> 00:25:50.510
the true function f, then it
would actually fit perfectly.
00:25:50.510 --> 00:25:54.100
But most of the time, not.
00:25:54.100 --> 00:25:57.740
By Murphy's Law, in fact, it's
going to be completely wrong.
00:25:57.740 --> 00:26:00.440
But this is not such a bad idea.
00:26:00.440 --> 00:26:03.020
We may come back to this later.
00:26:03.020 --> 00:26:06.890
Another possibility is
to choose k equals to n.
00:26:06.890 --> 00:26:21.390
And then usually, we can force
f tilde of tn to equal f of tn.
00:26:21.390 --> 00:26:24.310
So at the points that we have,
if we have enough parameters
00:26:24.310 --> 00:26:26.860
equal to the number
of points, then often,
00:26:26.860 --> 00:26:27.880
we can get an exact fit.
00:26:27.880 --> 00:26:31.510
We can force the function to
go through all the points.
00:26:31.510 --> 00:26:35.771
And when we do this, this
is called interpolation.
00:26:43.480 --> 00:26:46.620
Now, when should you do
the one or the other?
00:26:46.620 --> 00:26:48.310
When should I make
k as big as n?
00:26:48.310 --> 00:26:50.170
Or when should I
use k less than n?
00:26:50.170 --> 00:26:52.030
A really crucial thing
is whether you really
00:26:52.030 --> 00:26:53.820
believe your points.
00:26:53.820 --> 00:26:56.050
If you really believe
you know f of tn
00:26:56.050 --> 00:26:58.780
to a very, very high number
of significant figures,
00:26:58.780 --> 00:26:59.872
then it make sense.
00:26:59.872 --> 00:27:01.330
You know that's
true, so you should
00:27:01.330 --> 00:27:05.170
force your fitting function
to actually match that.
00:27:05.170 --> 00:27:08.650
However, if you get the f of tn,
for example, from an experiment
00:27:08.650 --> 00:27:11.890
or from a stochastic simulation
that has some noise in it,
00:27:11.890 --> 00:27:14.410
then you might not
really want to do this.
00:27:14.410 --> 00:27:16.630
Because you'll be
fitting the noise
00:27:16.630 --> 00:27:18.530
in the experiment as well.
00:27:18.530 --> 00:27:24.280
And then you might want to use
a smaller number of parameters
00:27:24.280 --> 00:27:25.642
to make the fits.
00:27:25.642 --> 00:27:26.850
OK, so you can do either one.
00:27:26.850 --> 00:27:28.240
You're in charge.
00:27:28.240 --> 00:27:29.674
Right, you get the problem.
00:27:29.674 --> 00:27:31.840
If somebody tells you,
please give me this integral,
00:27:31.840 --> 00:27:33.610
they don't care how you do it.
00:27:33.610 --> 00:27:35.485
They just want you to
come back and give them
00:27:35.485 --> 00:27:38.986
the true value of the integral
to some significant figures.
00:27:38.986 --> 00:27:40.610
So it's up to you to
decide what to do.
00:27:40.610 --> 00:27:42.252
So you can decide what to do.
00:27:42.252 --> 00:27:43.960
What a lot of people
do in the math world
00:27:43.960 --> 00:27:46.630
is this one,
interpolation, where
00:27:46.630 --> 00:27:50.490
you're assuming that the f
of tn's are known exactly.
00:27:50.490 --> 00:27:53.560
And so your function
evaluating function is great.
00:27:53.560 --> 00:27:57.920
And in that case, you can write
down this nice linear equation.
00:27:57.920 --> 00:28:01.440
So we're trying to make f
of tn with a true function
00:28:01.440 --> 00:28:06.006
to be equal to ck phi n of t.
00:28:06.006 --> 00:28:07.469
Sorry, phi k of t.
00:28:12.260 --> 00:28:14.790
Which n in [INAUDIBLE]
k is equal to n.
00:28:14.790 --> 00:28:17.900
And I can rewrite this this way.
00:28:17.900 --> 00:28:20.446
This is like a set of numbers.
00:28:20.446 --> 00:28:21.570
So I can write as a vector.
00:28:21.570 --> 00:28:23.540
I'll call it y.
00:28:23.540 --> 00:28:25.520
And this thing has
two indices on it,
00:28:25.520 --> 00:28:27.300
so it looks like a matrix.
00:28:27.300 --> 00:28:30.384
So I'll call that m.
00:28:30.384 --> 00:28:34.620
And c is my vector of
unknown parameters.
00:28:34.620 --> 00:28:37.560
And wow, I know how
to solve this one.
00:28:37.560 --> 00:28:40.730
So c is equal to m backslash y.
00:28:46.200 --> 00:28:49.770
Now, for this to work and
have a unique solution,
00:28:49.770 --> 00:28:52.200
I need m not to be singular.
00:28:52.200 --> 00:28:54.090
And what that means is
that the columns of m
00:28:54.090 --> 00:28:55.530
have to be linearly independent.
00:28:55.530 --> 00:28:57.120
So what are the columns of m?
00:28:57.120 --> 00:29:10.330
They are phi k of t1,
phi k of t2, phi k of tn.
00:29:14.160 --> 00:29:21.315
And I want that this column is
not a sum of the other columns.
00:29:36.030 --> 00:29:38.100
So this is kind of
requirement for j
00:29:38.100 --> 00:29:40.542
not equal [? to k ?] This
is kind of a requirement
00:29:40.542 --> 00:29:42.250
if I'm going to do
this kind of procedure
00:29:42.250 --> 00:29:46.390
is I can't choose a
set of basis functions
00:29:46.390 --> 00:29:48.850
that's going to have one of
them be a linear combination
00:29:48.850 --> 00:29:49.810
of the other ones.
00:29:49.810 --> 00:29:50.980
Because then I'm going
to get a singular matrix.
00:29:50.980 --> 00:29:52.210
My whole process is going
to have a big problem
00:29:52.210 --> 00:29:53.540
right from the beginning.
00:29:53.540 --> 00:29:56.290
So I have to ensure
that this is not true.
00:29:56.290 --> 00:29:57.550
Yeah, that this is true.
00:29:57.550 --> 00:29:59.290
These are not equal
to some of them.
00:30:02.400 --> 00:30:14.570
Another way to write that is to
choose the the columns that phi
00:30:14.570 --> 00:30:16.869
to be orthogonal to each other.
00:30:16.869 --> 00:30:18.160
And so I can write it this way.
00:30:34.760 --> 00:30:37.502
OK, so I want them to be
orthogonal to each other.
00:30:37.502 --> 00:30:39.585
And then you can see how
this kind of generalizes.
00:30:39.585 --> 00:30:42.280
If I had a lot of t's,
this would sort of
00:30:42.280 --> 00:30:49.432
look like the integral of phi k
of t, phi j of t is equal to 0.
00:30:55.770 --> 00:30:58.770
And so this gives the idea
of orthogonal functions.
00:30:58.770 --> 00:31:03.490
So this is orthogonal functions,
and my discrete point's tn.
00:31:03.490 --> 00:31:07.410
And this is the general concept
of orthogonal functions.
00:31:07.410 --> 00:31:08.880
And then when I
define it this way,
00:31:08.880 --> 00:31:12.849
I need to actually decide my own
a's and b's that I want to use.
00:31:12.849 --> 00:31:14.640
And if you use different
ones, you actually
00:31:14.640 --> 00:31:16.015
have different
rules about what's
00:31:16.015 --> 00:31:17.280
orthogonal to each other.
00:31:17.280 --> 00:31:20.040
And also, sometimes,
people are fishy,
00:31:20.040 --> 00:31:22.030
and they [? spit ?]
a little w of t
00:31:22.030 --> 00:31:25.249
in here just to do, because they
have some special problem where
00:31:25.249 --> 00:31:26.790
this w keeps showing
up or something,
00:31:26.790 --> 00:31:27.831
and they want to do that.
00:31:27.831 --> 00:31:29.425
But that's up to them.
00:31:29.425 --> 00:31:30.800
All right, we
won't specify that.
00:31:30.800 --> 00:31:33.200
But this is the general idea
of making basis functions
00:31:33.200 --> 00:31:35.540
orthogonal to each other.
00:31:35.540 --> 00:31:36.980
And it's very
analogous to making
00:31:36.980 --> 00:31:38.300
vectors orthogonal
to each other, which
00:31:38.300 --> 00:31:39.133
is what we did here.
00:31:39.133 --> 00:31:41.140
These are actually vectors.
00:31:41.140 --> 00:31:43.430
Right, it's just as the
vector gets longer and longer,
00:31:43.430 --> 00:31:44.990
it gets more and more like
the continuous function,
00:31:44.990 --> 00:31:46.580
and we have more t points.
00:31:46.580 --> 00:31:48.420
And eventually, it
looks like that integral
00:31:48.420 --> 00:31:49.660
as you go to infinity.
00:31:52.040 --> 00:31:53.540
So anyway, so a lot
of times, people
00:31:53.540 --> 00:31:57.290
will choose functions which
are designed to be orthogonal.
00:31:57.290 --> 00:31:59.720
And one way to do it is just
to go with a continuous one
00:31:59.720 --> 00:32:00.219
to start.
00:32:02.690 --> 00:32:05.570
Now, we can choose any
basis functions we want.
00:32:05.570 --> 00:32:07.430
Over here, we already
chose a whole bunch.
00:32:07.430 --> 00:32:09.180
You guys have done this already.
00:32:09.180 --> 00:32:10.070
I don't know if you
knew you were doing it.
00:32:10.070 --> 00:32:11.570
But you chose a lot of different
basis functions completely
00:32:11.570 --> 00:32:13.260
different from each other.
00:32:13.260 --> 00:32:16.170
You can choose the
different ones if you want.
00:32:16.170 --> 00:32:21.680
And one popular
choice is polynomials.
00:32:21.680 --> 00:32:23.180
But it's not the only choice.
00:32:23.180 --> 00:32:25.755
And in fact, a little bit
later, we're going to,
00:32:25.755 --> 00:32:27.410
in the homework
problem six, I think
00:32:27.410 --> 00:32:28.640
we're going to have
a problem where
00:32:28.640 --> 00:32:30.530
you have to use a different
kind of basis functions.
00:32:30.530 --> 00:32:32.238
You can use any basis
functions you want.
00:32:32.238 --> 00:32:35.110
And there's kind of ones that
match naturally to the problem.
00:32:35.110 --> 00:32:38.930
However, people know a
lot about polynomials.
00:32:38.930 --> 00:32:42.770
And so they've decided
to use polynomials a lot.
00:32:42.770 --> 00:32:44.530
And so that's the
most common choice
00:32:44.530 --> 00:32:46.500
for numerical integration
is polynomials.
00:32:46.500 --> 00:32:48.870
One thing is polynomials
are super easy to integrate,
00:32:48.870 --> 00:32:51.036
because that was the first
thing you learned, right,
00:32:51.036 --> 00:32:52.700
how to do the intervals
of polynomials.
00:32:52.700 --> 00:32:54.507
And also, they're
easy to evaluate.
00:32:54.507 --> 00:32:56.090
It only takes n
operations to evaluate
00:32:56.090 --> 00:32:58.340
a polynomial of nth order.
00:32:58.340 --> 00:33:01.304
And there's a ton of
theorems, and things
00:33:01.304 --> 00:33:02.470
are known about polynomials.
00:33:02.470 --> 00:33:03.800
We have all kinds of
special properties.
00:33:03.800 --> 00:33:04.850
And we have all kinds
of convenient tools
00:33:04.850 --> 00:33:05.600
for handling them.
00:33:05.600 --> 00:33:08.710
So let's go with
polynomials for now.
00:33:08.710 --> 00:33:12.490
Now, I want to warn you, though,
that how you order polynomials
00:33:12.490 --> 00:33:17.230
are really not very
good fitting functions
00:33:17.230 --> 00:33:20.020
or interpolating functions
for a lot of problems.
00:33:20.020 --> 00:33:22.420
And in particular,
for this one I
00:33:22.420 --> 00:33:25.620
showed here where
this is the function.
00:33:25.620 --> 00:33:31.682
Here, that function, you can
try to fit that with higher
00:33:31.682 --> 00:33:32.890
and higher order polynomials.
00:33:32.890 --> 00:33:33.973
So I can have my function.
00:33:38.670 --> 00:33:41.700
The function value at 0 is 1.
00:33:41.700 --> 00:33:44.094
So the function looks
kind of reasonable.
00:33:47.060 --> 00:33:50.650
Like that, symmetrical,
better than I drew it.
00:33:50.650 --> 00:33:56.300
And you can put some
number of points in here.
00:33:56.300 --> 00:33:58.550
And then you can do an
interpolating polynomial that
00:33:58.550 --> 00:34:01.640
goes through all the points.
00:34:01.640 --> 00:34:07.400
And if you do it for this
one with five points,
00:34:07.400 --> 00:34:11.409
you'll get something
that looks like this.
00:34:11.409 --> 00:34:18.150
It goes negative, comes up
like this, goes negative again.
00:34:18.150 --> 00:34:21.570
So that's one you
get if you use five
00:34:21.570 --> 00:34:23.969
points across this interval.
00:34:23.969 --> 00:34:26.340
If you use 10
points, then you get
00:34:26.340 --> 00:34:41.960
one that looks like this
where the peak value here
00:34:41.960 --> 00:34:45.280
is two where as the max of the
original function's just one.
00:34:45.280 --> 00:34:47.170
And the original
function's always positive,
00:34:47.170 --> 00:34:49.520
but they both go negative.
00:34:49.520 --> 00:34:54.130
So even within the range
of the interpolation,
00:34:54.130 --> 00:34:56.590
the polynomials can
have wild swings.
00:34:56.590 --> 00:35:00.740
And so this first one was
for a fifth order polynomial,
00:35:00.740 --> 00:35:03.280
or maybe fourth
order, five points.
00:35:03.280 --> 00:35:06.630
And this one is for the
next 10th order or 9th order
00:35:06.630 --> 00:35:07.350
polynomial.
00:35:07.350 --> 00:35:09.760
I'm not sure which one.
00:35:09.760 --> 00:35:10.260
Yes.
00:35:10.260 --> 00:35:12.760
AUDIENCE: What points are these
polynomials trying to fit?
00:35:12.760 --> 00:35:15.430
Or does it seem
like [INAUDIBLE]??
00:35:15.430 --> 00:35:18.220
PROFESSOR: So I
didn't draw very well.
00:35:18.220 --> 00:35:19.960
When you guys go home, do this.
00:35:19.960 --> 00:35:23.100
Right, on Matlab it
will take you a minute
00:35:23.100 --> 00:35:26.720
to solve this equation where
you put in the function values
00:35:26.720 --> 00:35:27.870
here.
00:35:27.870 --> 00:35:33.670
And put in the polynomial
values that [INAUDIBLE]
00:35:33.670 --> 00:35:37.215
here for the monomials,
t to the first power,
00:35:37.215 --> 00:35:39.340
to the second power, to
the third power, like that.
00:35:39.340 --> 00:35:43.230
And just do it, and you can
get the [? pot ?] yourself.
00:35:43.230 --> 00:35:45.230
And just to convince
yourself about the problem,
00:35:45.230 --> 00:35:47.042
and this is very
important to keep in mind.
00:35:47.042 --> 00:35:49.250
Because we're so used to
using polynomials that think
00:35:49.250 --> 00:35:50.708
they're the solution
to everything.
00:35:50.708 --> 00:35:51.830
And Excel does them great.
00:35:51.830 --> 00:35:55.010
But it's not this
solution to everything.
00:35:55.010 --> 00:35:56.890
Nonetheless, we're going
to plow on and keep
00:35:56.890 --> 00:35:58.098
using them for a few minutes.
00:36:02.790 --> 00:36:08.200
So one possibility when
I'm choosing this basis set
00:36:08.200 --> 00:36:10.280
is I could choose
what are called
00:36:10.280 --> 00:36:17.260
monomials, which is phi k is
equal to x to the k minus 1.
00:36:21.370 --> 00:36:23.575
And this one is super simple.
00:36:26.200 --> 00:36:28.390
And so a lot of
times it satisfies
00:36:28.390 --> 00:36:32.170
this orthogonality thing,
which I showed you before.
00:36:32.170 --> 00:36:36.280
But it's a really
bad choice usually.
00:36:36.280 --> 00:36:38.480
And what it is is
that the matrix
00:36:38.480 --> 00:36:41.170
you have to solve for
this interpolation,
00:36:41.170 --> 00:36:46.960
if you have evenly-spaced
points and you have x to the k,
00:36:46.960 --> 00:36:49.330
it turns out that
this matrix usually
00:36:49.330 --> 00:36:52.120
has a condition number
that grows exponentially
00:36:52.120 --> 00:36:55.430
with the number of
terms in your basis,
00:36:55.430 --> 00:36:57.880
the number of how many
monomials you include.
00:36:57.880 --> 00:37:00.010
So condition numbers
that grow exponentially
00:37:00.010 --> 00:37:03.400
is really bad idea.
00:37:03.400 --> 00:37:05.252
So this is not what people use.
00:37:05.252 --> 00:37:06.710
So you might be
tempted to do this,
00:37:06.710 --> 00:37:08.540
and maybe you did it
once in your life.
00:37:08.540 --> 00:37:10.180
I won't blame you if you did.
00:37:10.180 --> 00:37:11.980
But it's not a
really smart idea.
00:37:11.980 --> 00:37:14.350
So instead, a lot of
the mathematicians
00:37:14.350 --> 00:37:16.270
of the 1700s and
1800s spent a lot
00:37:16.270 --> 00:37:19.040
of time worrying about this.
00:37:19.040 --> 00:37:24.250
And they'd end up deciding
to use other polynomials that
00:37:24.250 --> 00:37:26.050
have people's names on them.
00:37:26.050 --> 00:37:35.020
So a guy named
Lagrange, the same guy
00:37:35.020 --> 00:37:38.340
who did the multiplier.
00:37:38.340 --> 00:37:40.594
He suggested some polynomials.
00:37:40.594 --> 00:37:41.760
They're kind of complicated.
00:37:41.760 --> 00:37:43.140
They're given in the book.
00:37:43.140 --> 00:37:47.440
But what they really are is I'm
taking the monomial basis set,
00:37:47.440 --> 00:37:49.650
and I'm making a new basis set.
00:37:49.650 --> 00:37:58.430
So I want phi l to be equal
to the sum dlk of phi k,
00:37:58.430 --> 00:38:02.640
where this phi k
it the monomials.
00:38:02.640 --> 00:38:04.927
So really, you can
choose any d you want.
00:38:04.927 --> 00:38:06.510
You're just taking
linear combinations
00:38:06.510 --> 00:38:10.920
of the original monomial basis
set to make up polynomials.
00:38:10.920 --> 00:38:12.690
And if you cleverly
choose the d,
00:38:12.690 --> 00:38:14.610
you can choose it to
have special properties.
00:38:14.610 --> 00:38:16.860
So Lagrange gave
a certain choice
00:38:16.860 --> 00:38:20.220
of the d that is
brilliant in the sense
00:38:20.220 --> 00:38:23.130
that, when you
solve the problem m
00:38:23.130 --> 00:38:27.720
times c is equal to y to try to
do your interpolation, it turns
00:38:27.720 --> 00:38:34.440
out that the c's you want,
ck, is equal to f of tk.
00:38:34.440 --> 00:38:37.310
So you don't have to do any
work to solve what the c's are.
00:38:37.310 --> 00:38:39.930
It's just the f values directly.
00:38:39.930 --> 00:38:44.530
No chance for condition
numbers, no problems, that's it.
00:38:44.530 --> 00:38:47.790
So that's the
Lagrange polynomials.
00:38:47.790 --> 00:38:51.500
Now, Newton had
a different idea.
00:38:51.500 --> 00:38:55.105
So this is for
Lagrange polynomials.
00:39:00.352 --> 00:39:04.510
If you choose a Lagrangian
[? polynomial ?] basis,
00:39:04.510 --> 00:39:05.200
then this is it.
00:39:05.200 --> 00:39:06.650
And you do a square system.
00:39:06.650 --> 00:39:07.950
That's it.
00:39:07.950 --> 00:39:12.065
If you choose the Newton, even
he got into this business too.
00:39:12.065 --> 00:39:13.770
And he made some polynomials.
00:39:13.770 --> 00:39:16.800
And he made ones that
have a special property
00:39:16.800 --> 00:39:19.830
that if you solve it--
00:39:19.830 --> 00:39:21.270
so he chose a different d.
00:39:21.270 --> 00:39:23.460
You still have to do some solve.
00:39:23.460 --> 00:39:25.620
You get the solution.
00:39:25.620 --> 00:39:28.500
And then you decide that I
don't like that interpretation.
00:39:28.500 --> 00:39:30.060
I want to put some
more points in.
00:39:30.060 --> 00:39:32.130
I'm going to do a few
more function evaluations.
00:39:32.130 --> 00:39:33.420
I'm going to add some more.
00:39:33.420 --> 00:39:38.010
So from the first solve, I
have the values of c1, c2, c3,
00:39:38.010 --> 00:39:41.550
up to c the number of
points I had initially.
00:39:41.550 --> 00:39:44.510
And now, I want to add
one more data point,
00:39:44.510 --> 00:39:47.811
or add one more f of tn, one
more function evaluation.
00:39:47.811 --> 00:39:49.435
I want to figure out
with c of n plus 1
00:39:49.435 --> 00:39:51.330
is, what the next one is.
00:39:51.330 --> 00:39:54.420
It turns out that
it's his polynomials.
00:39:54.420 --> 00:39:56.310
The values these
c's don't change
00:39:56.310 --> 00:39:59.140
when I make this matrix a little
bit bigger by adding one more
00:39:59.140 --> 00:40:01.060
y.
00:40:01.060 --> 00:40:04.060
And I add one more parameter.
00:40:04.060 --> 00:40:05.830
And I make one more
row in the matrix,
00:40:05.830 --> 00:40:08.070
one more column in the matrix.
00:40:08.070 --> 00:40:10.360
When I do that, it's
a unique property
00:40:10.360 --> 00:40:12.675
that none of these
other c's change.
00:40:12.675 --> 00:40:15.300
The solution is going to be the
same, except for the n plus one
00:40:15.300 --> 00:40:17.530
[? level. ?] So then I can write
a neat, cute little equation
00:40:17.530 --> 00:40:18.613
just for the n plus 1 one.
00:40:18.613 --> 00:40:20.920
I don't have to mess
with the previous ones.
00:40:20.920 --> 00:40:23.461
And so I'm just adding a basis
function that sort of polishes
00:40:23.461 --> 00:40:25.350
up my original fit.
00:40:25.350 --> 00:40:28.180
So that's a nice property too.
00:40:28.180 --> 00:40:30.550
Whereas with the
Lagrange ones, you
00:40:30.550 --> 00:40:32.110
have to recompute everything.
00:40:32.110 --> 00:40:33.276
And all the c's will change.
00:40:38.739 --> 00:40:40.030
Actually, the c's won't change.
00:40:40.030 --> 00:40:40.738
I take that back.
00:40:40.738 --> 00:40:41.440
You add points.
00:40:41.440 --> 00:40:44.010
The c's won't change, but the
actual polynomials will change.
00:40:44.010 --> 00:40:45.720
The d will change.
00:40:45.720 --> 00:40:49.200
Because the d depends on
the choice of the time
00:40:49.200 --> 00:40:51.750
points in a complicated
way in Lagrange.
00:40:56.046 --> 00:40:57.920
So now, you start to
think, oh my goodness, I
00:40:57.920 --> 00:40:59.450
have all these choices,
I can use Lagrange.
00:40:59.450 --> 00:41:00.241
I could use Newton.
00:41:00.241 --> 00:41:01.550
I could use monomials.
00:41:01.550 --> 00:41:03.740
Maybe, I should go back to bed.
00:41:03.740 --> 00:41:06.912
And you think, well,
could someone just
00:41:06.912 --> 00:41:08.120
tell me what the best one is?
00:41:08.120 --> 00:41:08.810
And I'll just do that one.
00:41:08.810 --> 00:41:11.018
I don't need to learn about
all this historical stuff
00:41:11.018 --> 00:41:13.220
about development of
polynomial mathematics.
00:41:13.220 --> 00:41:15.450
I'd just like you to just
tell me what the answer is.
00:41:15.450 --> 00:41:20.940
So Gauss, typically, figured
out what the best on is.
00:41:20.940 --> 00:41:26.050
And he said, you know, if
we're going to do integrals,
00:41:26.050 --> 00:41:28.340
we'd really like to just
write our integrals this way.
00:41:28.340 --> 00:41:32.600
So our approximate value
to be equal to some weight
00:41:32.600 --> 00:41:38.340
functions times the function
evaluations at some times.
00:41:38.340 --> 00:41:40.170
Right, that's easy to evaluate.
00:41:40.170 --> 00:41:41.322
That would be good.
00:41:41.322 --> 00:41:43.780
If I pre-computed the weight
functions [? and ?] the times,
00:41:43.780 --> 00:41:45.560
then I put any
function f in here,
00:41:45.560 --> 00:41:47.530
I could just do this like zip.
00:41:47.530 --> 00:41:51.486
It's one dot product of the f
of tn's I'll have the solution.
00:41:51.486 --> 00:41:52.360
So I want to do that.
00:41:52.360 --> 00:41:55.390
So what's the best possible
choices of the w's and the t's
00:41:55.390 --> 00:41:57.660
to make this really be close?
00:41:57.660 --> 00:41:59.370
Well, the problem
is that there's
00:41:59.370 --> 00:42:00.420
an infinite number of
functions that somebody
00:42:00.420 --> 00:42:01.860
might want to integrate.
00:42:01.860 --> 00:42:02.910
I can't do them all.
00:42:02.910 --> 00:42:05.280
So he said, let's
just find the one that
00:42:05.280 --> 00:42:12.220
works for as many monomials as
possible to the highest order,
00:42:12.220 --> 00:42:14.790
as far out as I can go.
00:42:14.790 --> 00:42:20.310
And it turns out I have nw's.
00:42:20.310 --> 00:42:21.750
I have nt's.
00:42:21.750 --> 00:42:25.280
So I have two n parameters
that I can fiddle around with.
00:42:25.280 --> 00:42:28.270
And so I should get something
like two n's in my polynomials,
00:42:28.270 --> 00:42:30.560
maybe 2n minus 1.
00:42:30.560 --> 00:42:34.630
So what I want is that
I want to choose this
00:42:34.630 --> 00:42:39.060
so that for the monomial
1, and the monomial t,
00:42:39.060 --> 00:42:42.880
and the monomial t squared,
that this i is actually
00:42:42.880 --> 00:42:45.270
going to be exact.
00:42:45.270 --> 00:42:49.190
So what I want is
for the monomial 1,
00:42:49.190 --> 00:42:52.970
I want the summation of wn.
00:42:52.970 --> 00:42:55.760
My f for the monomial
1 is always 1.
00:42:55.760 --> 00:42:57.620
So it's times 1.
00:42:57.620 --> 00:42:59.600
I want this to equal
the integral from a
00:42:59.600 --> 00:43:06.330
to b of 1dt, which is b minus a.
00:43:06.330 --> 00:43:09.730
OK, so that gives me one
equation for the w's.
00:43:09.730 --> 00:43:15.310
Then for the next one,
I want the summation wn
00:43:15.310 --> 00:43:21.478
times tn to equal the
integral of a to b of t.
00:43:21.478 --> 00:43:28.570
dt, which is equal to b squared
over a number 2 minus a squared
00:43:28.570 --> 00:43:32.390
over 2 if I did my
calculus correctly.
00:43:32.390 --> 00:43:33.796
And so there's another question.
00:43:33.796 --> 00:43:35.420
I have another question
for these guys.
00:43:35.420 --> 00:43:37.919
And I keep going until I get
enough equations that determine
00:43:37.919 --> 00:43:39.395
all the w's and c's.
00:43:39.395 --> 00:43:40.895
All right, so I do
it for t squared.
00:43:40.895 --> 00:43:49.410
So I can wn tn squared is equal
to [INAUDIBLE] b t squared dt,
00:43:49.410 --> 00:43:52.260
and so on, OK.
00:43:52.260 --> 00:43:54.510
And so if you do
this process, you
00:43:54.510 --> 00:43:57.840
get the weights and
the points to use
00:43:57.840 --> 00:44:00.950
that will give you what's called
the Gaussian quadrature rule.
00:44:00.950 --> 00:44:05.470
Now, the brilliant thing about
this is that by this setup,
00:44:05.470 --> 00:44:08.770
it's going to be exact for
integrating any polynomial up
00:44:08.770 --> 00:44:10.641
to order 2n minus 1.
00:44:13.347 --> 00:44:17.990
So if I choose
three, here I'll be
00:44:17.990 --> 00:44:20.570
able to determine-- if
I choose three points,
00:44:20.570 --> 00:44:22.550
I'll have six parameters.
00:44:22.550 --> 00:44:27.898
So I'll be able
determine up to the 11--
00:44:27.898 --> 00:44:29.814
11th order polynomial
would integrate exactly.
00:44:33.610 --> 00:44:36.482
No, that's wrong.
00:44:36.482 --> 00:44:38.480
n's only 3.
00:44:38.480 --> 00:44:41.370
I have six parameters.
00:44:41.370 --> 00:44:42.420
I can't do this.
00:44:42.420 --> 00:44:44.601
2 times 3 minus 1, what's that?
00:44:44.601 --> 00:44:46.800
5, there we go.
00:44:46.800 --> 00:44:50.180
To the fifth order
polynomial integrate exactly,
00:44:50.180 --> 00:44:52.580
I'll have an [? error ?]
for sixth order polynomial.
00:44:52.580 --> 00:44:56.280
So this is a rule.
00:44:56.280 --> 00:45:02.160
So before, if I gave you three
points like here, most of you
00:45:02.160 --> 00:45:03.720
probably would use
Simpson's rule.
00:45:03.720 --> 00:45:05.295
Fit this to a parabola.
00:45:05.295 --> 00:45:09.030
It has three parameters, a,
ax squared plus bx plus c,
00:45:09.030 --> 00:45:10.955
right, for a second
order polynomial.
00:45:10.955 --> 00:45:11.580
I would fit it.
00:45:11.580 --> 00:45:12.330
Then I'd integrate that.
00:45:12.330 --> 00:45:14.710
That would be the normal thing
that a lot of people would do.
00:45:14.710 --> 00:45:15.510
I don't know if
you would do that.
00:45:15.510 --> 00:45:17.402
But a lot of people
would do that.
00:45:17.402 --> 00:45:19.110
So that would be that
what you would get.
00:45:19.110 --> 00:45:20.100
And what [? error ?]
would you get?
00:45:20.100 --> 00:45:21.766
You'd have error that
would be the scale
00:45:21.766 --> 00:45:25.110
of the width of the intervals
to the fourth power.
00:45:25.110 --> 00:45:28.430
Right, because Simpson's
rule exact up to the cubic.
00:45:28.430 --> 00:45:31.035
You guys remember
this from high school?
00:45:31.035 --> 00:45:33.160
I don't know if you remember
delta t to the fourth.
00:45:33.160 --> 00:45:34.170
Yeah, OK.
00:45:34.170 --> 00:45:41.260
So Simpson's rule has an order
of delta t to the fourth.
00:45:41.260 --> 00:45:48.530
But Gauss's method has order
of delta t to the sixth.
00:45:48.530 --> 00:45:51.810
So this is for three points.
00:45:51.810 --> 00:45:54.060
In Gauss's method, you can
extend it to as many points
00:45:54.060 --> 00:45:56.550
as you want.
00:45:56.550 --> 00:45:58.220
All right, so
that's kind of nice.
00:45:58.220 --> 00:45:59.890
You only do three
points, you get
00:45:59.890 --> 00:46:02.160
[? error ?] of the sixth order.
00:46:02.160 --> 00:46:04.030
So that seems pretty efficient.
00:46:04.030 --> 00:46:05.305
So people kept going.
00:46:05.305 --> 00:46:08.490
And so the very,
very popular one
00:46:08.490 --> 00:46:10.860
use, actually, seven
points for Gauss.
00:46:14.940 --> 00:46:20.270
So you'd be able to
get to the n minus 1.
00:46:20.270 --> 00:46:23.220
So you'd be exact for 13.
00:46:23.220 --> 00:46:26.832
So if you have delta t to
the 14th power as your error,
00:46:26.832 --> 00:46:28.450
that's pretty good.
00:46:28.450 --> 00:46:36.760
And then what people often do is
they'll add eight more points.
00:46:36.760 --> 00:46:39.220
And they make another method
that uses all these points.
00:46:39.220 --> 00:46:44.894
It's called
[? Conrad's ?] method.
00:46:44.894 --> 00:46:46.810
He made another method
that gives you the best
00:46:46.810 --> 00:46:49.476
possible solution if you had the
15 points, given that 7 of them
00:46:49.476 --> 00:46:52.090
are already fixed by
the Gauss procedure.
00:46:52.090 --> 00:46:56.340
And so you can get the
integral using Gauss's method
00:46:56.340 --> 00:46:57.740
and using [? Conrad's ?] method.
00:47:07.480 --> 00:47:11.590
So you can compute the i using
Gauss's method, the seven
00:47:11.590 --> 00:47:14.122
points, the i with the
[? Conrad ?] method.
00:47:14.122 --> 00:47:15.580
I don't know how
to spell his name,
00:47:15.580 --> 00:47:18.260
Actually I think it's like this.
00:47:18.260 --> 00:47:19.930
15 points.
00:47:19.930 --> 00:47:22.580
And then you can take the
difference between these two.
00:47:22.580 --> 00:47:24.080
And if these two
are really similar,
00:47:24.080 --> 00:47:26.121
then you might think your
integral's pretty good.
00:47:26.121 --> 00:47:27.896
Because they're pretty
different methods,
00:47:27.896 --> 00:47:29.770
and if they get the same
number, you're good.
00:47:29.770 --> 00:47:31.690
And people have done a
lot of numerical studies
00:47:31.690 --> 00:47:33.231
in this particular
case, because it's
00:47:33.231 --> 00:47:35.350
used a lot in actual practice.
00:47:35.350 --> 00:47:39.750
And it turns out
that the [? error ?]
00:47:39.750 --> 00:47:45.760
is almost always less than this
quantity to the 3 Gauss power.
00:47:45.760 --> 00:47:48.480
So you have an error estimate.
00:47:48.480 --> 00:47:51.760
So you do 15
function evaluations.
00:47:51.760 --> 00:47:54.770
You get a pretty accurate
value for your integral,
00:47:54.770 --> 00:47:59.480
maybe, if you can be fit well
with whatever order polynomial,
00:47:59.480 --> 00:48:02.120
13th order polynomial.
00:48:02.120 --> 00:48:04.850
13th order polynomial, yeah.
00:48:04.850 --> 00:48:08.190
And so any function
that can be fit pretty
00:48:08.190 --> 00:48:09.690
well with the 13th
order polynomial,
00:48:09.690 --> 00:48:11.700
you should get pretty good
value in the integral.
00:48:11.700 --> 00:48:14.210
And you can check by comparing
to the [? Conrad ?] estimate,
00:48:14.210 --> 00:48:15.550
which uses a lot more points.
00:48:15.550 --> 00:48:17.950
And if two of these guys
give the same number
00:48:17.950 --> 00:48:19.550
or pretty close,
then you're pretty
00:48:19.550 --> 00:48:21.027
confident that you're good.
00:48:21.027 --> 00:48:23.360
And there's even a math proof
about this particular one.
00:48:23.360 --> 00:48:26.810
So this is like really
popular, and it's
00:48:26.810 --> 00:48:28.310
a lot better than
Simpson's rule,
00:48:28.310 --> 00:48:29.600
but it's kind of complicated.
00:48:29.600 --> 00:48:31.100
And part of it you
have to recompute
00:48:31.100 --> 00:48:35.000
these tn's and wn's for the
kind of problem you have.
00:48:35.000 --> 00:48:35.666
Yeah.
00:48:35.666 --> 00:48:41.360
AUDIENCE: If you're
using 15 points already,
00:48:41.360 --> 00:48:45.440
why not just split
your step size smaller?
00:48:45.440 --> 00:48:48.185
And if you did that, how
would the error [INAUDIBLE]
00:48:48.185 --> 00:48:49.075
the simplest?
00:48:49.075 --> 00:48:50.117
Use a trapezoid.
00:48:50.117 --> 00:48:51.700
PROFESSOR: Trapezoid
or Simpson's rule
00:48:51.700 --> 00:48:52.370
of something like that.
00:48:52.370 --> 00:48:52.590
Yes.
00:48:52.590 --> 00:48:54.381
AUDIENCE: How would
the error with those 15
00:48:54.381 --> 00:48:59.450
steps using trapezoid compare
to [? Conrad ?] and Gauss?
00:48:59.450 --> 00:49:02.190
PROFESSOR: So we should
probably do the math.
00:49:02.190 --> 00:49:06.410
So say we did Simpson's
rule, and we have delta t,
00:49:06.410 --> 00:49:11.960
and we divide it
by 15 or something.
00:49:11.960 --> 00:49:14.120
Then we can figure
out what this is,
00:49:14.120 --> 00:49:17.919
but it's still going to have
the scaling to the fourth power.
00:49:17.919 --> 00:49:19.460
So the prefactor
will be really tiny,
00:49:19.460 --> 00:49:22.320
because it will be 1/15
to the fourth power.
00:49:22.320 --> 00:49:24.239
But the scaling is bad.
00:49:24.239 --> 00:49:26.530
So that if I decide that my
integral's not good enough,
00:49:26.530 --> 00:49:29.780
and I need to do more points as
I double the number of points
00:49:29.780 --> 00:49:33.180
to get better accuracy, my
thing won't improve that much.
00:49:33.180 --> 00:49:35.430
Because it's still only going
to use the fourth power,
00:49:35.430 --> 00:49:36.590
whereas the Gauss'
one say it's going
00:49:36.590 --> 00:49:38.330
use the sixth power if I
only use the three points.
00:49:38.330 --> 00:49:40.204
Where we use 15 points,
is going to get scale
00:49:40.204 --> 00:49:43.110
to some 16th power or 15th
power or some really high power.
00:49:43.110 --> 00:49:47.880
So for doing sequences, the
Gauss' one is a lot better.
00:49:47.880 --> 00:49:51.110
However, what people
do is they usually
00:49:51.110 --> 00:49:52.730
pre-solve the
system of equations.
00:49:52.730 --> 00:49:55.064
Because this system of
equation does not depend.
00:49:55.064 --> 00:49:56.480
When you do it
with the monomials,
00:49:56.480 --> 00:49:58.310
it doesn't depend what f is.
00:49:58.310 --> 00:50:01.527
You can figure out the optimal
t's and w's ahead of time.
00:50:01.527 --> 00:50:03.110
And so they [INAUDIBLE]
tables of them
00:50:03.110 --> 00:50:05.900
depending on where
you want to go.
00:50:05.900 --> 00:50:08.510
People recomputed them, and
then you can get them that way.
00:50:11.217 --> 00:50:12.050
Time is running out.
00:50:12.050 --> 00:50:14.060
I just want to jump up.
00:50:16.820 --> 00:50:19.020
Now, for single integrals,
I didn't have to tell you
00:50:19.020 --> 00:50:19.520
and of this.
00:50:19.520 --> 00:50:22.960
Because you could have just
solved it with your ODE solver.
00:50:22.960 --> 00:50:26.870
So you can just go
to ode15s or ode45,
00:50:26.870 --> 00:50:28.400
and you'll get good enough.
00:50:28.400 --> 00:50:30.080
It won't be as
efficient as Gaussian.
00:50:30.080 --> 00:50:31.044
But who cares.
00:50:31.044 --> 00:50:32.210
You still get a good answer.
00:50:32.210 --> 00:50:33.440
You can report to your boss.
00:50:33.440 --> 00:50:35.810
It'll take you three
minutes to call
00:50:35.810 --> 00:50:38.099
it e45 to integrate the
one-dimensional integral.
00:50:38.099 --> 00:50:40.640
The real challenge comes with
the multidimensional integrals.
00:50:40.640 --> 00:50:44.580
So in a lot of problems, we
have more than one variable
00:50:44.580 --> 00:50:47.600
we want to integrate over.
00:50:47.600 --> 00:50:51.380
And so if you have an
integral, say, of a
00:50:51.380 --> 00:50:58.620
to b from the lower bound
of x to upper bound of x.
00:50:58.620 --> 00:51:03.209
So this is dx, dy, f of xy.
00:51:03.209 --> 00:51:04.750
This is the kind of
integral that you
00:51:04.750 --> 00:51:06.489
might need to do sometimes.
00:51:06.489 --> 00:51:08.530
OK, because you might have
two variables that you
00:51:08.530 --> 00:51:10.180
want to integrate over.
00:51:10.180 --> 00:51:15.070
And so how we do this, the best
way to think about it, I think,
00:51:15.070 --> 00:51:22.362
is that f of x is equal
to the integral of l
00:51:22.362 --> 00:51:27.695
of x u of x dy f of xy.
00:51:30.600 --> 00:51:32.524
And then, if I
knew what that was,
00:51:32.524 --> 00:51:34.690
I would evaluate that say
with a Gaussian quadrature
00:51:34.690 --> 00:51:35.231
or something.
00:51:35.231 --> 00:51:40.070
So I'd say that this i 2d is
approximately equal to the sum
00:51:40.070 --> 00:51:41.540
of some weights--
00:51:41.540 --> 00:51:44.150
I'll label them x just to remind
you where they came from--
00:51:46.680 --> 00:51:52.190
times f of xn.
00:51:52.190 --> 00:51:53.910
And I would use, say,
Gaussian quadrature
00:51:53.910 --> 00:51:56.243
to figure out what the best
value of the x's and the y's
00:51:56.243 --> 00:51:57.110
are to use.
00:51:57.110 --> 00:51:59.780
And then I would actually
evaluate this guy
00:51:59.780 --> 00:52:01.730
with a quadrature two.
00:52:01.730 --> 00:52:08.800
So I'd say f of xk or xn is
equal to the integral from l
00:52:08.800 --> 00:52:17.100
of xn to u of xn dy f of xn y.
00:52:17.100 --> 00:52:20.260
And so this itself, I would
give it another quadrature.
00:52:20.260 --> 00:52:24.390
This is going to be
equal some other weights.
00:52:24.390 --> 00:52:32.495
nj f of xn yj.
00:52:32.495 --> 00:52:34.910
Is that all right?
00:52:34.910 --> 00:52:37.130
And so this is
nested inside this.
00:52:37.130 --> 00:52:40.350
So now, I have a double loop.
00:52:40.350 --> 00:52:43.240
So if I took me, I
don't know, 15 points,
00:52:43.240 --> 00:52:45.740
15 function evaluations to get
a good number for my integral
00:52:45.740 --> 00:52:47.930
in one dimension,
then it might take
00:52:47.930 --> 00:52:53.240
me 15 squared function
evaluations to get
00:52:53.240 --> 00:52:55.240
a pretty good integral for two.
00:52:55.240 --> 00:52:57.260
15 might be a little
optimistic, so maybe 100
00:52:57.260 --> 00:52:59.210
is more like if you
really want a good number.
00:52:59.210 --> 00:53:02.577
So I need 101 dimension, 100 in
the second dimension squared.
00:53:02.577 --> 00:53:04.910
Because I subdivided the
integral into six equal pieces,
00:53:04.910 --> 00:53:06.020
say.
00:53:06.020 --> 00:53:08.030
And then now, so this is OK.
00:53:08.030 --> 00:53:08.980
So you can do this.
00:53:08.980 --> 00:53:10.640
It's only 10,000
function evaluations.
00:53:10.640 --> 00:53:13.630
No problem, you've got a
good gigahertz computer.
00:53:13.630 --> 00:53:15.260
Now, in reality,
in a lot of them,
00:53:15.260 --> 00:53:17.100
we do have a lot more
dimensions than that.
00:53:17.100 --> 00:53:20.420
So in my life, I do a lot of
electronic structure integrals.
00:53:20.420 --> 00:53:25.937
We have integrals that are dx 1d
y1 d z1 for the first electron
00:53:25.937 --> 00:53:27.020
interacting with a second.
00:53:30.670 --> 00:53:34.820
And then something that
I'm trying to integrate.
00:53:34.820 --> 00:53:37.920
It depends on, say, the
distance between electron one
00:53:37.920 --> 00:53:38.730
and electron two.
00:53:41.250 --> 00:53:46.386
OK, so this is
actually six integrals.
00:53:46.386 --> 00:53:48.720
And so instead of being
a 2 in the exponent,
00:53:48.720 --> 00:53:51.700
this is going to be like
a 6 in the exponent.
00:53:51.700 --> 00:53:54.350
So now, this is 10 to the 12th.
00:53:54.350 --> 00:53:57.270
That means that just evaluating
one integral this way
00:53:57.270 --> 00:53:59.740
might take me 100
seconds if I can evaluate
00:53:59.740 --> 00:54:02.026
the functions like that.
00:54:02.026 --> 00:54:04.150
But in fact, the functions
I want to evaluate here,
00:54:04.150 --> 00:54:06.100
these guys will take
multiple operations
00:54:06.100 --> 00:54:07.910
to evaluate the function.
00:54:07.910 --> 00:54:09.790
And so this turns
out I can't do.
00:54:09.790 --> 00:54:12.980
I can't even do one integral
this way with its six integrals
00:54:12.980 --> 00:54:14.270
deep.
00:54:14.270 --> 00:54:18.490
And so this is the
general problem
00:54:18.490 --> 00:54:19.834
of multidimensional integration.
00:54:19.834 --> 00:54:21.500
It's called the curse
of dimensionality.
00:54:21.500 --> 00:54:24.680
Each time you add one more
integral, you're nested there.
00:54:24.680 --> 00:54:28.220
All of a sudden, the effort
goes up tremendously and more,
00:54:28.220 --> 00:54:30.410
a couple of orders
of magnitude more.
00:54:30.410 --> 00:54:33.050
And so we'll come back
later in the class
00:54:33.050 --> 00:54:35.150
about how to deal
with cases like this
00:54:35.150 --> 00:54:37.070
when you have high
dimensionality cases to try
00:54:37.070 --> 00:54:38.390
to get better scaling.
00:54:38.390 --> 00:54:40.240
Instead of having it
being exponential,
00:54:40.240 --> 00:54:47.060
it was really 100 to the n right
now, 100 to the n dimensions.
00:54:47.060 --> 00:54:50.070
[? Can we ?] figure out
some better way to do it.
00:54:50.070 --> 00:54:52.330
All right, thank you.