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: OK, so I'll start with this.
This is a method for a particular class of problem
and maybe moving interfaces, moving curves, curves that are
-- or surfaces in higher dimensions that are being
transported, and it's their shape that matters.
So it's a sort of topic within shape differential equations
rather than particle differential equations.
And some cases of this problem can be studied by a
particularly fast method, developed here at MIT in
Course 6 actually, called the fast marching method.
It's really thanks to Per-Olof Persson that there's now quite
a bit on the 18.086 web page.
The basic equations that'll be in the lecture and also
demonstrations, quick demos that show the point.
And I better mention Sethian and Osher, who jointly
proposed the method, conceived of this approach in the late
'80s, and each has written a book, or coauthored a book, on
the level set method.
So those books are, so to speak, the
authorities on the method.
OK, so what's the idea?
Well, here is the curve, the shape, the
interface that's moving.
And let's suppose --
well, imagine it's sort of -- there's a wall a fire here,
and think of it as a moving fire front.
Then that fire front will move, normal to itself, to
some later position at time 1.
We could normalize --
in a perfect case when it was homogeneous and the speed of
movement didn't change with position, we could imagine
that the speed is 1.
So where is this fire front here at time 1?
It contains all the points that are reached from the
initial fire in unit time.
And since we're going at speed 1, let's say, this is all the
points at distance 1.
This curve, in this best, simplest case, will be the
curve that's everywhere a distance 1
away from this curve.
Now, how to describe the curves is the key point here.
Shall I describe the curve by particles?
I could imagine these points as being my unknowns, my
positions, and some equation for following them.
So I follow those ordinary differential equations, one
for every point, and they reach new positions at time 1.
And then I interpolate the next curve, the new
curve at time 1.
Well, that's not what we do.
That's not the idea of the level set method, and I'll
explain why that idea of following particles, which is
a totally natural one, comes to grief
numerically too often.
This picture, of course, was like an ideal case, and the
particles separated a little bit, but not too much.
Just to say in the word, the trouble with the particle idea
is that particles could separate wildly, so you had a
big, big gap to fill in like a fan.
Or they could come very close together so that you had total
instability numerically.
You couldn't compute any differences because the
distance was too small, and that would be like a shock.
So we're really connecting this type of problem.
It connects well with the conservation laws that we've
been studying.
Fans and shocks appear, and I'll show you right away.
Can I just say, what is the approach?
What does this word "level set" mean?
Well, the level set of a function --
I'm saying what this word level set is, and that's going
to be the way the curve is described, sort of implicitly.
Instead of saying where the points are that are on the
curve, the curve is being given by an equation like phi
of x, y equals 0, or phi of x, y equal 1, or phi
of x, y equals 7.
Those three curves would be level sets of phi.
Do you see the idea of levels sets?
They're the solutions to the equation phi
of x, y equal constant.
Let me take a simple example here on the
board that I can --
so the idea of level sets, first of all.
Let me take an easy -- a particularly nice function.
phi of x, y equals, say, x squared plus y squared.
OK, so I want to ask, what are the level
sets of that function?
So I set phi of x, y equals a constant to get the c level
set, you could say.
And, of course, we recognize that as a circle when c is
positive, and there is no set when c is negative.
So this would be -- the level sets would be circles.
So that would be one description.
I mean, what we're really doing is we're introducing a
new dimension to the problem.
Here, it's a 2D problem with x square and y squared, but I'm
bringing in the third dimension z and thinking of
the curve in the x, y plane as being produced by -- so what
would be the graph of z equals x squared plus y squared?
That's some surface now in 3D.
So I have x, y, and now I've introduced this third
dimension in order to see a whole surface.
And now I cut that surface by the plane z equal to c.
So z equal to c cuts through that surface by a plane here.
And that cuts out a circle, which I can drop down to the
x, y plane if I want to see it down there, and that's the
level set at height c.
So that's a way to describe the curve, but it's sort of --
it's implicit.
You'll see that it has major advantages when -- maybe I can
just describe them without drawing them.
Suppose I have a problem with two fires, say, two circles of
fire, spreading.
What's going to happen as time evolves with these
two circles of fire?
Well, having taken circles, they're just going to expand,
so after a little time, I've got two bigger circles.
But after a little more time, the circles will meet and form
one noncircle, right?
A sort of figure-eight shape.
Now that's a change in topology, and the world of
shape calculations is always having to ask how am I going
to deal with changes in topology when curves meets,
when unconnected things become connected, or when connected
things disconnect?
And the one great value of the level set approach is it deals
with that easily.
I could imagine a double -- so I don't know that I'll succeed
to do it here, but I could imagine a function, which --
OK, maybe this example, this picture, isn't too bad.
So that's intended to be a surface.
Well, let me make the surface go up that way.
OK, so when I cut it at z equals c, that might cut the
surface in two different planes.
So the level set at a certain level is two circles.
There's two circles.
But now, let z increase, and what happens?
You see how naturally, as z increases, these two circles
are coming closer and closer, and at a
certain point, they merge.
And phi of x and y remains a totally nice function.
I mean, I could create some probably fourth-degree
function curve that would do that fine.
I could graph it.
I could visualize it.
And the point is that by going up that one dimension and then
slicing at different levels, you allow these things to
happen smoothly.
OK, so that's one feature of level sets, that description.
OK, so that's the way curves are going to be defined.
So our unknown becomes this function, which will become a
function of x, y and t with some initial condition.
OK, so now, let me evolve the curve.
Well, OK, yeah, let me evolve a curve.
So we already took the case of a circle, which just grows,
and that's a piece that it.
Let me do a couple more example so that you see --
how does a corner evolve?
Suppose that's a t equals 0.
We have a corner.
We have sort of a fire in here, and
it's growing outward.
And we're looking for all the points at
distance 1 from this curve.
This is our curve at time 0.
OK, what do the points look like at time 1?
Well, these points here will certainly be a distance 1 away
from our curve.
These points here will be a distance 1
away from our curve.
But then there's a little circular arc, right -- maybe I
should make that dash line, too -- which are additional
points at distance 1.
So in this case, the evolving curve has smoothed out.
That corner, which we might describe as a jump in the
slope if we were thinking about jumps from the last
week's topic, has -- well, I guess I'd use the word fan.
That jump in slope has sort of the characteristics or some
key lines here, and here was an area where none of these
lines are matching, but we have distance 1 still.
Somehow, the equation is -- the slope is becoming smooth,
which it wasn't initially.
Now the opposite could also happen.
A corner could appear from a smooth start.
So that would correspond, of course, in our thinking.
We immediately think, oh, we started with something smooth,
evolved for a little while, and a shock appeared.
So a shock would correspond to the -- it's a shock in the
slope, so to speak.
So the curve itself will be continuous.
Let me take it -- well, let me just take it straight.
I mean -- oops!
Well, I was hoping to make it tangent so let me make it more
obviously tangent, so a little -- maybe an arc of a circle or
something, and then straight again, just as a model that
you can see what's happening.
OK, what will happen now, so that's again t equals 0.
So that, you see, has --
these are now --
the fire is moving normal to the surface, normal to the
fire front.
And you see what happens now is, after a certain time, say,
t equal to 1, we're going to have a corner appear.
This'll get steeper and steeper.
The slope will change faster and faster.
Here it takes a little while to change from that slope.
It curves around and slowly gets to that slope.
But that is going to get squeezed in as we move forward
in time and -- well, it's essentially time running in
the opposite direction from this box.
So this would be a t equal to 1 with a corner.
OK, and maybe those examples will explain why the idea of a
particle method isn't great because if we put particles --
I'll just say that again.
If we tried to follow this, which would be a completely
natural idea in discretizing the differential problem,
would be to follow particles there.
They travel perpendicular to the fire.
They reach here and here, but nobody reaches there, so we
don't have any information there.
And if we did try the particle method here, we'd have a bunch
of particles that start following the law of moving --
well, we have to figure out where the tangent direction
is, and perpendicular to that is the normal direction.
The particles move along, and you see that all these
particles are going to be really close together after a
little time, and then there, they meet.
So the particle method just basically fails on a
problem like this.
On a problem that stays beautiful and smooth and keeps
it shape, probably the particle -- the following
ordinary differential equation for each particle and then
interpolating to find a curve is reasonable.
But level set method is a different way.
OK, so let me stay with the level set idea, which
involves, of course, remembering if this is the
curve P of x, y equals 0, then the normal direction is the
same direction as the gradient.
I mean, that's the one bit of calculus that we absolutely
need, is that the gradient of a function points in the
normal direction.
And if we want a unit vector n, the unit normal vector, we
have to divide by its length.
So that's the basic formula from
calculus of two variables.
And then there's a special case of a function when the
gradient is 1.
So we have -- that formula simplifies because the
denominator is just a 1 for these functions.
OK, so these are distance functions.
Functions with gradient everywhere 1 are kind of
special and nice, very special and nice.
And I guess that if everything is --
so I spoke about distance function here.
So I allowed myself to take -- to imagine that the function
phi was the distance, that everything was traveling with
speed F equal 1.
So let me introduce that letter F for the speed, and in
this case, it's 1.
So add time t.
The function that gives the distance has gradient 1.
Maybe an example would serve.
Let me go back to this circle here.
OK, suppose these are my --
I'm not doing this tricky case with two circles, but just
back to a single circle.
So if I take the gradient of that function, I get -- it's a
vector, right?
The x-derivative and then the y-derivative, it's 2x and 2y.
And that's a vector with some magnitude.
Probably the magnitude is about 2r, the distance from
the center.
It's not 1.
This is not a distance function.
The distance function that has the same level curves but has
gradient equal 1 would just be the square root of x squared
plus y squared.
That's a distance.
OK, so that distance function -- say, and let me measure
distance from the center.
OK, so the distance from the origin is 0,
when x and y are 0.
I claim that this distance function is -- that that
function is this perfect kind for the level set method where
its level sets are still those same circles, but the gradient
of that thing, the gradient of this -- of course,
this is just r.
And the gradient of it, if I take the x-derivative of this
-- so let me take the gradient of the distance, just in case
you were really remembering what we easily mastered in
calculus of several variables.
The x-derivative of that is, well, I get 2x, and then I get
a half from the square root, and then the square root goes
in the denominator.
So that's just x over r, x over the square root, if you
like. x over the square root is the first component.
That's the x-derivative.
The y-derivative of this is 2y from the chain rule, and then
the whole square root produces 1/2 square root to the --
square root in the denominator.
In other words, it produces y over the square root.
And sure enough -- so that's the gradient and sure enough,
if I take its magnitude, I take this squared plus this
squared, square root, but I get 1.
OK, so distance functions are the good ones.
So in working with the level set method, we will follow
equations, numerically, of course, and phi of x and y
might start out -- phi of x and y and 0 might be a
distance function, but later on in time after we've
followed phi, it probably changes
from a distance function.
And if it changes a lot, if its gradient is getting out of
hand, for example, the gradient of
this is getting bigger.
And after awhile, we can stop, reinitialize to make the
gradient 1 again.
Effectively, we would come back to this
function and proceed.
So this reinitialization is part of the numerical
computation that I won't say more about.
Really what I now have to do is say
what's the key equation?
What's the equation that we solve numerically in x, y, t?
Remember, it's going to be a partial differential equation.
OK, so first we have to know what's the rule for the
evolution of a curve.
So this decides the movement of the curve.
This decides the movement of the curve, this velocity
vector, velocity field, velocity vector, v, v1, v2.
It's a vector.
All right, then the neat thing is that we get by using this
level set description of the curve to this equation.
The convection equation, you could say.
It's like our one-way wave equation.
If v is a constant, it's just our equation, but now we're in
two space dimensions, instead where before, we
were just in one.
But if v is not a constant, then, for example, v may
depend on phi.
The velocity might depend, and it often does depend on the
curvature of phi in some nonlinear problem, then our
equation becomes nonlinear, but it still has a nice form.
So that expresses exactly what our thinking is that
the level sets --
that the function, the evolution of the function,
tells us how the level sets evolve, and then, at any time,
we read off the position of the curve from the level set.
And then one more step to that equation, which is often
easier computationally than this one.
So I'm going to take a simple step between there and there
in the case when v is in the normal direction.
So let me look at this v dot grad phi that came into the
convection equation.
I'll divide and multiply by the magnitude of the gradient,
and then I recognize this as the normal vector.
I recognize that as n.
So this is v dot n, and I'm going to give a name to v dot
n and call it F. That's the speed, speed in the normal
direction, v dot n.
And then there's this factor magnitude of
grad phi coming along.
OK, so that's the equation.
That's sort of the famous level set equation where if F
is a constant, that's the case we're speaking
about of a fire moving.
But if F depends on the --
F could be a function of -- if I put parentheses in there, F
could be a function of gradient phi, a more
complicated function.
For example, I wrote down, because I would otherwise
forget, this could be our --
this is a possible F of grad phi.
Actually, it depends on more than the magnitude so I'd
better -- yeah.
And certainly not a linear one because of that denominator.
So a big set of examples in this subject, not linear, are
curvature shapes.
Actually, the curvature problem is quite interesting.
If the movement depends on the curvature --
so suppose I start with an elliptical-looking shape, the
curvature is smaller here, right?
Oh, I'm sorry, smaller curvature there, larger
curvature there, and if we choose the right evolution,
often it will come -- maybe this is a case where the
movement is inward because the curvature is somehow the --
they're probably vectors pointing in.
Here it would move in slowly because the
curvature is small.
Here it would move in faster because the
curvature is larger.
The curvature is tending to even out.
And, in fact, after a certain -- as time goes on, any shape
approaches a circle, and the circle shrinks and shrinks and
disappears.
So that this class of -- this sort of nonlinear group of
problems has a nice geometric picture, that whatever shape
you begin with, you could even begin with a really weird
shape like a barbell or something.
So this has zero curvature.
It wouldn't move at the start.
This would move in, so after a little while it would be like
this, then the corner would move in.
So it seems a little unlikely it's true.
After awhile, the shape would become nearly circular, get
smaller and smaller, and the shape would vanish.
Anyway, that a curvature problem.
OK, so I guess following the main theme of the course would
be write down a finite difference method for the
level set equation.
That's the natural thing to do.
Write down a finite difference method for
the level set equation.
Solve that numerically.
So that would be on the web page.
The difference equation is going to be a sort of upwind
equation, but it's a little messier in two variables to
write out in detail here.
You'll see it there on the web.
OK, so you solve it, and proceed in time, and
reinitialize when necessary and/or when it seems
advisable, and you get a great picture of this evolving,
moving interface.
So problems of that kind.
Or optimizing the interface.
There are other problems where you might -- ones on the web
where you might be look --
the unknown might be a shape.
And you're looking for the shape that maybe makes an
eigenvalue as small as possible, subject to some
constraint.
So problems for shapes.
That's really the type of problem that level set methods
are created for.
OK, and these example show the connection --
I'd say that somehow the -- yeah, it's really that if I
tried to connect it to conservation laws, to the
previous lectures, this is our earlier guy, would be, in one
variable, let's say, this would correspond to the
derivative of maybe phi or capital U. So the conservation
laws for little u correspond to equations
like this for phi.
Well, or really more rad phi in more dimensions.
OK, well, all right.
So what else remains to say about this topic?
Really I guess my main point is to introduce you to the
idea of level sets and level set method for when problems
of moving interfaces and shape optimization come up.
I guess I also wanted to say something about this fast
marching method.
OK, so what's that?
That's a case where there's a very fast
solution to this problem.
And that case occurs when the wave front is always moving
the same way, when the F holds its sign.
So it would apply here.
So I could use the fast marching method.
Rather than the level set equation and a finite
difference method for it, I could use a fast marching
method for this problem because I have no danger of
the shape turning, changing direction, coming back in on
itself, all the crazy things that can happen, the things
that are illustrated on the web.
OK, so the fast marching method, well, maybe I -- so
what's the problem in that nice case?
I want to find --
I could put on a mesh.
Yeah, let me take as a problem the computation of the
distance function from this curve.
OK, so I could create a mesh.
So I'm going to take that problem and turn it into a --
discretize it in a natural way, and I'm searching for the
distance function.
So d of x and y is the distance that
I now want to compute.
It's the distance from the curve C. And let me make an
assigned distance function so I'll make d negative in here
and positive outside of the curve.
OK.
All right, so the idea is we could -- this fits our level
set picture, but there's a direct way to do it.
And it's a nice -- it's an algorithm that computer
scientists probably created and appreciated, that we want
to know the distance.
We want to compute d at the mesh points.
Let me put another mesh point.
OK, so how do I -- so the question is, how do I get from
-- let me put another y in there, too, and that created
some more points.
So how to -- what path -- you know, how to get the distance
from the curve to each mesh point?
Well, here's the idea, just in a nutshell.
You keep a list of tentative distances to all the mesh
points, and you look for the one that's closest, like it's
probably maybe that guy.
So maybe that's right on the curve, let's say.
So the idea is I can -- that will be --
to get to this point, I'm never going to take a path
that passes through the others.
To get here, I might take a path that goes through an
earlier point.
But the nearest point, the work is done.
So you successively sweep through the mesh points, pick
the one whose current distance, tentative distance,
is the smallest, and you say that is the distance.
That's the right number.
Then now you know how far it is to get to that point.
So now you could, if you had to, update all the previous --
all the other distances by knowing the correct one, how
distant it is, how distant that point is.
OK, so then that one's settled.
Now you look at your list, your updated list.
You pick the smallest one in that, the closest one in that.
You found the best way because things are never going to
reverse in sign.
You're never going to travel out and come back, so the
closest one, you've got it.
So that might be this or something, and then we would
update how far it takes to get to those.
So that fast marching method is a, you might say, fairly
natural idea of how to compute distances when everything is
measured, when our problem is --
has a fixed sign.
OK, so that's my short comment on the fast marching method.
My longer comment is on the level set method and level set
equation where I haven't written out the finite
difference approximation that people typically use and note
on the web.
All right, so that's level sets, which I probably will
not come back to.
It's kind of a one-shot lecture.
So let me take the remaining time to get
feedback on the homework.
So were there four problems, possible candidate problems,
for homework?
Let me put those down again.
I'm just going to ask first, who did which problem.
So the first problem was the conservation law -- no, was
the one way --
OK, remind me.
I don't remember the order I gave.
So was Number 4 the conservation law?
The nonlinear one?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah, OK.
So the conservation law.
I'll just write down ut plus uux equals 0.
Oh yeah, Number 1 was the basic question
of ut equals cux.
Yeah, that's right.
So that's a very important and not too obvious question, a
good numerical question, of what to do.
And Number 2 was the wave equation; was that right?
AUDIENCE: Schrodinger's.
PROFESSOR: Schrodinger's, Number 2.
Oh, boy.
OK.
So iuxx, which I haven't discussed in
class and I had --
the good question came up, how do you deal with these complex
numbers because, of course, it's going to get complex
right away.
Well, one way that just occurs to me is take the real and
imaginary parts so that we get two
equations, two real equations.
PROFESSOR: OK, Number 3 was -- what did I have --
AUDIENCE: The wave equation.
PROFESSOR: Was the wave equation?
OK, second order?
AUDIENCE: First order.
PROFESSOR: Oh, first order wave equation?
Oh, OK.
So what was I --
OK, first order wave equation like so.
I'm sorry?
AUDIENCE: [INAUDIBLE]
PROFESSOR: cux, OK, yep.
All right, so was the question there to go back to the--
AUDIENCE: Compare.
PROFESSOR: Compare, right.
Compare the methods we've already established.
And then Number 4 was to tackle some case of nonlinear
equation where shocks and fans could appear.
OK, maybe I just get a sense of who did what.
Number 1 how many?
Wow!
Number 1 was popular and I think it's justified.
So this is certainly the biggest choice.
So let me come back and get some
ideas of what you concluded.
Number 2, did anybody tackle Schrodinger?
Yes.
A couple of complex guys.
Well, I know you asked about it.
Yep, OK, good.
Anybody go back to these?
For this?
You did?
OK, you may have done a little bit of both.
And did anybody venture to nonlinear world?
One.
Good, OK.
Nonlinear person.
I'm glad to know that, right.
So this is like -- this homework you could think of us
as like prep for a possible project.
Your project could grow out of this homework, or it could
grow from other things that are coming.
OK.
All right, let's stick with this one then since a lot of
people will have input on it.
Let me ask -- well, I'll just ask naive questions because
you've actually done the calculation.
Is there a method you'd recommend?
If you had a problem like that or, I mean, like for real, and
you had to invest in a code and use it,
just give an opinion.
Or did you only try one method?
Tell me something, positive or negative, about your
experience with this.
Yeah, thanks.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Maybe you push -- in theory, there's some mike
you push and this will be the first time we've ever done it.
OK, terrific.
Go ahead.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Centered explicit.
Let me write that down.
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK, right.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes.
AUDIENCE: [INAUDIBLE]
PROFESSOR: That's right.
So, OK, the centered explicit one is computing this from
three old values.
We certainly need three because we've got a second
derivative there.
But then we have a certain amount of that second
derivative and a certain amount of this one, and that
can make one of those three numbers go negative.
They always add the one because the constant is
certainly a solution we want to keep.
OK, so when a coefficient went negative, that doesn't mean
instability; is that right?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah.
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Fine.
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Implicit had no trouble with that.
Yeah, OK.
Or it didn't matter what the Peclet number, or the cell
Peclet number was, or whatever we called it.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Uh-huh.
AUDIENCE: [INAUDIBLE]
PROFESSOR: I see.
AUDIENCE: [INAUDIBLE]
PROFESSOR: OK, if you went --
when you say implicit, it's just the diffusion term that's
going implicit?
OK, OK.
So that should be a very stable process, but, of
course, now we've got this term affecting it.
So you actually found out that even the
implicit one could break?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yes.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Really?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Huh!
AUDIENCE: [INAUDIBLE]
PROFESSOR: Even the implicit one?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah, yeah.
So on the implicit part -- so this was method one and this
was method two implicit for the diffusion, but you found
that that could still be unstable.
OK, yeah.
Anybody else do this comparison?
Can you push your mike button just to --
this is a hotshot room that we haven't taken advantage of,
but it's our chance.
AUDIENCE: I actually looked at the -- when I first did the --
looking at just the growth factors of explicit and
implicit --
PROFESSOR: OK, looking at that g --
AUDIENCE: --on the convection term as well?
PROFESSOR: Yeah.
AUDIENCE: And I used implicit for all the diffusion.
And I found the most stable using upwind, whereas when
both were implicit, the growth factor was actually able to
stay below 1.
PROFESSOR: So you moved the first derivative also up to
the new time?
OK, and that was -- then you don't have a problem.
AUDIENCE: Yeah, exactly.
So stable for pretty much all big R, little r.
But the explicit upwind for the convection term was stable
for most as long as the little r was below 0.5.
PROFESSOR: I see.
So the numbers could tell you the--
AUDIENCE: Sure.
PROFESSOR: --limits on little r and big R?
AUDIENCE: Yep.
And then I did the explicit convection upwind with the
implicit center diffusion, and that seemed to be pretty
stable unless I really started cranking the numbers out.
PROFESSOR: I see.
Oh, that's the one that he had also done?
AUDIENCE: Yep.
PROFESSOR: Yeah, OK.
And now so I guess I remember in the figure that isn't yet
drawn that might emerge from your calculations a
possibility where -- depending on the Peclet number, where it
was stable.
It didn't blow up, but oscillations appeared and the
calculations -- you wouldn't do it even so --
AUDIENCE: [INAUDIBLE]
PROFESSOR: That's what you meant, yeah.
OK, so there are cases where the method isn't good even
though it's stable, which is like important to realize that
that can happen.
Because up to now, stability has been, you
know, the holy grail.
But it's really a good solution
that's the holy grail.
Yes?
Go ahead.
AUDIENCE: I did kind of a bigger picture where I
basically tried to build the pieces so just to isolate the
advection term and the diffusion term and then build
their matrices.
PROFESSOR: Yes.
AUDIENCE: And then look at the combination rule.
So when you went to where you had both terms present--
PROFESSOR: Right.
AUDIENCE: --how the matrices were actually -- came together
to form the--
PROFESSOR: So you took a matrix picture of the --
with all the mesh points, mesh values, at once rather than--
AUDIENCE: And then I simulated with them just convection,
just diffusion and then the coupled problem of both of
them to see how that -- just one method may be unstable for
a given formulation--
PROFESSOR: Right.
AUDIENCE: --but yet when you saw the coupled problem with
those two pieces, it actually is stable.
So you could look at like the -- when you build the
matrices, you can look at the eignevalues and you can make
sure that they're--
PROFESSOR: Let me comment a little on that.
I'll do that, but I don't want to cut off anybody who has
another comment to throw in here before I say something
about -- follow up on matrices.
Anybody got some urgent thing?
Let me just, in the remaining minute, just say a word
because, of course, you'll know that I think of a matrix
description quickly and fondly.
So just to write down the sort of matrices that are in here,
from the second difference, the second difference matrix
is going to be -- it's going to have our --
well, I guess it'll be 1 minus 2 1 1 in a typical row.
We haven't written any matrices down.
It's amazing that --
we'd be at about the, I don't know, eighth or ninth lecture
now, and not have seen this matrix, but here it is.
OK, so that's the matrix where I'm at a fixed time level.
I'm taking the second differences of all the points
along that level.
And now the main point for the moment is to contrast that
with a first difference matrix which, well, that depends
whether I take it centered.
Maybe if it's a centered difference, it would
be a 1 0 minus 1.
And, of course, there's a 2 delta x in the denominator.
so a minus 1 0 1.
So that would be a centered difference that doesn't do
anything, doesn't use the middle value, but takes the
difference between the value at the right and the value at
the left, or an upwind.
So delta upwind would be more of a 1's forward of the point
and minus 1's at ujn and 0's here.
OK.
So what you're saying, am I right in thinking you took
some combination of the--
AUDIENCE: An implicit method for one and explicit for the
other, you can do them independently by you can
invert the implicit matrix--
PROFESSOR: Right.
AUDIENCE: --and premultiply the other side and build.
PROFESSOR: So this is a way to see, and actually the next
months of our course are going to be about, how do you deal
with these large difference matrices.
These are in one dimension so they're tridiagonal, and
they're, of course, easy to deal with it.
But problems, serious problems,
come up in more dimension.
But here, can I just express one thought here?
The eigenvalues of these matrices and their
combinations are what's going to control, as eigenvalues
always do, growth of solution.
Now we already had numbers g, growth factors, that we
computed by the von Neumann idea of following
exponentials, e to the ikx.
So these numbers too were -- their purpose was to track the
growth of solution.
And what I want to say is that very often, I mean, that these
numbers and the eigenvalues of the matrices, they're very
closely related ideas.
In fact, they'll be the same idea if the eigenvector is a
pure exponential.
So that's really the link between eigenvalues of
matrices and growth factors g that depended on frequency k.
If the eigenvector was a pure e to the ikx on the mesh, then
that growth factor is identical to the eigenvalue
that we get it one step of the difference equation.
So the idea was that the eigenvectors, the details of
the eigenvectors, will be affected by boundary
conditions.
So von Neumann thought, all right, I'll cut through that
Gordian knot and just pretend I have a pure e to the ikx.
I look inside the region, not at the boundary where we may
be doing something different, and we get a g, a
growth factor g.
And again, if the eigenfunction was really an
eigenvector, then that g would be the same as lambda.
OK, we're past time.
I always appreciate that you are patient about that.
So I'll collect those that are complete.
And if you wanted to speak to Mr. Cho or have a little more
time to do more, Monday is good, too.
OK, thanks.