NARRATOR: The
following content is
provided by MIT OpenCourseWare
under a Creative Commons
license.
Additional information
about our license
and MIT OpenCourseWare
in general
is available at ocw.mit.edu.
PROFESSOR: Finite difference
methods for initial value
problems that we're
coming to the end of,
and the solving large
systems that we're
coming to the beginning of, was
to talk today about matrices,
because that language is
just useful for everything.
So about the homeworks, Mr.
Cho put in a long weekend
grading the homeworks
that were turned in Friday
and preparing code
and output to go up
on the web page, probably
late tonight or tomorrow.
So have a look to see and I
hope that those codes will
be useful for the future too.
Right.
So we'll maybe say
more about the outputs.
And about the projects,
which by the way,
could grow out of
that homework or could
go in a totally
different direction,
I'm thinking that the right
time to say projects due
would be after spring break.
So pretty much nearly
immediately after the spring
break would be -- and
we'll talk about it more,
but just so you have
an idea of what's --
what timetable I had in mind.
So, matrices then.
In particular, finite
difference matrices.
So that's the second
difference matrix, K,
and I'll frequently use that
letter K as I did in 18.085
for that very, very
important and useful matrix.
So what are its properties?
It's tridiagonal.
That's a very important
property, which we'll see,
because that means
that computations
solving linear systems are very
fast with a tridiagonal matrix.
It's symmetric.
All its eigenvalues
are positive,
so I would say it's
symmetric, positive definite.
And those eigenvalues
and eigenvectors --
so the eigenvectors for this
matrix turn out to be discrete
-- not quite discrete
exponentials, discrete sines,
discrete sine function.
If I change the
boundary conditions,
I can get discrete cosines
as the eigenvectors.
Or I could get
discrete exponentials
as the eigenvectors
by making it periodic.
Maybe I'll mention how
to make it periodic
and then I'll erase it again.
To make it periodic means that
this second different centered
at point 1 should
look ahead to point 2
and look behind to point 0, but
that'll be the same as point n,
so I would put a minus 1
in that corner and this one
similarly, it looks ahead,
which really brings it around
again, since we're
sort of on a circle,
brings it around again here.
So that matrix now, with
minus 1's added in the corner,
I would call C, a circular.
So I'll leave the
letter K there,
but the right letter is C
while these minus 1's are here
to make it periodic.
By the way, they mess
up that the tridiagonal.
It's no longer tridiagonal.
Also, it certainly got
-- it's very sparse.
Mentioning sparse
reminds me, if you're
coding large matrices
that are sparse,
you should let MATLAB
know that they're sparse.
So MATLAB has a whole -- it
carries out operations --
if you tell it it's
a sparse matrix,
then it only operates where the
non-zeros are located and it
doesn't waste its time looking
at these 0's, these 0's through
all the matrix steps.
So using sparse MATLAB is
important thing to know about.
It's just typically
an s or an sp
will appear in MATLAB commands.
For example, I'll just
maybe fill it in here.
What's the sparse
identity matrix?
The normal identity
matrix would be
eye of n and the sparse identity
matrix is sp, speye of n.
Similarly, we would create K --
we could use 2 times speye of n
as the diagonal
of K and the two,
the upper and lower diagonals --
and we could tell it these two
entries.
Another thing to say -- what
I said I wasn't going to do,
I'll do because I hate to see
a K up there while it's not
right.
Two more points
about this matrix.
It's singular now.
The determinant is 0.
Now, I'm never going to take the
determinant of a giant matrix.
That's a bad thing to do.
Much better to recognize that
there's a vector, x, let's say,
and it's the vector of all 1's.
It's the vector of n 1's.
If you imagine multiplying this
matrix by the vector of n 1's,
what do you get?
You get all 0's.
So that vector of all
1's is in the null space,
I would say, of a matrix.
Null space is just
the vectors that
get wiped out by the matrix.
C*x is all 0's.
So that vector -- this matrix
has some nonzero vectors in its
null space.
I know right away then,
its determinant is 0.
So the determinant
of that is 0 and now
that would tell me something
about the eigenvalues
of the matrix.
It tells me about
one eigenvalue.
It's 0.
A matrix, like C, that
has C*x equals 0, well,
I could also say that
that's C*x equals 0*x.
That would actually be better,
so that both sides are vectors.
So I'm realizing that that
vector x in the null space
is an eigenvector and the
corresponding eigenvalue is 0.
Otherwise, the eigenvalues
will all still be positive.
So this would be a
positive semidefinite.
I would call that matrix
positive semidefinite,
the semi telling me that it
isn't quite definite, that it
gets down and has 0,
the matrix is singular.
But still, it's good to know
where all the other n minus 1
eigenvalues are.
They're all positive.
About the bandwidth --
let me go back to K now.
Because the bandwidth, strictly
speaking, the bandwidth of C
is very large.
The bandwidth is, if I'm
looking for only one number,
that number tells me
how many diagonals
I have to grow until I
reach the last nonzero.
So the bandwidth here is large.
And in general, the
operation count,
the amount of work
to do, is going
to grow with the bandwidth.
Of course, that
having full bandwidth
isn't quite the full
story for this matrix,
because it's so sparse.
I've got hundreds
of zero diagonals
in between and just this one.
Anyway, this is a matrix
with a large bandwidth,
but a little deceptive.
Now here's a matrix with
a -- back to K again.
I call the bandwidth
just 1 here.
Really, maybe half bandwidth
would be a better word.
The bandwidth is the number
of diagonals above or below --
take the maximum of the count
above and the count below
and in this case,
both counts are 1.
One diagonal above, one
diagonal below, so really,
half bandwidth, I
would say, is 1.
Just some convention
is needed there.
The crucial point is that
the bandwidth measures
the amount of work to do
when you do elimination,
as MATLAB will do, of course.
One other thing about MATLAB
-- so I'm often referring
to MATLAB and I'm thinking
of its backslash command.
The backslash, which solves A*x
equal b by just A backslash b.
So if it doesn't know
these matrices are sparse,
it will go through
all the steps,
not taking advantage of the fact
that we've got all these 0's
here.
If it does know
that they're sparse,
then it's tremendously fast.
Let me come back to that point.
Maybe actually backslash
is smart enough
to look to see
whether the matrix,
whether sparseness is available.
So I shouldn't have said --
I think maybe there's a lot
engineered into backslash.
Actually, backslash,
also engineered in there
is the least squares solution.
If you give it a
rectangular problem
and, say, too many equations,
so that you can't expect
to have an exact
solution, backslash
will pick the least
squares solution.
And much more.
That's probably the
most used operation.
So I said something about
the eigenvalues and everybody
sees that if I multiply this
matrix by the values of u
at successive mesh points, I'll
get the second difference that
corresponds to u_xx.
Now this would correspond
-- this matrix --
so that tells me I'm
taking a finite difference.
That plus tells me I'm going
in the forward direction,
so that's 1 at the mesh
value to the right minus 1
of the mesh value
at the center point.
Of course, this has
to be divided by delta
x and that by delta x squared.
So these are matrices, along
with the backward difference
and the centered
difference, out of which
you build the basic
finite difference
equation, as you've done.
I want to make a comment
here though, now,
on this topic of eigenvalues.
Eigenvalues for
K really tell you
the truth about the
matrix K. The eigenvalues
of that matrix K start
just a little above 0
and they go to a
little before 4.
So the eigenvalues for K
-- shall I put that here?
The eigenvalues of K
are between 0 and 4.
They come very close
to 4 and quite close
to 0, depending on the size
of the matrix of course.
Let me just do -- if I
took the one by one case,
its eigenvalue is 2.
If I took the two by two case,
its eigenvalues are 1 and 3.
Notice nice properties there.
The eigenvalues 1 and 3
are positive, as we said.
This matrix K has positive
eigenvalues, whatever size.
What's more, the 1 and 3
kind of interlace the 2.
So what I'm saying is,
the eigenvalue for that
is in between the eigenvalue
for the two by two,
and the two
eigenvalues 1 and 3 are
in between the three eigenvalues
that I would get for the three
by three case.
Maybe I'll just
write those down.
They are useful numbers.
So K_2 has lambda equal 1 and 3.
The three by three one,
I think the eigenvalues
are 2 minus root 2,
which is smaller than 1;
2, which is in between;
and 2 plus root 2, which
is larger than 3.
And of course, they're
all between 0 and 4.
Just a comment.
How do I know they're
between 0 and 4?
There's a somewhat
handy little rule
for getting the location
of eigenvalues, that's
just worth knowing as a
sort of general principle,
but of course it can't tell
you exactly where they are.
First of all, the fact that
the matrix is symmetric
tells us what about
the eigenvalues?
So we learn a very,
very important fact
about the eigenvalues from
just looking at the matrix
and observing that
it's symmetric.
That tells us that the
eigenvalues are real.
They're real numbers.
Actually, it tells us
something equally important
about the eigenvectors.
The eigenvectors of the
matrix are orthogonal.
The symmetric matrix
has real eigenvalues,
orthogonal eigenvectors.
That's a little bit of
linear algebra to know.
Now, why between
0 and 4, though?
There's there's this
-- what's his name?
Gershgorin.
Gershgorin pointed out --
and it's a two-line proof --
that every eigenvalue
is in one of his --
one or more of his circles.
So where are the
Gershgorin circles?
A typical Gershgorin circle is
centered at the number -- here,
2 -- that's on the diagonal,
and its radius is the sum off
the diagonal, but you
take absolute values.
Gershgorin wasn't doing
any of the fine points that
really locate the eigenvalue.
So in this matrix,
all the centers
are at 2, the
diagonals, and the radii
are 2, also 2, because that and
that in absolute value make 2.
This first row makes a 1 and
actually, that's what -- that,
by a little careful argument,
is what gets the eigenvalues not
actually touching
0 or touching 4.
It takes a little patience.
All you could say from here --
all Gershgorin could say would
be, the eigenvalues are
in a circle centered at 2,
its radius is 2, so
they're between 0 and 4.
Then it's that first
and last row --
either the first or the last
row would do to say that
they're strictly positive.
Of course, the true
second difference,
I should divide K
by delta x squared.
So the eigenvalues will be
divided by delta x squared,
divided by a small number,
so they will really
be much larger.
The ratio of the
biggest to the smallest
isn't affected by
the delta x squared
and that's a key number.
Let me put that maybe on
the next board, for K.
So lambda_max is about 4.
Lambda_min, the smallest
eigenvalue, is close to 0.
How close?
It's of the order of some
constant over n squared.
Now I use the word
condition number and that --
so let me write that down.
Condition number -- say c
of K for condition number --
is the ratio lambda_max
over lambda_min.
Of course, I'm
using here the fact
that K is a symmetric matrix.
Symmetric matrices
are so beautiful
that their eigenvalues give
you a reliable story here.
So 4 divided by this -- the
main point is that it's O of --
it's of order 1 over n -- sorry.
When I divide by
lambda_min, that
puts the n squared
up in the numerator.
It's O of n squared,
growing like n squared.
And that condition
number, somehow it
measures how close the
matrix is to being singular,
because it involves
this lambda_min, which
is the smallest eigenvalue and
would be 0 if it were singular,
and it's scaled by lambda_max.
So if I'm multiply
the matrix by 100,
what happens to the
condition number?
No change.
No change, because I'm doing
lambda_max over lambda_min.
Both lambda_max and lambda_min
would be multiplied by 100,
but the ratio
wouldn't be different
So this is a useful measure.
And O of n squared, that's
a big number if n is large.
A big numbers if n is large.
And that condition number, where
does it come in elimination?
In elimination, the
round-off error -- roughly,
the rule of thumb
is that you would --
if the condition number
is 10 to the eighth,
you might lose eight significant
bits in the back slide.
You could.
So this condition number
measures how sensitive
your matrix is to round-off.
So that's a few
thoughts about matrices
and that matrix K in particular.
Now what about the other
matrix, the first difference?
The point I want to
make about that matrix
is, what about its eigenvalues?
What are the eigenvalues of
that upper triangular matrix?
They are, if you remember
linear algebra --
but I can just tell you
quickly the main point --
for a triangular matrix,
the values are sitting
on the diagonal.
So this matrix
has the eigenvalue
minus 1 repeated n times.
That true fact is totally
misleading, totally misleading.
The eigenvalues for
this triangular matrix
don't even notice what
I've got above the diagonal
and somehow they can't give
a reasonable picture of what
the matrix is actually doing.
So maybe that's my
warning here, that
for a matrix which is
absolutely not symmetric, right?
I mean, not at all symmetric.
For the centered
difference, which is --
what's the centered difference?
I was going to say symmetric,
but it's the opposite.
Centered difference would be --
let's put delta_centered down
here.
Centered difference
would be, I'd have a 1.
0's would go on the point
that -- for the central value.
1 would multiply the
forward value and minus 1
would multiply that,
and then I'd have 1/2
and then I'd probably
have a delta x.
But the main point is, my matrix
would look something like this.
Minus 1's and 1's
on two diagonals.
Now we could find the
eigenvalues of that matrix.
Do you know anything
about the eigenvalues?
This is a chance
for me just to speak
for a few minutes
about useful facts
that you can tell about a
matrix just by looking at it.
So what do I see when
I look at that matrix?
Is it symmetric?
It's the opposite of symmetric,
because the symmetric matrix
up there, if I transpose
it, it doesn't change.
If I transpose this
one, I'll get some kind
of a backward difference.
If I transpose this one, then
the 1's and the minus 1's will
reverse.
I'll get the negative.
So the rule for this centered
difference is -- so shall I --
how am I going to call
centered difference?
Del centered, for the moment.
The transpose of
that is minus itself.
So what does that tell me?
First, that matrix is -- it's
the opposite of symmetric,
but it's actually OK.
What I mean by OK is,
its eigenvectors --
we're back to
orthogonal eigenvectors.
I didn't say anything about
the eigenvectors of del plus,
but actually, that was
the biggest problem.
This matrix del plus has one
eigenvalue repeated n times,
and it has only one
eigenvector, not --
it doesn't even have a
full set of eigenvectors,
much less orthogonal ones.
So that matrix is like --
you don't want to trust
the eigenvalue picture that you
get from a matrix like that.
Here this anti-symmetric
matrix can be trusted.
Its eigenvalue
picture is reliable.
It does tell you
what's going on.
The eigenvectors are orthogonal.
They're complex, actually.
Actually, they'll look a lot
like our e to the i*k*x's So we
don't panic when we see
complex eigenvectors.
The eigenvalues are -- do you
know what the eigenvalues looks
like for an
anti-symmetric matrix?
They're pure imaginary, just
the way that when we took second
differences -- maybe I'll
just put here the centered
difference, the centered
first difference,
when we applied it to -- I want
apply to e to the i*k*x to find
what factor comes out.
So I get e to the i k plus 1
delta x from the plus side.
This is the matrix I'm
doing, with 1 and minus 1.
So the minus 1 will give me
minus e to the i k minus 1
delta x and of course, I
factor out of that the e
to the i*k*x's so I'm
left with e to the --
the thing that factors out
is e to the i delta x minus e
to the minus i delta
x -- and what's that?
That multiplies the e to
the i*k*x, the eigenvector.
This is like the
eigenvalue, and what
do I say about that quantity?
Of course, it's 2*i sine
delta x, pure imaginary.
And I should have
divided by the 2, which
would take away that 2.
So it's pure imaginary.
In reality and in
this Fourier analysis,
both are giving this
understanding of what i is.
So that when we did this sort
of operation, von Neumann's
rule of following
the exponential,
we got something reasonable.
When we do it with
this one, it's
von Neumann that's reliable
and the eigenvalues
that are not reliable.
So the eigenvalues of this
being all minus 1's is nonsense,
doesn't tell us what the forward
difference operator is really
doing.
But von Neumann tells us
what is truly going on.
Of course, that would
be the same thing
in which this minus
1 wouldn't appear
and the 2 wouldn't appear.
So what would we get
out of von Neumann?
So this is for delta plus.
I would factor out an e
to the i delta x minus 1.
That's what would multiply
it e to the i*k*x.
Oh, e to the i -- sorry.
Yes, that's right.
That would multiple
e to the i*k delta x.
Sorry, should've
had delta x there,
but I wasn't paying
attention to that.
Here I was paying attention to
this and it was pure imaginary.
Here I'm paying
attention -- here,
von Neumann at least is paying
attention to this and what's he
seeing?
Not pure imaginary
or purely real,
off in the complex
plane and that's
really what the right growth
factor or the right number
to associate with frequency k
for this forward difference.
So I guess I'm saying
that von Neumann does --
did the right thing to
come up with these pictures
for the growth factors.
Eigenvalues confirm that
and really pin it down
when they're reliable.
Eigenvectors and -- eigenvalues
are reliable when eigenvector
are orthogonal and that's for
matrices that are symmetric
or anti-symmetric or -- there's
a little bit larger class
of matrices that includes
orthogonal matrices,
but beyond that --
actually, there's been a lot
of discussion over many years
of eigenvalues and the --
for problems that are not
controlled by symmetric
or anti-symmetric matrices.
The alternative, the
more refined idea
of pseudo-eigenvalues is now
appearing in an important book
by Trefethen and Embree,
with many examples -- OK,
I won't pursue that.
Right.
So this is some basic fact about
those matrices, all of which
are one-dimensional.
Now this board prepares the
way to get into 2D and 3D.
So I just want to
ask, what does the --
we didn't really do the heat
equation or the wave equation
in 2D, but we could have.
The von Neumann test
would be straightforward,
but now I want to think
about the matrices.
What does the two-dimensional
second difference matrix
look like?
What I'm going to do,
just to look ahead,
I'm going to use MATLAB's
operation called kron,
short for Kronecker, to
create a 2D matrix out
of this 1D centered difference.
So if I think now about
centered differences,
second differences -- so I'm
approximating u_xx plus u_yy --
or really, I'm approximating
minus u_xx minus u_yy,
because that K approximates
minus the second difference.
What do I -- what do
my matrices look like?
What's their bandwidth?
How expensive is
it to invert them?
These are the key questions.
What's the matrix K2D
that corresponds to --
gives me second differences
in the x direction plus second
differences in the y direction.
So I'll write K2D
for that matrix.
Let's get some picture of it.
First of all, let
me imagine that I'm
on a square with a square grid.
Delta x in both directions.
Square grid, let me say
N mesh points each way.
So N, I don't know whether
I'm counting -- right now,
I won't worry whether I'm
counting the boundary ones
or not.
Probably not.
So N -- I'll say N point, point
N, N points -- so N delta x,
and in this direction N delta x.
So my matrix is a border.
It's of size N squared,
the number of unknowns.
Now I will be a
little more careful.
Here, let me take the
boundary values as given.
They're not unknown.
So N in this picture is 4.
One, two, three, four unknowns
there on a typical row.
Now I have to give them a new
number -- five, six, seven,
eight, nine, 10, 11, 12, 13,
14, 15, 16 -- and N being 4,
N squared is 16 for
that particular square.
So my matrix is 16 by 16.
But somehow I want to be able
to create it out of 4 by 4
matrices like K. So K1D, which
I'm just going to call K --
K_N will be the 4 by 4 one.
2 minus 1 -- so that's the
matrix that gives me second
differences along a typical
row or down a typical column.
But now what am I looking for?
I'm looking to do both --
second differences in a row
and a column.
So if I pick a typical
mesh point, like number 11.
Mesh point 11 -- let me
blow up this picture here.
It's going to be influenced,
mesh point 11, by 10 and 12,
the second differences in the
x direction, and by 7 and 15 --
notice those are not
so close to 10 or 11 --
the second differences
in the y direction.
So let me blow that up.
So here's mesh point number 11
corresponding to row 11 of K2D.
So I guess I'm asking what
row 11 of K2D will look like.
So here's mesh points 10 and 12,
so I have a second difference
in the x direction from
u_xx, for minus u_xx, that
means I can put a minus 1 there,
a 2 there and a minus 1 there.
Now I have the same in the
y direction with 15 and 7,
so I have a minus 1 in column
15, a minus 1 in column 7,
so I have four minus 1's and
then 2 more for the center
gives me a 4.
So a typical row will have --
it'll be sparse, of course --
it'll have a minus 1 in position
7, a minus 1 in position 10,
a 4, a minus 1, and a minus
1 over there in position 15.
That's a typical row of K2D.
It adds to 0.
If Gershgorin got his
hands on this matrix,
he would say that since all
the rows -- the 4 will be on --
will the 4 be on the diagonal?
Yes, the 4 will be on
the diagonal all the way.
I guess -- let's see a
little bit more clearly what
the matrix K2D looks like.
This is probably the
most studied matrix
in numerical analysis because
it's a model of what stays nice
and what gets more difficult as
you move into two dimensions.
So some things
certainly stay nice.
Symmetries -- so properties
of K2D, K2D will be --
it'll be symmetric again.
It'll be positive
definite again.
So what does K2D -- it's 16 by
16 and a typical row looks like
that.
That's a typical interior row.
What does maybe -- if
I took the first row,
what does the very
first row look like?
If I put row one above it.
Let me draw row one.
So what's the
difference with row one?
It's these are boundary values.
These are not -- these are going
to show up on the right-hand
side of the equation.
They're known numbers.
They're not -- they don't
involve unknown things,
so the only neighbors
are 2 and 5.
So I'll still have the second
difference, minus 1, 2,
minus 1, and minus
1, 2, minus 1 --
still this five-point
molecule with 4 at the center,
but there won't be
so many neighbors.
There'll be the neighbor at
the -- just to the right,
neighbor number two, and
there'll be neighbor number
five a little further
along and that's it.
It's like the 2
minus 1 boundary row.
It's a boundary row
and a boundary row
hasn't got as many minus
1's because it hasn't got
as many neighboring unknowns.
That boundary row would
have three neighbors.
Row two would have a 4,
minus 1, minus 1, but not --
nobody there.
So now I'm going
to try to draw K2D.
Let me try to draw K2D.
I can do it.
K2D will have, from the u_xx,
the second differences along
the rows, it will have --
I'll have a row, row one,
it'll have another K on row two.
It'll have a K on
-- a K for each row.
These are blocks now.
That's of size N by
N. All of those are,
so the whole thing is
N squared by N squared.
Actually, I'll stop here.
If I wanted to create that
matrix with the same K on --
it's somehow the identity is
in there and the matrix K is
in there and that's
-- let's see.
It's one or the other of those.
I guess it's the first one.
So what's this
Kronecker product?
Now I'm saying what
this construction is.
It's very valuable, because it
allows you to create matrices.
If you created K as
a sparse matrix and I
as a sparse matrix, then
the Kronecker product
would be automatically dealt
with as a sparse matrix.
So what's the rule for
a Kronecker product?
You take the first matrix,
which is the identity.
And the rest 0's.
So that's I. Then
-- so that's 1D --
and then you multiple each --
each entry in I becomes a block
with entry a_(i, j) times this
guy, K. This'll be the 0 block,
this'll be the K block,
this'll be the K block.
I take those 1's and multiply
K and those 0's and multiple K
and that's what I get.
So you see now, the Kronecker
product is a bigger guy.
If matrix A was p by p
and matrix B was q by q,
then the Kronecker product
would be p times q by p times q.
It would be square again.
It would be symmetric if
A and B are symmetric.
Actually, it would have
various nice properties.
Its eigenvalues
would be the products
of the eigenvalues of A
times the eigenvalues of B.
It's a very handy
construction and here we
saw it in a pretty easy case
where A was only the identity.
Well, let me do this case.
What does the Kronecker
rule produce here?
So I take this
first matrix and I
put it in here, 2, minus 1,
minus 1, 2, minus 1, minus 1,
2.
So that's the matrix K. And
now, each of those numbers
multiplies the second thing,
which here is I. So it's 2I,
minus I, minus I, 2I,
minus I, minus I, minus I,
minus I, 2I, and 2I.
Now that is the second
difference matrix
that does all the columns.
When I add those -- so now
I'm going to add that to that.
They're both size n squared.
They give me K2D.
So the neat construction of
K2D is Kronecker product kron
of I, K plus kron of K, I.
That's the matrix and let's look
at what it actually looks like.
We're seeing it block-wise here.
We saw it row-wise here.
And maybe now we can take one
more look at it, assemble it.
So K2D -- I plan to add that
matrix to that matrix to get
K2D.
So it will have -- since they
both have 2's on the diagonal,
it'll have 4's on the
diagonal, so 4's all the way,
but it'll be block-wise.
Up here, this
block is K plus 2I.
So that block is
-- goes down to 4.
Now, the K contributed minus
1's next to the diagonal,
and I guess that's it, right?
K plus 2I has that block as
-- that's the one, one block.
Why is it -- so we're seeing
-- did I get it right?
Yes, so we're seeing the
neighbors to the right
and left, but now
let me bring in --
only now comes the neighbor
above or below and that comes
from this off-diagonal
block, minus I,
which is then minus
1's to minus 1's.
These are all 0 blocks.
Here will be a minus
the identity block.
Then another block,
the diagonal blocks
are all the same,
and another one,
another minus the identity.
Now we're getting to
interior mesh points,
where we see typical rows.
So a typical row has
the 4's on the diagonal,
the minus 1's left and right,
and the minus 1's far left
and far right -- and of course,
what I'm going to ask you is,
what's the bandwidth?
We're coming to
the key point now.
What's the bandwidth
of this matrix?
I only have two nonzeros above
the diagonal on a typical row,
but I have to wait n diagonals
before I get to the second one.
So the bandwidth is n, because
I have to wait that long.
Then the operation count -- if
I just do ordinary elimination
on this matrix, the operation
will be the size of the matrix
times the bandwidth squared.
We can easily check that that's
the right count of operations
for a banded matrix.
This is operations
on a banded matrix.
So what do we get?
The size of the
matrix is N squared.
The bandwidth is N,
and it gets squared,
so we get N to the fourth.
It's getting up there.
If N is 1,000, we've got a
serious-sized matrix here.
Still a very sparse matrix,
so it's not like give up,
but the question is, does
the matrix stay sparse
as we do elimination?
That's the center
of the next lecture.
Can we organize elimination
-- how closely can we organize
elimination to preserve all
these zillions of 0's that are
between here?
They're easy to preserve
down here, way up there,
but in this intermediate, those
diagonals tend to fill in,
and that's not a
happy experience.
And it's even less happy in 3D.
So let me just do this
same calculation for 3D
and then we're done for today.
So K3D -- you might like to
think how K3D could be created
by the kron operation.
Let me just imagine it.
So what's K3D looking like?
I've got three directions,
so I have a 6 in the center
and six minus 1's
beside it and so 6 goes
on the diagonal of K3D.
Six minus 1's go
on a typical row
and how long do I have to wait
until I reach the last one?
So I'm again going to do
the size -- which is --
the matrix size is
going to be like N cube,
and what's the bandwidth
going to be like?
Maybe I'll ask you now.
What do you think
for the bandwidth?
How long -- if I just number
it in the simplest way --
and that'll be
the key next time,
is there a better numbering?
But if I just number along rows
until the rows fill up a plane
of rows and then I
move up to the next,
in the z direction to the plane
above, then I have to wait --
this was the x and y, so this
gave me a bandwidth of N,
but that's not the bandwidth
because I have to wait until I
finish the whole plane, until
I go up to the next plane
and catch this chap, so
the bandwidth is N squared.
Then the operations, which are
size times bandwidth squared,
we're up to N^7 and that is
a really horrifying exponent
to see, N to the seventh.
That means that even
for a moderate --
a problem in 3D with
a moderate number,
say 1,000 unknowns at
each direction, we have --
we're looking, in theory, at
a cost that we can't afford.
Of course, there's
a lot to do here.
So there's a lot to do if
I stay with direct methods
and that's what I'll do for
the next couple of lectures.
Then it's really in 3D
that iterative methods
become essential.
I mean, this one, I could
do those by direct methods,
2D, by a smarter direct method
than any I've tried here,
but in 3D, even
smarter elimination
is facing a serious exponent
and loses to iterative methods.
So that's what's
coming, actually.
So today's lecture, you
see, with the simplest
possible matrices, what
the central questions are.
Thanks and I've got homeworks
coming back from Mr. Cho
and I'll collect any that are
ready to come in and see you
Wednesday.
