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 pf 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 a 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 xy-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 xy-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 non-circle, 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 meet, 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 of it.
Let me do a couple more example so that you see --
how does a corner evolve?
Suppose that's at 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 dashed 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 at 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 phi 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 at 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 a while, 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, (v_1, v_2).
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 appropriate 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 that after a while,
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 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 -- actually,
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, grad 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 it a signed 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 appreciate, 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 u_t plus u*u_x equals 0.
Oh yeah, Number 1 was the basic question of u_t equals c*u_x.
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 i*u_xx, 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: c*u_x, 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,
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 centered 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 eigenvalues 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 the 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, and 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 u_(j, n) 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 pre-multiply 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, as you say,
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 i*k*x.
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 i*k*x on the mesh,
then that growth factor is identical to the eigenvalue
that we get in 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 i*k*x.
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.