Lecture 11: Level Set Method

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

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.

Free Downloads

Video


Audio


Caption

  • English-US (SRT)