1 00:00:00,499 --> 00:00:01,992 The following content is provided 2 00:00:01,992 --> 00:00:06,060 by MIT OpenCourseWare under a Creative Commons license. 3 00:00:06,060 --> 00:00:08,110 Additional information about our license 4 00:00:08,110 --> 00:00:10,710 and MIT OpenCourseWare in general 5 00:00:10,710 --> 00:00:11,930 is available at ocw.mit.edu. 6 00:00:15,470 --> 00:00:18,900 PROFESSOR: OK, so I'll start with this. 7 00:00:18,900 --> 00:00:25,150 This is a method for a particular class of problem 8 00:00:25,150 --> 00:00:30,490 and maybe moving interfaces, moving curves, 9 00:00:30,490 --> 00:00:34,010 curves that are -- or surfaces, in higher dimensions, 10 00:00:34,010 --> 00:00:40,910 that are being transported, and it's their shape that matters. 11 00:00:40,910 --> 00:00:48,010 So it's a sort of topic within shape differential equations 12 00:00:48,010 --> 00:00:51,970 rather than particle differential equations. 13 00:00:51,970 --> 00:00:57,600 And some cases of this problem can 14 00:00:57,600 --> 00:01:00,500 be studied by a particularly fast method, 15 00:01:00,500 --> 00:01:02,490 developed here at MIT in Course 6 16 00:01:02,490 --> 00:01:06,730 actually, called the fast marching method. 17 00:01:06,730 --> 00:01:09,940 It's really thanks to Per-Olof Persson 18 00:01:09,940 --> 00:01:15,400 that there's now quite a bit on the 18.086 web page. 19 00:01:15,400 --> 00:01:21,770 The basic equations that'll be in the lecture and also 20 00:01:21,770 --> 00:01:26,450 demonstrations, quick demos that show the point. 21 00:01:26,450 --> 00:01:34,790 And I better mention Sethian and Osher, who jointly proposed 22 00:01:34,790 --> 00:01:41,760 the method, conceived of this approach in the late '80s, 23 00:01:41,760 --> 00:01:44,990 and each has written a book, or coauthored a book, 24 00:01:44,990 --> 00:01:47,410 on the level set method. 25 00:01:47,410 --> 00:01:52,320 So those books are, so to speak, the authorities on the method. 26 00:01:52,320 --> 00:01:56,040 OK, so what's the idea? 27 00:01:56,040 --> 00:02:04,040 Well, here is the curve, the shape, 28 00:02:04,040 --> 00:02:05,760 the interface that's moving. 29 00:02:05,760 --> 00:02:11,800 And let's suppose -- well, imagine it's sort of -- 30 00:02:11,800 --> 00:02:17,660 there's a wall pf fire here, and think of it as a moving fire 31 00:02:17,660 --> 00:02:18,970 front. 32 00:02:18,970 --> 00:02:26,100 Then that fire front will move, normal to itself, 33 00:02:26,100 --> 00:02:30,780 to some later position at time 1. 34 00:02:30,780 --> 00:02:38,480 We could normalize -- in a perfect case when it was 35 00:02:38,480 --> 00:02:43,050 homogeneous and the speed of movement didn't change with 36 00:02:43,050 --> 00:02:46,550 position, we could imagine that the speed is 1. 37 00:02:49,350 --> 00:02:53,020 So where is this fire front here at time 1? 38 00:02:53,020 --> 00:02:57,050 It contains all the points that are 39 00:02:57,050 --> 00:03:01,710 reached from the initial fire in unit time. 40 00:03:01,710 --> 00:03:04,110 And since we're going at speed 1, 41 00:03:04,110 --> 00:03:11,820 let's say, this is all the points at distance 1. 42 00:03:11,820 --> 00:03:16,210 This curve, in this best, simplest case, 43 00:03:16,210 --> 00:03:21,020 will be the curve that's everywhere 44 00:03:21,020 --> 00:03:23,050 a distance 1 away from this curve. 45 00:03:25,650 --> 00:03:30,350 Now, how to describe the curves is the key point here. 46 00:03:30,350 --> 00:03:35,910 Shall I describe the curve by particles? 47 00:03:35,910 --> 00:03:39,800 I could imagine these points as being 48 00:03:39,800 --> 00:03:46,470 my unknowns, my positions, and some equation for following 49 00:03:46,470 --> 00:03:47,340 them. 50 00:03:47,340 --> 00:03:50,410 So I follow those -- ordinary differential equations, 51 00:03:50,410 --> 00:03:58,760 one for every point -- and they reach new positions at time 1. 52 00:03:58,760 --> 00:04:05,970 And then I interpolate the next curve, the new curve at time 1. 53 00:04:05,970 --> 00:04:08,090 Well, that's not what we do. 54 00:04:08,090 --> 00:04:10,340 That's not the idea of the level set method, 55 00:04:10,340 --> 00:04:16,280 and I'll explain why that idea of following particles, 56 00:04:16,280 --> 00:04:21,760 which is a totally natural one, comes to grief, numerically, 57 00:04:21,760 --> 00:04:22,340 too often. 58 00:04:26,470 --> 00:04:31,750 This picture, of course, was like an ideal case, 59 00:04:31,750 --> 00:04:37,270 and the particles separated a little bit, but not too much. 60 00:04:37,270 --> 00:04:40,290 Just to say in a word, the trouble with the particle idea 61 00:04:40,290 --> 00:04:43,240 is that particles could separate wildly, 62 00:04:43,240 --> 00:04:48,570 so you had a big, big gap to fill in, like a fan. 63 00:04:48,570 --> 00:04:52,600 Or they could come very close together so that you had 64 00:04:52,600 --> 00:04:55,040 total instability numerically. 65 00:04:55,040 --> 00:04:56,850 You couldn't compute any differences 66 00:04:56,850 --> 00:04:59,170 because the distance was too small, 67 00:04:59,170 --> 00:05:01,970 and that would be like a shock. 68 00:05:01,970 --> 00:05:05,070 So we're really connecting this type of problem. 69 00:05:05,070 --> 00:05:08,070 It connects well with the conservation laws 70 00:05:08,070 --> 00:05:09,540 that we've been studying. 71 00:05:09,540 --> 00:05:12,650 Fans and shocks appear, and I'll show you right away. 72 00:05:12,650 --> 00:05:16,500 Can I just say, what is the approach? 73 00:05:16,500 --> 00:05:19,080 What does this word "level set" mean? 74 00:05:19,080 --> 00:05:23,530 Well, the level set of a function -- 75 00:05:23,530 --> 00:05:25,980 I'm saying what this word level set is, 76 00:05:25,980 --> 00:05:29,520 and that's going to be the way the curve is described, 77 00:05:29,520 --> 00:05:30,810 sort of implicitly. 78 00:05:30,810 --> 00:05:33,340 Instead of saying where the points are 79 00:05:33,340 --> 00:05:35,580 that are on the curve, the curve is 80 00:05:35,580 --> 00:05:42,570 being given by an equation like phi of x, y equals 0, 81 00:05:42,570 --> 00:05:46,470 or phi of x, y equal 1, or phi of x, y equals 7. 82 00:05:46,470 --> 00:05:53,870 Those three curves would be level sets of phi. 83 00:05:53,870 --> 00:05:56,270 Do you see the idea of levels sets? 84 00:05:56,270 --> 00:05:59,760 They're the solutions to the equation 85 00:05:59,760 --> 00:06:01,480 phi of x, y equal constant. 86 00:06:01,480 --> 00:06:08,220 Let me take a simple example here on the board that I can -- 87 00:06:08,220 --> 00:06:13,370 so the idea of level sets, first of all. 88 00:06:13,370 --> 00:06:16,840 Let me take an easy -- a particularly nice function. 89 00:06:16,840 --> 00:06:22,200 phi of x, y equals, say, x squared plus y squared. 90 00:06:22,200 --> 00:06:27,660 OK, so I want to ask, what are the level 91 00:06:27,660 --> 00:06:29,980 sets of that function? 92 00:06:29,980 --> 00:06:36,550 So I set phi of x, y equals a constant, to get the c level 93 00:06:36,550 --> 00:06:38,360 set, you could say. 94 00:06:38,360 --> 00:06:39,920 And, of course, we recognize that 95 00:06:39,920 --> 00:06:43,580 as a circle when c is positive, and there 96 00:06:43,580 --> 00:06:45,970 is no set when c is negative. 97 00:06:45,970 --> 00:06:50,360 So this would be -- the level sets would be circles. 98 00:06:59,140 --> 00:07:00,590 So that would be one description. 99 00:07:00,590 --> 00:07:02,690 I mean, what we're really doing is 100 00:07:02,690 --> 00:07:09,930 we're introducing a new dimension to the problem. 101 00:07:09,930 --> 00:07:13,940 Here, it's a 2D problem with x square and y squared, 102 00:07:13,940 --> 00:07:18,860 but I'm bringing in the third dimension, z, 103 00:07:18,860 --> 00:07:25,310 and thinking of the curve in the xy-plane as being produced 104 00:07:25,310 --> 00:07:29,290 by -- so what would be the graph of z equals x squared plus y 105 00:07:29,290 --> 00:07:29,790 squared? 106 00:07:29,790 --> 00:07:32,530 That's some surface now in 3D. 107 00:07:32,530 --> 00:07:40,390 So I have x, y, and now I've introduced this third dimension 108 00:07:40,390 --> 00:07:43,880 in order to see a whole surface. 109 00:07:43,880 --> 00:07:49,160 And now I cut that surface by the plane z equal to c. 110 00:07:49,160 --> 00:07:56,280 So z equal to c cuts through that surface by a plane here. 111 00:07:56,280 --> 00:08:00,030 And that cuts out a circle, which 112 00:08:00,030 --> 00:08:05,180 I can drop down to the xy-plane if I want to see it down there, 113 00:08:05,180 --> 00:08:08,500 and that's the level set at height c. 114 00:08:11,120 --> 00:08:18,690 So that's a way to describe the curve, but it's sort of -- 115 00:08:18,690 --> 00:08:20,950 it's implicit. 116 00:08:20,950 --> 00:08:24,780 You'll see that it has major advantages when -- 117 00:08:24,780 --> 00:08:28,250 maybe I can just describe them without drawing them. 118 00:08:28,250 --> 00:08:36,820 Suppose I have a problem with two fires, 119 00:08:36,820 --> 00:08:40,560 say, two circles of fire, spreading. 120 00:08:40,560 --> 00:08:43,380 What's going to happen, as time evolves, 121 00:08:43,380 --> 00:08:45,470 with these two circles of fire? 122 00:08:45,470 --> 00:08:50,490 Well, having taken circles, they're just going to expand, 123 00:08:50,490 --> 00:08:55,350 so after a little time, I've got two bigger circles. 124 00:08:55,350 --> 00:08:58,790 But after a little more time, the circles 125 00:08:58,790 --> 00:09:03,110 will meet and form one non-circle, right? 126 00:09:03,110 --> 00:09:09,360 A sort of figure-eight shape. 127 00:09:09,360 --> 00:09:13,260 Now that's a change in topology, and the world 128 00:09:13,260 --> 00:09:17,420 of shape calculations is always having 129 00:09:17,420 --> 00:09:20,450 to ask, how am I going to deal with changes in topology 130 00:09:20,450 --> 00:09:25,340 when curves meet, when unconnected things become 131 00:09:25,340 --> 00:09:28,740 connected, or when connected things disconnect? 132 00:09:28,740 --> 00:09:32,495 And the one great value of the level set approach 133 00:09:32,495 --> 00:09:35,860 is it deals with that easily. 134 00:09:35,860 --> 00:09:39,930 I could imagine a double -- so I don't know that I'll succeed 135 00:09:39,930 --> 00:09:47,500 to do it here, but I could imagine a function, which -- 136 00:09:47,500 --> 00:09:52,100 OK, maybe this example, this picture, isn't too bad. 137 00:09:52,100 --> 00:09:54,620 So that's intended to be a surface. 138 00:09:57,250 --> 00:09:59,950 Well, let me make the surface go up that way. 139 00:09:59,950 --> 00:10:03,210 OK, so when I cut it at z equals c, 140 00:10:03,210 --> 00:10:08,320 that might cut the surface in two different planes. 141 00:10:08,320 --> 00:10:14,190 So the level set at a certain level is two circles. 142 00:10:14,190 --> 00:10:15,740 There's two circles. 143 00:10:15,740 --> 00:10:19,550 But now, let z increase, and what happens? 144 00:10:19,550 --> 00:10:23,280 You see how naturally, as z increases, 145 00:10:23,280 --> 00:10:25,550 these two circles are coming closer and closer, 146 00:10:25,550 --> 00:10:28,980 and at a certain point, they merge. 147 00:10:28,980 --> 00:10:35,200 And phi of x and y remains a totally nice function. 148 00:10:35,200 --> 00:10:40,030 I mean, I could create some probably fourth-degree function 149 00:10:40,030 --> 00:10:41,500 curve that would do that fine. 150 00:10:41,500 --> 00:10:42,730 I could graph it. 151 00:10:42,730 --> 00:10:44,500 I could visualize it. 152 00:10:44,500 --> 00:10:47,430 And the point is that by going up that one dimension 153 00:10:47,430 --> 00:10:50,580 and then slicing at different levels, 154 00:10:50,580 --> 00:10:56,070 you allow these things to happen smoothly. 155 00:10:56,070 --> 00:11:00,300 OK, so that's one feature of level sets, that description. 156 00:11:00,300 --> 00:11:06,170 OK, so that's the way curves are going to be defined. 157 00:11:06,170 --> 00:11:09,280 So our unknown becomes this function, 158 00:11:09,280 --> 00:11:12,920 which will become a function of x, y, and t 159 00:11:12,920 --> 00:11:15,840 with some initial condition. 160 00:11:15,840 --> 00:11:24,390 OK, so now, let me evolve the curve. 161 00:11:24,390 --> 00:11:29,010 Well, OK, yeah, let me evolve a curve. 162 00:11:29,010 --> 00:11:35,470 So we already took the case of a circle, which just grows, 163 00:11:35,470 --> 00:11:38,040 and that's a piece of it. 164 00:11:38,040 --> 00:11:42,780 Let me do a couple more example so that you see -- 165 00:11:42,780 --> 00:11:47,430 how does a corner evolve? 166 00:11:47,430 --> 00:11:50,250 Suppose that's at t equals 0, we have a corner. 167 00:11:52,800 --> 00:11:57,980 We have sort of a fire in here, and it's growing outward. 168 00:11:57,980 --> 00:12:02,010 And we're looking for all the points at distance 1 169 00:12:02,010 --> 00:12:02,980 from this curve. 170 00:12:02,980 --> 00:12:05,370 This is our curve at time 0. 171 00:12:05,370 --> 00:12:09,370 OK, what do the points look like at time 1? 172 00:12:09,370 --> 00:12:14,210 Well, these points here will certainly be a distance 173 00:12:14,210 --> 00:12:16,050 1 away from our curve. 174 00:12:16,050 --> 00:12:19,930 These points here will be a distance 1 away from our curve. 175 00:12:19,930 --> 00:12:24,820 But then there's a little circular arc, right -- 176 00:12:24,820 --> 00:12:27,780 maybe I should make that dashed line, too -- 177 00:12:27,780 --> 00:12:34,770 which are additional points at distance 1. 178 00:12:34,770 --> 00:12:41,490 So in this case, the evolving curve has smoothed out. 179 00:12:41,490 --> 00:12:48,110 That corner, which we might describe as a jump in the slope 180 00:12:48,110 --> 00:12:54,800 if we were thinking about jumps from the last week's topic, 181 00:12:54,800 --> 00:12:58,770 has -- well, I guess I'd use the word fan. 182 00:12:58,770 --> 00:13:05,280 That jump in slope has -- sort of the characteristics or some 183 00:13:05,280 --> 00:13:12,520 key lines here, and here was an area where none of these lines 184 00:13:12,520 --> 00:13:16,100 are matching, but we have distance 1 still. 185 00:13:18,970 --> 00:13:27,630 Somehow, the equation is -- the slope is becoming smooth, 186 00:13:27,630 --> 00:13:29,830 which it wasn't initially. 187 00:13:29,830 --> 00:13:34,610 Now the opposite could also happen. 188 00:13:34,610 --> 00:13:41,360 A corner could appear from a smooth start. 189 00:13:41,360 --> 00:13:45,000 So that would correspond, of course, in our thinking. 190 00:13:45,000 --> 00:13:47,780 We immediately think, oh, we started 191 00:13:47,780 --> 00:13:51,460 with something smooth, evolved for a little while, 192 00:13:51,460 --> 00:13:52,930 and a shock appeared. 193 00:13:52,930 --> 00:13:55,740 So a shock would correspond to the -- 194 00:13:55,740 --> 00:13:58,450 it's a shock in the slope, so to speak. 195 00:13:58,450 --> 00:14:01,200 So the curve itself will be continuous. 196 00:14:01,200 --> 00:14:04,020 Let me take it -- well, let me just take it straight. 197 00:14:04,020 --> 00:14:05,860 I mean -- oops! 198 00:14:05,860 --> 00:14:08,000 Well, I was hoping to make it tangent, 199 00:14:08,000 --> 00:14:11,580 so let me make it more obviously tangent, so a little -- 200 00:14:11,580 --> 00:14:13,670 maybe an arc of a circle or something, 201 00:14:13,670 --> 00:14:17,890 and then straight again, just as a model that you can see 202 00:14:17,890 --> 00:14:19,400 what's happening. 203 00:14:19,400 --> 00:14:22,080 OK, what will happen now? 204 00:14:22,080 --> 00:14:23,920 So that's again t equals 0. 205 00:14:26,650 --> 00:14:36,320 So that, you see, has -- these are now -- 206 00:14:36,320 --> 00:14:42,700 the fire is moving normal to the surface, 207 00:14:42,700 --> 00:14:44,690 normal to the fire front. 208 00:14:44,690 --> 00:14:49,740 And you see what happens now is, after a certain time, 209 00:14:49,740 --> 00:14:55,980 say t equal to 1, we're going to have a corner appear. 210 00:14:55,980 --> 00:15:03,820 This'll get steeper and steeper. 211 00:15:03,820 --> 00:15:05,970 The slope will change faster and faster. 212 00:15:05,970 --> 00:15:08,810 Here it takes a little while to change from that slope. 213 00:15:08,810 --> 00:15:12,130 It curves around and slowly gets to that slope. 214 00:15:12,130 --> 00:15:16,070 But that is going to get squeezed in as we move forward 215 00:15:16,070 --> 00:15:21,530 in time and -- well, it's essentially time running 216 00:15:21,530 --> 00:15:24,780 in the opposite direction from this box. 217 00:15:24,780 --> 00:15:28,550 So this would be at t equal to 1, with a corner. 218 00:15:31,140 --> 00:15:39,840 OK, and maybe those examples will explain why the idea 219 00:15:39,840 --> 00:15:46,400 of a particle method isn't great because if we put particles -- 220 00:15:46,400 --> 00:15:47,760 I'll just say that again. 221 00:15:47,760 --> 00:15:51,430 If we tried to follow this, which 222 00:15:51,430 --> 00:15:55,200 would be a completely natural idea in discretizing 223 00:15:55,200 --> 00:16:00,180 the differential problem, would be to follow particles there. 224 00:16:00,180 --> 00:16:06,000 They travel perpendicular to the fire. 225 00:16:06,000 --> 00:16:10,070 They reach here and here, but nobody reaches there, 226 00:16:10,070 --> 00:16:13,520 so we don't have any information there. 227 00:16:13,520 --> 00:16:16,920 And if we did try the particle method here, 228 00:16:16,920 --> 00:16:24,080 we'd have a bunch of particles that start following the law 229 00:16:24,080 --> 00:16:28,590 of moving -- well, we have to figure out where the tangent 230 00:16:28,590 --> 00:16:31,701 direction is, and perpendicular to that is the normal 231 00:16:31,701 --> 00:16:32,200 direction. 232 00:16:32,200 --> 00:16:33,970 The particles move along, and you 233 00:16:33,970 --> 00:16:35,595 see that all these particles are going 234 00:16:35,595 --> 00:16:40,850 to be really close together after a little time, 235 00:16:40,850 --> 00:16:45,460 and then there, they meet. 236 00:16:45,460 --> 00:16:49,340 So the particle method just basically fails 237 00:16:49,340 --> 00:16:51,180 on a problem like this. 238 00:16:51,180 --> 00:16:55,570 On a problem that stays beautiful and smooth and keeps 239 00:16:55,570 --> 00:16:58,460 it shape, probably the particle -- 240 00:16:58,460 --> 00:17:01,670 the following ordinary differential equation for each 241 00:17:01,670 --> 00:17:05,430 particle and then interpolating to find a curve is reasonable. 242 00:17:05,430 --> 00:17:09,790 But level set method is a different way. 243 00:17:09,790 --> 00:17:17,510 OK, so let me stay with the level set idea, which involves, 244 00:17:17,510 --> 00:17:25,990 of course, remembering if this is the curve phi of x, y 245 00:17:25,990 --> 00:17:29,990 equals 0, then the normal direction 246 00:17:29,990 --> 00:17:32,510 is the same direction as the gradient. 247 00:17:32,510 --> 00:17:34,680 I mean, that's the one bit of calculus 248 00:17:34,680 --> 00:17:38,370 that we absolutely need, is that the gradient of a function 249 00:17:38,370 --> 00:17:42,020 points in the normal direction. 250 00:17:42,020 --> 00:17:45,920 And if we want a unit vector n, the unit normal vector, 251 00:17:45,920 --> 00:17:49,280 we have to divide by its length. 252 00:17:49,280 --> 00:17:57,930 So that's the basic formula from calculus of two variables. 253 00:17:57,930 --> 00:18:01,980 And then there's a special case of a function 254 00:18:01,980 --> 00:18:04,840 when the gradient is 1. 255 00:18:04,840 --> 00:18:07,950 So we have -- that formula simplifies 256 00:18:07,950 --> 00:18:11,310 because the denominator is just a 1 for these functions. 257 00:18:11,310 --> 00:18:14,190 OK, so these are distance functions. 258 00:18:14,190 --> 00:18:16,690 Functions with gradient everywhere 1 259 00:18:16,690 --> 00:18:21,050 are kind of special and nice, very special and nice. 260 00:18:21,050 --> 00:18:24,860 And I guess that if everything is -- 261 00:18:24,860 --> 00:18:26,960 so I spoke about distance function here. 262 00:18:26,960 --> 00:18:31,880 So I allowed myself to take -- to imagine that the function 263 00:18:31,880 --> 00:18:35,280 phi was the distance, that everything was traveling with 264 00:18:35,280 --> 00:18:37,330 speed F equal 1. 265 00:18:37,330 --> 00:18:43,250 So let me introduce that letter F for the speed, 266 00:18:43,250 --> 00:18:45,420 and in this case, it's 1. 267 00:18:45,420 --> 00:18:51,620 So at time t -- the function that gives the distance has 268 00:18:51,620 --> 00:18:52,310 gradient 1. 269 00:18:52,310 --> 00:18:54,570 Maybe an example would serve. 270 00:18:54,570 --> 00:18:59,250 Let me go back to this circle here. 271 00:18:59,250 --> 00:19:06,630 OK, suppose these are my -- I'm not doing this tricky case with 272 00:19:06,630 --> 00:19:09,160 two circles, but just back to a single circle. 273 00:19:11,690 --> 00:19:18,000 So if I take the gradient of that function, I get -- 274 00:19:18,000 --> 00:19:19,380 it's a vector, right? 275 00:19:19,380 --> 00:19:24,910 The x-derivative and then the y-derivative, it's 2x and 2y. 276 00:19:24,910 --> 00:19:28,980 And that's a vector with some magnitude. 277 00:19:28,980 --> 00:19:30,890 Probably the magnitude is about 2r, 278 00:19:30,890 --> 00:19:33,310 the distance from the center. 279 00:19:33,310 --> 00:19:34,170 It's not 1. 280 00:19:34,170 --> 00:19:37,320 This is not a distance function. 281 00:19:37,320 --> 00:19:41,510 The distance function that has the same level curves 282 00:19:41,510 --> 00:19:44,260 but has gradient equal 1 would just 283 00:19:44,260 --> 00:19:48,850 be the square root of x squared plus y squared. 284 00:19:48,850 --> 00:19:51,430 That's a distance. 285 00:19:51,430 --> 00:19:54,930 OK, so that distance function -- say, 286 00:19:54,930 --> 00:19:58,070 and let me measure distance from the center. 287 00:19:58,070 --> 00:20:01,390 OK, so the distance from the origin 288 00:20:01,390 --> 00:20:04,480 is 0, when x and y are 0. 289 00:20:04,480 --> 00:20:08,840 I claim that this distance function is -- 290 00:20:08,840 --> 00:20:14,080 that that function is this perfect kind for the level set 291 00:20:14,080 --> 00:20:18,400 method where its level sets are still those same circles, 292 00:20:18,400 --> 00:20:23,140 but the gradient of that thing, the gradient of this -- 293 00:20:23,140 --> 00:20:25,860 of course, this is just r. 294 00:20:25,860 --> 00:20:31,545 And the gradient of it, if I take the x-derivative of this 295 00:20:31,545 --> 00:20:35,970 -- so let me take the gradient of the distance, 296 00:20:35,970 --> 00:20:41,940 just in case you were really remembering what we easily 297 00:20:41,940 --> 00:20:48,340 mastered in calculus of several variables. 298 00:20:48,340 --> 00:20:53,740 The x-derivative of that is, well, I get 2x, 299 00:20:53,740 --> 00:20:55,770 and then I get a half from the square root, 300 00:20:55,770 --> 00:20:57,970 and then the square root goes in the denominator. 301 00:20:57,970 --> 00:21:03,990 So that's just x over r, x over the square root, if you like. 302 00:21:03,990 --> 00:21:07,980 x over the square root is the first component. 303 00:21:07,980 --> 00:21:09,760 That's the x-derivative. 304 00:21:09,760 --> 00:21:15,980 The y-derivative of this is 2y from the chain rule, 305 00:21:15,980 --> 00:21:21,080 and then the whole square root produces 1/2 square root 306 00:21:21,080 --> 00:21:23,630 to the -- square root in the denominator. 307 00:21:23,630 --> 00:21:28,400 In other words, it produces y over the square root. 308 00:21:28,400 --> 00:21:32,520 And sure enough -- so that's the gradient and sure enough, 309 00:21:32,520 --> 00:21:39,090 if I take its magnitude, I take this squared plus this squared, 310 00:21:39,090 --> 00:21:42,660 square root, but I get 1. 311 00:21:42,660 --> 00:21:47,420 OK, so distance functions are the good ones. 312 00:21:47,420 --> 00:21:51,190 So in working with the level set method, 313 00:21:51,190 --> 00:21:59,110 we will follow equations, numerically, of course, 314 00:21:59,110 --> 00:22:03,050 and phi of x and y might start out -- 315 00:22:03,050 --> 00:22:08,820 phi of x and y and 0 might be a distance function, 316 00:22:08,820 --> 00:22:13,670 but later on in time after we've followed phi, 317 00:22:13,670 --> 00:22:15,900 it probably changes from a distance function. 318 00:22:15,900 --> 00:22:21,330 And if it changes a lot, if its gradient is getting out of hand 319 00:22:21,330 --> 00:22:27,220 -- for example, the gradient of this is getting bigger. 320 00:22:29,780 --> 00:22:36,980 And after a while, we can stop, reinitialize to make 321 00:22:36,980 --> 00:22:38,790 the gradient 1 again. 322 00:22:38,790 --> 00:22:43,210 Effectively, we would come back to this function and proceed. 323 00:22:45,810 --> 00:22:51,220 So this reinitialization is part of the numerical computation 324 00:22:51,220 --> 00:22:53,060 that I won't say more about. 325 00:22:53,060 --> 00:23:00,170 Really what I now have to do is say what's the key equation? 326 00:23:00,170 --> 00:23:04,000 What's the equation that we solve numerically in x, y, t? 327 00:23:04,000 --> 00:23:11,050 Remember, it's going to be a partial differential equation. 328 00:23:11,050 --> 00:23:14,040 OK, so first we have to know what's 329 00:23:14,040 --> 00:23:18,360 the rule for the evolution of a curve. 330 00:23:18,360 --> 00:23:22,710 So this decides the movement of the curve. 331 00:23:22,710 --> 00:23:30,640 This decides the movement of the curve, this velocity 332 00:23:30,640 --> 00:23:36,877 vector, velocity field, velocity vector, v, (v_1, v_2). 333 00:23:36,877 --> 00:23:37,460 It's a vector. 334 00:23:41,330 --> 00:23:44,680 All right, then the neat thing is 335 00:23:44,680 --> 00:23:51,810 that we get, by using this level set description of the curve, 336 00:23:51,810 --> 00:23:54,750 to this equation. 337 00:23:54,750 --> 00:23:56,470 The convection equation, you could say. 338 00:23:56,470 --> 00:23:59,230 It's like our one-way wave equation. 339 00:23:59,230 --> 00:24:03,950 If v is a constant, it's just our equation, 340 00:24:03,950 --> 00:24:07,740 but now we're in two space dimensions 341 00:24:07,740 --> 00:24:10,760 instead, where before, we were just in one. 342 00:24:10,760 --> 00:24:15,990 But if v is not a constant, then, for example, 343 00:24:15,990 --> 00:24:19,860 v may depend on phi. 344 00:24:19,860 --> 00:24:22,230 The velocity might depend, and it often 345 00:24:22,230 --> 00:24:29,010 does depend on the curvature of phi in some nonlinear problem, 346 00:24:29,010 --> 00:24:32,000 then our equation becomes nonlinear, 347 00:24:32,000 --> 00:24:34,880 but it still has a nice form. 348 00:24:34,880 --> 00:24:40,280 So that expresses exactly what our thinking is that the level 349 00:24:40,280 --> 00:24:45,440 sets -- that the function, the evolution of the function, 350 00:24:45,440 --> 00:24:51,410 tells us how the level sets evolve, and then, at any time, 351 00:24:51,410 --> 00:24:55,760 we read off the position of the curve from the level set. 352 00:24:55,760 --> 00:25:01,440 And then one more step to that equation, 353 00:25:01,440 --> 00:25:06,550 which is often easier computationally than this one. 354 00:25:06,550 --> 00:25:09,300 So I'm going to take a simple step between there 355 00:25:09,300 --> 00:25:22,270 and there in the case when v is in the normal direction. 356 00:25:22,270 --> 00:25:28,160 So let me look at this v dot grad phi 357 00:25:28,160 --> 00:25:30,180 that came into the convection equation. 358 00:25:30,180 --> 00:25:36,050 I'll divide and multiply by the magnitude of the gradient, 359 00:25:36,050 --> 00:25:39,160 and then I recognize this as the normal vector. 360 00:25:39,160 --> 00:25:41,250 I recognize that as n. 361 00:25:41,250 --> 00:25:46,160 So this is v dot n, and I'm going to give a name to v dot n 362 00:25:46,160 --> 00:25:50,790 and call it F. That's the speed, speed in the normal direction, 363 00:25:50,790 --> 00:25:52,290 v dot n. 364 00:25:52,290 --> 00:25:54,740 And then there's this factor magnitude 365 00:25:54,740 --> 00:25:57,590 of grad phi coming along. 366 00:25:57,590 --> 00:25:59,510 OK, so that's the equation. 367 00:25:59,510 --> 00:26:03,410 That's sort of the famous level set equation 368 00:26:03,410 --> 00:26:05,990 where, if F is a constant, that's the case 369 00:26:05,990 --> 00:26:09,510 we're speaking about of a fire moving. 370 00:26:09,510 --> 00:26:13,940 But if F depends on the -- F could be a function of -- 371 00:26:13,940 --> 00:26:18,270 if I put parentheses in there, F could be a function of gradient 372 00:26:18,270 --> 00:26:22,650 phi, a more complicated function. 373 00:26:22,650 --> 00:26:26,540 For example, I wrote down, because I would otherwise 374 00:26:26,540 --> 00:26:36,690 forget, this could be our -- this is a possible F of grad 375 00:26:36,690 --> 00:26:39,680 phi. 376 00:26:39,680 --> 00:26:42,730 Actually, it depends on more than the magnitude 377 00:26:42,730 --> 00:26:48,140 so I'd better -- yeah. 378 00:26:48,140 --> 00:26:53,070 And certainly not a linear one because of that denominator. 379 00:26:53,070 --> 00:26:59,890 So a big set of examples in this subject, not linear, 380 00:26:59,890 --> 00:27:03,940 are curvature shapes. 381 00:27:03,940 --> 00:27:06,950 Actually, the curvature problem is quite interesting. 382 00:27:06,950 --> 00:27:11,240 If the movement depends on the curvature -- 383 00:27:11,240 --> 00:27:17,690 so suppose I start with an elliptical-looking shape, 384 00:27:17,690 --> 00:27:21,320 the curvature is smaller here, right? 385 00:27:21,320 --> 00:27:26,600 Oh, I'm sorry, smaller curvature there, larger curvature there, 386 00:27:26,600 --> 00:27:32,640 and if we choose the right evolution, 387 00:27:32,640 --> 00:27:37,200 often it will come -- maybe this is a case where the movement is 388 00:27:37,200 --> 00:27:46,770 inward because the curvature is somehow the appropriate vectors 389 00:27:46,770 --> 00:27:48,010 pointing in. 390 00:27:48,010 --> 00:27:52,060 Here it would move in slowly because the curvature is small. 391 00:27:52,060 --> 00:27:54,710 Here it would move in faster because the curvature is 392 00:27:54,710 --> 00:27:57,280 larger. 393 00:27:57,280 --> 00:27:59,270 The curvature is tending to even out. 394 00:27:59,270 --> 00:28:03,850 And, in fact, after a certain -- as time goes on, any shape 395 00:28:03,850 --> 00:28:07,960 approaches a circle, and the circle shrinks and shrinks 396 00:28:07,960 --> 00:28:09,390 and disappears. 397 00:28:09,390 --> 00:28:15,000 So that this class of -- this sort of nonlinear group 398 00:28:15,000 --> 00:28:18,560 of problems has a nice geometric picture, 399 00:28:18,560 --> 00:28:21,550 that whatever shape you begin with, 400 00:28:21,550 --> 00:28:25,150 you could even begin with a really weird shape like 401 00:28:25,150 --> 00:28:29,840 a barbell or something. 402 00:28:29,840 --> 00:28:31,580 So this has zero curvature. 403 00:28:31,580 --> 00:28:34,570 It wouldn't move at the start. 404 00:28:34,570 --> 00:28:38,410 This would move in, so after a little while 405 00:28:38,410 --> 00:28:43,950 it would be like this, then the corner would move in. 406 00:28:43,950 --> 00:28:48,800 So it seems a little unlikely it's true that after a while, 407 00:28:48,800 --> 00:28:50,780 the shape would become nearly circular, 408 00:28:50,780 --> 00:28:54,020 get smaller and smaller, and the shape would vanish. 409 00:28:54,020 --> 00:28:56,580 Anyway, that a curvature problem. 410 00:28:56,580 --> 00:29:04,020 OK, so I guess following the main theme of the course 411 00:29:04,020 --> 00:29:10,370 would be write down a finite difference method for the level 412 00:29:10,370 --> 00:29:12,580 set equation. 413 00:29:12,580 --> 00:29:14,220 That's the natural thing to do. 414 00:29:14,220 --> 00:29:15,800 Write down a finite difference method 415 00:29:15,800 --> 00:29:17,440 for the level set equation. 416 00:29:17,440 --> 00:29:18,940 Solve that numerically. 417 00:29:18,940 --> 00:29:21,500 So that would be on the web page. 418 00:29:21,500 --> 00:29:23,280 The difference equation is going to be 419 00:29:23,280 --> 00:29:28,820 a sort of upwind equation, but it's 420 00:29:28,820 --> 00:29:32,740 a little messier in two variables 421 00:29:32,740 --> 00:29:36,820 to write out in detail here. 422 00:29:36,820 --> 00:29:38,720 You'll see it there on the web. 423 00:29:38,720 --> 00:29:44,650 OK, so you solve it, and proceed in time, 424 00:29:44,650 --> 00:29:50,420 and reinitialize when necessary or when it seems advisable, 425 00:29:50,420 --> 00:29:54,340 and you get a great picture of this evolving, 426 00:29:54,340 --> 00:29:57,010 moving interface. 427 00:29:57,010 --> 00:29:58,880 So problems of that kind. 428 00:29:58,880 --> 00:30:01,300 Or optimizing the interface. 429 00:30:01,300 --> 00:30:04,070 There are other problems where you might -- 430 00:30:04,070 --> 00:30:06,810 ones on the web where-- you might be, look, 431 00:30:06,810 --> 00:30:09,380 the unknown might be a shape. 432 00:30:09,380 --> 00:30:11,510 And you're looking for the shape that 433 00:30:11,510 --> 00:30:16,400 maybe makes an eigenvalue as small as possible, subject 434 00:30:16,400 --> 00:30:17,650 to some constraint. 435 00:30:17,650 --> 00:30:21,460 So problems for shapes. 436 00:30:21,460 --> 00:30:24,080 That's really the type of problem 437 00:30:24,080 --> 00:30:27,900 that level set methods are created for. 438 00:30:27,900 --> 00:30:35,790 OK, and these example show the connection -- actually, 439 00:30:35,790 --> 00:30:40,000 somehow the -- yeah, it's really that if I tried to connect it 440 00:30:40,000 --> 00:30:47,160 to conservation laws, to the previous lectures, 441 00:30:47,160 --> 00:30:58,880 this is our earlier guy, would be -- in one variable, 442 00:30:58,880 --> 00:31:06,250 let's say, this would correspond to the derivative of maybe phi 443 00:31:06,250 --> 00:31:12,550 or capital U. So the conservation laws for little u 444 00:31:12,550 --> 00:31:16,150 correspond to equations like this for phi. 445 00:31:19,480 --> 00:31:29,070 Well, or really, grad phi in more dimensions. 446 00:31:32,020 --> 00:31:35,040 OK, well, all right. 447 00:31:35,040 --> 00:31:39,480 So what else remains to say about this topic? 448 00:31:39,480 --> 00:31:43,400 Really, I guess my main point is to introduce you 449 00:31:43,400 --> 00:31:48,540 to the idea of level sets and level set method for when 450 00:31:48,540 --> 00:31:54,390 problems of moving interfaces and shape optimization come up. 451 00:31:58,330 --> 00:32:00,620 I guess I also wanted to say something 452 00:32:00,620 --> 00:32:02,940 about this fast marching method. 453 00:32:02,940 --> 00:32:06,710 OK, so what's that? 454 00:32:06,710 --> 00:32:10,640 That's a case where there's a very fast solution 455 00:32:10,640 --> 00:32:14,250 to this problem. 456 00:32:14,250 --> 00:32:19,030 And that case occurs when the wave front is always 457 00:32:19,030 --> 00:32:26,530 moving the same way, when the F holds its sign. 458 00:32:26,530 --> 00:32:28,170 So it would apply here. 459 00:32:28,170 --> 00:32:30,650 So I could use the fast marching method. 460 00:32:30,650 --> 00:32:34,410 Rather than the level set equation 461 00:32:34,410 --> 00:32:36,370 and a finite difference method for it, 462 00:32:36,370 --> 00:32:40,210 I could use a fast marching method for this problem 463 00:32:40,210 --> 00:32:47,160 because I have no danger of the shape turning, 464 00:32:47,160 --> 00:32:49,950 changing direction, coming back in on itself, 465 00:32:49,950 --> 00:32:53,880 all the crazy things that can happen, 466 00:32:53,880 --> 00:32:57,710 the things that are illustrated on the web. 467 00:32:57,710 --> 00:33:02,900 OK, so the fast marching method, well, maybe I -- 468 00:33:02,900 --> 00:33:06,010 so what's the problem in that nice case? 469 00:33:09,070 --> 00:33:15,350 I want to find -- I could put on a mesh. 470 00:33:15,350 --> 00:33:19,120 Yeah, let me take as a problem the computation of the distance 471 00:33:19,120 --> 00:33:22,720 function from this curve. 472 00:33:22,720 --> 00:33:25,340 OK, so I could create a mesh. 473 00:33:25,340 --> 00:33:34,090 So I'm going to take that problem and turn it into a -- 474 00:33:34,090 --> 00:33:42,660 discretize it in a natural way, and I'm searching 475 00:33:42,660 --> 00:33:44,210 for the distance function. 476 00:33:44,210 --> 00:33:51,010 So d of x and y is the distance that I now want to compute. 477 00:33:51,010 --> 00:33:55,790 It's the distance from the curve C. 478 00:33:55,790 --> 00:33:58,150 And let me make it a signed distance function 479 00:33:58,150 --> 00:34:03,320 so I'll make d negative in here and positive outside 480 00:34:03,320 --> 00:34:05,340 of the curve. 481 00:34:05,340 --> 00:34:05,840 OK. 482 00:34:10,750 --> 00:34:14,900 All right, so the idea is we could -- 483 00:34:14,900 --> 00:34:25,112 this fits our level set picture, but there's a direct way to do 484 00:34:25,112 --> 00:34:27,370 it. 485 00:34:27,370 --> 00:34:31,210 And it's a nice -- it's an algorithm that computer 486 00:34:31,210 --> 00:34:35,080 scientists probably created and appreciate, that -- 487 00:34:35,080 --> 00:34:37,680 we want to know the distance. 488 00:34:37,680 --> 00:34:40,640 We want to compute d at the mesh points. 489 00:34:40,640 --> 00:34:42,690 Let me put another mesh point. 490 00:34:42,690 --> 00:34:47,120 OK, so how do I -- so the question is, 491 00:34:47,120 --> 00:34:53,240 how do I get from -- let me put another y in there, too, 492 00:34:53,240 --> 00:34:56,140 and that created some more points. 493 00:34:56,140 --> 00:35:02,830 So how to -- what path -- you know, 494 00:35:02,830 --> 00:35:07,990 how to get the distance from the curve to each mesh point? 495 00:35:07,990 --> 00:35:11,450 Well, here's the idea, just in a nutshell. 496 00:35:11,450 --> 00:35:18,000 You keep a list of tentative distances 497 00:35:18,000 --> 00:35:20,330 to all the mesh points, and you look 498 00:35:20,330 --> 00:35:27,280 for the one that's closest, like it's probably maybe that guy. 499 00:35:27,280 --> 00:35:30,660 So maybe that's right on the curve, let's say. 500 00:35:30,660 --> 00:35:37,880 So the idea is I can -- that will be -- 501 00:35:37,880 --> 00:35:42,460 to get to this point, I'm never going to take a path that 502 00:35:42,460 --> 00:35:45,010 passes through the others. 503 00:35:45,010 --> 00:35:47,630 To get here, I might take a path that 504 00:35:47,630 --> 00:35:49,990 goes through an earlier point. 505 00:35:49,990 --> 00:35:54,740 But the nearest point, the work is done. 506 00:35:54,740 --> 00:36:00,990 So you successively sweep through the mesh points, 507 00:36:00,990 --> 00:36:05,970 pick the one whose current distance, tentative distance, 508 00:36:05,970 --> 00:36:09,150 is the smallest, and you say that is the distance. 509 00:36:09,150 --> 00:36:11,310 That's the right number. 510 00:36:11,310 --> 00:36:15,350 Then now you know how far it is to get to that point. 511 00:36:15,350 --> 00:36:20,230 So now you could, if you had to, update all the previous -- 512 00:36:20,230 --> 00:36:24,270 all the other distances by knowing the correct one, 513 00:36:24,270 --> 00:36:30,230 how distant it is, how distant that point is. 514 00:36:30,230 --> 00:36:32,270 OK, so then that one's settled. 515 00:36:32,270 --> 00:36:36,510 Now you look at your list, your updated list. 516 00:36:36,510 --> 00:36:40,860 You pick the smallest one in that, the closest one in that. 517 00:36:40,860 --> 00:36:44,870 You found the best way because things are never 518 00:36:44,870 --> 00:36:46,630 going to reverse in sign; you're never 519 00:36:46,630 --> 00:36:49,420 going to travel out and come back. 520 00:36:49,420 --> 00:36:52,660 So the closest one, you've got it. 521 00:36:52,660 --> 00:36:55,140 So that might be this or something, 522 00:36:55,140 --> 00:36:58,640 and then we would update how far it takes to get to those. 523 00:36:58,640 --> 00:37:03,840 So that fast marching method is a, you might say, 524 00:37:03,840 --> 00:37:12,020 fairly natural idea of how to compute distances when 525 00:37:12,020 --> 00:37:20,560 everything is measured, when our problem is -- has a fixed sign. 526 00:37:20,560 --> 00:37:25,310 OK, so that's my short comment on the fast marching method. 527 00:37:25,310 --> 00:37:30,310 My longer comment is on the level set method and level set 528 00:37:30,310 --> 00:37:32,600 equation, where I haven't written out 529 00:37:32,600 --> 00:37:37,350 the finite difference approximation that people 530 00:37:37,350 --> 00:37:44,500 typically use, and note on the web. 531 00:37:44,500 --> 00:37:48,500 All right, so that's level sets, which I probably 532 00:37:48,500 --> 00:37:51,140 will not come back to. 533 00:37:51,140 --> 00:37:52,720 It's kind of a one-shot lecture. 534 00:37:56,570 --> 00:37:58,570 So let me take the remaining time 535 00:37:58,570 --> 00:38:05,190 to get feedback on the homework. 536 00:38:05,190 --> 00:38:08,510 So were there four problems, possible candidate problems, 537 00:38:08,510 --> 00:38:09,280 for homework? 538 00:38:09,280 --> 00:38:11,740 Let me put those down again. 539 00:38:11,740 --> 00:38:15,570 I'm just going to ask first, who did which problem. 540 00:38:15,570 --> 00:38:19,640 So the first problem was the conservation law -- no, 541 00:38:19,640 --> 00:38:23,790 was the one-way-- OK, remind me. 542 00:38:23,790 --> 00:38:27,560 I don't remember the order I gave. 543 00:38:27,560 --> 00:38:29,590 So was Number 4 the conservation law? 544 00:38:29,590 --> 00:38:30,610 The nonlinear one? 545 00:38:30,610 --> 00:38:31,190 AUDIENCE: [INAUDIBLE] 546 00:38:31,190 --> 00:38:32,030 PROFESSOR: Yeah, OK. 547 00:38:32,030 --> 00:38:33,430 So the conservation law. 548 00:38:33,430 --> 00:38:38,300 I'll just write down u_t plus u*u_x equals 0. 549 00:38:38,300 --> 00:38:42,540 Oh yeah, Number 1 was the basic question of u_t equals c*u_x. 550 00:38:49,179 --> 00:38:49,970 Yeah, that's right. 551 00:38:49,970 --> 00:38:53,430 So that's a very important and not too obvious 552 00:38:53,430 --> 00:38:56,450 question, a good numerical question, of what to do. 553 00:38:56,450 --> 00:38:59,320 And Number 2 was the wave equation; was that right? 554 00:38:59,320 --> 00:39:01,042 AUDIENCE: Schrodinger's. 555 00:39:01,042 --> 00:39:02,500 PROFESSOR: Schrodinger's, Number 2. 556 00:39:02,500 --> 00:39:03,110 Oh, boy. 557 00:39:03,110 --> 00:39:03,880 OK. 558 00:39:03,880 --> 00:39:10,000 So i*u_xx, which I haven't discussed in class and I had -- 559 00:39:10,000 --> 00:39:13,630 the good question came up, how do you deal with these complex 560 00:39:13,630 --> 00:39:16,600 numbers because, of course, it's going to get complex right 561 00:39:16,600 --> 00:39:17,510 away. 562 00:39:17,510 --> 00:39:19,700 Well, one way that just occurs to me 563 00:39:19,700 --> 00:39:22,320 is take the real and imaginary parts 564 00:39:22,320 --> 00:39:25,390 so that we get two equations, two real equations. 565 00:39:25,390 --> 00:39:30,344 PROFESSOR: OK, Number 3 was -- what did I have -- 566 00:39:30,344 --> 00:39:31,510 AUDIENCE: The wave equation. 567 00:39:31,510 --> 00:39:32,940 PROFESSOR: Was the wave equation? 568 00:39:32,940 --> 00:39:34,040 OK, second order? 569 00:39:34,040 --> 00:39:37,200 AUDIENCE: First order. 570 00:39:37,200 --> 00:39:40,940 PROFESSOR: Oh, first-order wave equation? 571 00:39:40,940 --> 00:39:43,200 Oh, OK. 572 00:39:43,200 --> 00:39:46,380 So what was I -- OK, first-order wave equation like so. 573 00:39:50,910 --> 00:39:51,410 I'm sorry? 574 00:39:51,410 --> 00:39:53,680 AUDIENCE: [INAUDIBLE] 575 00:39:53,680 --> 00:39:56,370 PROFESSOR: c*u_x, OK, yep. 576 00:39:56,370 --> 00:40:00,430 All right, so was the question there to go back to the-- 577 00:40:00,430 --> 00:40:01,680 AUDIENCE: Compare. 578 00:40:01,680 --> 00:40:02,830 PROFESSOR: Compare, right. 579 00:40:02,830 --> 00:40:05,710 Compare the methods we've already established. 580 00:40:08,450 --> 00:40:14,900 And then Number 4 was to tackle some case of nonlinear equation 581 00:40:14,900 --> 00:40:19,390 where shocks and fans could appear. 582 00:40:19,390 --> 00:40:23,690 OK, maybe I just get a sense of who did what. 583 00:40:23,690 --> 00:40:25,480 Number 1 how many? 584 00:40:25,480 --> 00:40:26,220 Wow! 585 00:40:26,220 --> 00:40:30,980 Number 1 was popular and I think it's justified. 586 00:40:30,980 --> 00:40:35,770 So this is certainly the biggest choice. 587 00:40:35,770 --> 00:40:41,780 So let me come back and get some ideas of what you concluded. 588 00:40:41,780 --> 00:40:43,780 Number 2, did anybody tackle Schrodinger? 589 00:40:43,780 --> 00:40:44,430 Yes. 590 00:40:44,430 --> 00:40:45,857 A couple of complex guys. 591 00:40:45,857 --> 00:40:47,190 Well, I know you asked about it. 592 00:40:47,190 --> 00:40:49,660 Yep, OK, good. 593 00:40:49,660 --> 00:40:51,760 Anybody go back to these? 594 00:40:51,760 --> 00:40:52,700 For this? 595 00:40:52,700 --> 00:40:53,290 You did? 596 00:40:53,290 --> 00:40:56,850 OK, you may have done a little bit of both. 597 00:40:56,850 --> 00:41:01,790 And did anybody venture to nonlinear world? 598 00:41:01,790 --> 00:41:02,350 One. 599 00:41:02,350 --> 00:41:03,610 Good, OK. 600 00:41:03,610 --> 00:41:05,600 Nonlinear person. 601 00:41:05,600 --> 00:41:07,920 I'm glad to know that, right. 602 00:41:10,920 --> 00:41:13,900 So this is like -- this homework you could think of us as like 603 00:41:13,900 --> 00:41:16,710 prep for a possible project. 604 00:41:16,710 --> 00:41:19,760 Your project could grow out of this homework, 605 00:41:19,760 --> 00:41:23,800 or it could grow from other things that are coming. 606 00:41:23,800 --> 00:41:24,760 OK. 607 00:41:24,760 --> 00:41:28,310 All right, let's stick with this one then since a lot of people 608 00:41:28,310 --> 00:41:31,770 will have input on it. 609 00:41:31,770 --> 00:41:34,900 Let me ask -- well, I'll just ask naive questions 610 00:41:34,900 --> 00:41:39,600 because you've actually done the calculation. 611 00:41:39,600 --> 00:41:43,230 Is there a method you'd recommend? 612 00:41:43,230 --> 00:41:48,350 If you had a problem like that or, I mean, like for real, 613 00:41:48,350 --> 00:41:54,340 and you had to invest in a code and use it, 614 00:41:54,340 --> 00:41:55,610 just give an opinion. 615 00:41:59,600 --> 00:42:02,030 Or did you only try one method? 616 00:42:02,030 --> 00:42:04,480 Tell me something, positive or negative, 617 00:42:04,480 --> 00:42:09,630 about your experience with this. 618 00:42:09,630 --> 00:42:11,430 Yeah, thanks. 619 00:42:11,430 --> 00:42:13,340 AUDIENCE: [INAUDIBLE] 620 00:42:13,340 --> 00:42:15,590 PROFESSOR: Maybe you push -- in theory, 621 00:42:15,590 --> 00:42:18,170 there's some mike you push and this will be the first time 622 00:42:18,170 --> 00:42:19,260 we've ever done it. 623 00:42:19,260 --> 00:42:20,260 OK, terrific. 624 00:42:20,260 --> 00:42:20,760 Go ahead. 625 00:42:20,760 --> 00:42:24,922 AUDIENCE: [INAUDIBLE] 626 00:42:24,922 --> 00:42:26,130 PROFESSOR: Centered explicit. 627 00:42:26,130 --> 00:42:27,720 Let me write that down. 628 00:42:27,720 --> 00:42:29,390 AUDIENCE: [INAUDIBLE] 629 00:42:29,390 --> 00:42:30,570 PROFESSOR: OK, right. 630 00:42:30,570 --> 00:42:43,900 AUDIENCE: [INAUDIBLE] 631 00:42:43,900 --> 00:42:44,660 PROFESSOR: Yes. 632 00:42:44,660 --> 00:42:51,500 AUDIENCE: [INAUDIBLE] 633 00:42:51,500 --> 00:42:52,500 PROFESSOR: That's right. 634 00:42:52,500 --> 00:42:54,520 So, OK, the centered explicit one 635 00:42:54,520 --> 00:42:58,440 is computing this from three old values. 636 00:42:58,440 --> 00:43:00,130 We certainly need three because we've 637 00:43:00,130 --> 00:43:02,500 got a second derivative there. 638 00:43:02,500 --> 00:43:04,840 But then we have a certain amount 639 00:43:04,840 --> 00:43:08,750 of that second derivative and a certain amount of this one, 640 00:43:08,750 --> 00:43:12,350 and that can make one of those three numbers go negative. 641 00:43:12,350 --> 00:43:16,350 They always add the one because the constant is certainly 642 00:43:16,350 --> 00:43:18,460 a solution we want to keep. 643 00:43:18,460 --> 00:43:22,540 OK, so when a coefficient went negative, 644 00:43:22,540 --> 00:43:27,580 that doesn't mean instability; is that right? 645 00:43:27,580 --> 00:43:29,524 AUDIENCE: [INAUDIBLE] 646 00:43:29,524 --> 00:43:30,190 PROFESSOR: Yeah. 647 00:43:30,190 --> 00:43:32,897 AUDIENCE: [INAUDIBLE] 648 00:43:32,897 --> 00:43:33,480 PROFESSOR: OK. 649 00:43:33,480 --> 00:43:35,924 AUDIENCE: [INAUDIBLE] 650 00:43:35,924 --> 00:43:36,590 PROFESSOR: Fine. 651 00:43:36,590 --> 00:43:39,527 AUDIENCE: [INAUDIBLE] 652 00:43:39,527 --> 00:43:40,110 PROFESSOR: OK. 653 00:43:40,110 --> 00:43:41,700 AUDIENCE: [INAUDIBLE] 654 00:43:41,700 --> 00:43:43,940 PROFESSOR: Implicit had no trouble with that. 655 00:43:43,940 --> 00:43:45,280 Yeah, OK. 656 00:43:45,280 --> 00:43:48,040 Or it didn't matter what the Peclet number, or the cell 657 00:43:48,040 --> 00:43:50,295 Peclet number was, or whatever we called it. 658 00:43:50,295 --> 00:43:51,170 AUDIENCE: [INAUDIBLE] 659 00:43:57,590 --> 00:43:58,340 PROFESSOR: Uh-huh. 660 00:43:58,340 --> 00:44:02,932 AUDIENCE: [INAUDIBLE] 661 00:44:02,932 --> 00:44:03,640 PROFESSOR: I see. 662 00:44:03,640 --> 00:44:05,500 AUDIENCE: [INAUDIBLE] 663 00:44:05,500 --> 00:44:09,210 PROFESSOR: OK, if you went -- when you say implicit, 664 00:44:09,210 --> 00:44:12,480 it's just the diffusion term that's going implicit? 665 00:44:12,480 --> 00:44:15,360 OK, OK. 666 00:44:15,360 --> 00:44:17,590 So that should be a very stable process, 667 00:44:17,590 --> 00:44:23,640 but, of course, now we've got this term affecting it. 668 00:44:23,640 --> 00:44:27,230 So you actually found out that even the implicit one 669 00:44:27,230 --> 00:44:29,230 could break? 670 00:44:29,230 --> 00:44:31,344 AUDIENCE: [INAUDIBLE] 671 00:44:31,344 --> 00:44:32,010 PROFESSOR: Yeah. 672 00:44:32,010 --> 00:44:34,546 AUDIENCE: [INAUDIBLE] 673 00:44:34,546 --> 00:44:35,170 PROFESSOR: Yes. 674 00:44:35,170 --> 00:44:36,470 AUDIENCE: [INAUDIBLE] 675 00:44:36,470 --> 00:44:37,220 PROFESSOR: Really? 676 00:44:37,220 --> 00:44:38,040 AUDIENCE: [INAUDIBLE] 677 00:44:38,040 --> 00:44:38,350 PROFESSOR: Huh! 678 00:44:38,350 --> 00:44:38,980 AUDIENCE: [INAUDIBLE] 679 00:44:38,980 --> 00:44:40,440 PROFESSOR: Even the implicit one? 680 00:44:40,440 --> 00:44:44,284 AUDIENCE: [INAUDIBLE] 681 00:44:44,284 --> 00:44:45,200 PROFESSOR: Yeah, yeah. 682 00:44:45,200 --> 00:44:49,910 So on the implicit part -- so this was method one and this 683 00:44:49,910 --> 00:44:53,550 was method two, implicit for the diffusion, 684 00:44:53,550 --> 00:45:00,200 but you found that that could still be unstable. 685 00:45:00,200 --> 00:45:01,830 OK, yeah. 686 00:45:01,830 --> 00:45:05,870 Anybody else do this comparison? 687 00:45:05,870 --> 00:45:10,600 Can you push your mike button just to -- 688 00:45:10,600 --> 00:45:13,280 this is a hotshot room that we haven't taken advantage 689 00:45:13,280 --> 00:45:14,409 of, but it's our chance. 690 00:45:14,409 --> 00:45:16,950 AUDIENCE: I actually looked at the -- when I first did the -- 691 00:45:16,950 --> 00:45:21,000 looking at just the growth factors of explicit 692 00:45:21,000 --> 00:45:21,780 and implicit -- 693 00:45:21,780 --> 00:45:23,900 PROFESSOR: OK, looking at that G -- 694 00:45:23,900 --> 00:45:26,034 AUDIENCE: --on the convection term as well? 695 00:45:26,034 --> 00:45:26,700 PROFESSOR: Yeah. 696 00:45:26,700 --> 00:45:30,060 AUDIENCE: And I used implicit for all the diffusion. 697 00:45:30,060 --> 00:45:34,640 And I found the most stable using upwind, 698 00:45:34,640 --> 00:45:37,290 when both were implicit, the growth factor 699 00:45:37,290 --> 00:45:39,560 was actually able to stay below 1. 700 00:45:39,560 --> 00:45:44,420 PROFESSOR: So you moved the first derivative also 701 00:45:44,420 --> 00:45:45,830 up to the new time? 702 00:45:45,830 --> 00:45:47,961 OK, and that was -- then you don't have a problem. 703 00:45:47,961 --> 00:45:48,960 AUDIENCE: Yeah, exactly. 704 00:45:48,960 --> 00:45:52,870 So stable for pretty much all big R, little r. 705 00:45:52,870 --> 00:45:55,970 But the explicit upwind for the convection term 706 00:45:55,970 --> 00:46:00,792 was stable for most, as long as the little r was below 0.5. 707 00:46:00,792 --> 00:46:01,500 PROFESSOR: I see. 708 00:46:01,500 --> 00:46:03,496 So the numbers could tell you the-- 709 00:46:03,496 --> 00:46:04,120 AUDIENCE: Sure. 710 00:46:04,120 --> 00:46:06,490 PROFESSOR: --limits on little r and big R? 711 00:46:06,490 --> 00:46:07,170 AUDIENCE: Yep. 712 00:46:07,170 --> 00:46:11,180 And then I did the explicit convection 713 00:46:11,180 --> 00:46:13,500 upwind with the implicit centered diffusion, 714 00:46:13,500 --> 00:46:15,980 and that seemed to be pretty stable unless I really started 715 00:46:15,980 --> 00:46:17,132 cranking the numbers out. 716 00:46:17,132 --> 00:46:17,840 PROFESSOR: I see. 717 00:46:17,840 --> 00:46:20,270 Oh, that's the one that he had also done? 718 00:46:20,270 --> 00:46:21,260 AUDIENCE: Yep. 719 00:46:21,260 --> 00:46:24,380 PROFESSOR: Yeah, OK. 720 00:46:24,380 --> 00:46:31,610 And now, so I guess I remember, in the figure that isn't yet 721 00:46:31,610 --> 00:46:34,620 drawn that might emerge from your calculations, 722 00:46:34,620 --> 00:46:39,990 a possibility where -- depending on the Peclet number, 723 00:46:39,990 --> 00:46:41,750 where it was stable. 724 00:46:41,750 --> 00:46:44,950 It didn't blow up, but oscillations appeared 725 00:46:44,950 --> 00:46:49,580 and the calculations -- you wouldn't do it even so -- 726 00:46:49,580 --> 00:46:50,795 AUDIENCE: [INAUDIBLE] 727 00:46:50,795 --> 00:46:52,420 PROFESSOR: That's what you meant, yeah. 728 00:46:52,420 --> 00:46:58,070 OK, so there are cases where the method isn't good 729 00:46:58,070 --> 00:47:04,030 even though it's stable, which is like important to realize 730 00:47:04,030 --> 00:47:07,280 that that can happen. 731 00:47:07,280 --> 00:47:11,900 Because up to now, stability has been, you know, the holy grail. 732 00:47:11,900 --> 00:47:15,210 But it's really a good solution that's the holy grail. 733 00:47:15,210 --> 00:47:17,090 Yes? 734 00:47:17,090 --> 00:47:17,600 Go ahead. 735 00:47:17,600 --> 00:47:21,770 AUDIENCE: I did kind of a bigger picture where I basically 736 00:47:21,770 --> 00:47:25,670 tried to build the pieces so just 737 00:47:25,670 --> 00:47:28,520 to isolate the advection term and the diffusion term 738 00:47:28,520 --> 00:47:29,905 and then build their matrices. 739 00:47:29,905 --> 00:47:30,530 PROFESSOR: Yes. 740 00:47:30,530 --> 00:47:32,530 AUDIENCE: And then look at the combination rule. 741 00:47:32,530 --> 00:47:34,902 So when you went to where you had both terms present-- 742 00:47:34,902 --> 00:47:35,610 PROFESSOR: Right. 743 00:47:35,610 --> 00:47:38,650 AUDIENCE: --how the matrices were actually -- 744 00:47:38,650 --> 00:47:40,650 came together to form the-- 745 00:47:40,650 --> 00:47:45,640 PROFESSOR: So you took a matrix picture of the -- 746 00:47:45,640 --> 00:47:49,930 with all the mesh points, mesh values, at once rather than-- 747 00:47:49,930 --> 00:47:52,800 AUDIENCE: And then I simulated with them just convection, 748 00:47:52,800 --> 00:47:55,840 just diffusion and then the coupled problem of both of them 749 00:47:55,840 --> 00:48:00,430 to see how that -- just one method may be unstable 750 00:48:00,430 --> 00:48:01,562 for a given formulation-- 751 00:48:01,562 --> 00:48:02,270 PROFESSOR: Right. 752 00:48:02,270 --> 00:48:04,436 AUDIENCE: --but yet when you saw the coupled problem 753 00:48:04,436 --> 00:48:06,540 with those two pieces, it actually is stable. 754 00:48:06,540 --> 00:48:09,280 So you could look at like the -- when you build the matrices, 755 00:48:09,280 --> 00:48:12,111 you can look at the eigenvalues and you can make sure that 756 00:48:12,111 --> 00:48:12,610 they're-- 757 00:48:12,610 --> 00:48:15,570 PROFESSOR: Let me comment a little on that. 758 00:48:15,570 --> 00:48:19,340 I'll do that, but I don't want to cut off anybody who has 759 00:48:19,340 --> 00:48:23,530 another comment to throw in here before I say something about -- 760 00:48:23,530 --> 00:48:25,920 follow up on matrices. 761 00:48:25,920 --> 00:48:29,190 Anybody got some urgent thing? 762 00:48:29,190 --> 00:48:31,130 Let me just, in the remaining minute, 763 00:48:31,130 --> 00:48:34,180 just say a word because, of course, 764 00:48:34,180 --> 00:48:38,070 you'll know that I think of the matrix description 765 00:48:38,070 --> 00:48:41,180 quickly and fondly. 766 00:48:41,180 --> 00:48:46,800 So just to write down the sort of matrices that are in here, 767 00:48:46,800 --> 00:48:51,410 from the second difference, the second difference matrix is 768 00:48:51,410 --> 00:48:56,120 going to be -- it's going to have our -- well, 769 00:48:56,120 --> 00:49:01,490 I guess it'll be 1, minus 2, and 1 in a typical row. 770 00:49:01,490 --> 00:49:04,040 We haven't written any matrices down. 771 00:49:04,040 --> 00:49:11,400 It's amazing that -- we'd be at about the, I don't know, 772 00:49:11,400 --> 00:49:15,460 eighth or ninth lecture now, and not have seen this matrix, 773 00:49:15,460 --> 00:49:16,810 but here it is. 774 00:49:16,810 --> 00:49:21,820 OK, so that's the matrix where I'm at a fixed time level. 775 00:49:21,820 --> 00:49:26,110 I'm taking the second differences of all the points 776 00:49:26,110 --> 00:49:27,570 along that level. 777 00:49:27,570 --> 00:49:31,150 And now, the main point for the moment 778 00:49:31,150 --> 00:49:35,850 is to contrast that with a first difference matrix 779 00:49:35,850 --> 00:49:40,040 which, well, that depends whether I take it centered. 780 00:49:40,040 --> 00:49:42,710 Maybe if it's a centered difference, 781 00:49:42,710 --> 00:49:45,260 it would be a 1, 0, minus 1. 782 00:49:45,260 --> 00:49:50,420 And, of course, there's a 2 delta x in the denominator. 783 00:49:50,420 --> 00:49:53,880 so a minus 1, 0, 1. 784 00:49:53,880 --> 00:49:59,370 So that would be a centered difference that 785 00:49:59,370 --> 00:50:02,900 doesn't do anything, doesn't use the middle value, 786 00:50:02,900 --> 00:50:06,360 but takes the difference between the value at the right 787 00:50:06,360 --> 00:50:07,570 and the value at the left. 788 00:50:07,570 --> 00:50:09,860 Or, an upwind. 789 00:50:09,860 --> 00:50:18,830 So delta_upwind would be more of a 1's forward of the point 790 00:50:18,830 --> 00:50:26,570 and minus 1's at u_(j, n) and 0's here. 791 00:50:26,570 --> 00:50:28,890 OK. 792 00:50:28,890 --> 00:50:31,420 So what you're saying, am I right 793 00:50:31,420 --> 00:50:36,000 in thinking you took some combination of the-- 794 00:50:39,350 --> 00:50:43,290 AUDIENCE: An implicit method for one and explicit for the other, 795 00:50:43,290 --> 00:50:47,300 you can do them independently by -- 796 00:50:47,300 --> 00:50:50,042 you can invert the implicit matrix-- 797 00:50:50,042 --> 00:50:50,750 PROFESSOR: Right. 798 00:50:50,750 --> 00:50:53,510 AUDIENCE: --and pre-multiply the other side and build-- 799 00:50:53,510 --> 00:50:55,150 PROFESSOR: So this is a way to see, 800 00:50:55,150 --> 00:50:58,510 and actually the next months of our course 801 00:50:58,510 --> 00:51:01,990 are going to be about, how do you 802 00:51:01,990 --> 00:51:04,920 deal with these large difference matrices. 803 00:51:04,920 --> 00:51:09,060 These are in one dimension so they're tridiagonal, 804 00:51:09,060 --> 00:51:11,090 and they're, of course, easy to deal with it. 805 00:51:11,090 --> 00:51:14,030 But problems, serious problems, come up in more dimension. 806 00:51:14,030 --> 00:51:20,590 But here, can I just express one thought here? 807 00:51:20,590 --> 00:51:24,340 The eigenvalues of these matrices 808 00:51:24,340 --> 00:51:26,510 and their combinations, as you say, 809 00:51:26,510 --> 00:51:30,310 are what's going to control, as eigenvalues always do, 810 00:51:30,310 --> 00:51:32,540 growth of solution. 811 00:51:32,540 --> 00:51:37,970 Now we already had numbers G, growth factors, 812 00:51:37,970 --> 00:51:43,750 that we computed by the von Neumann idea of following 813 00:51:43,750 --> 00:51:46,100 exponentials, e to the i*k*x. 814 00:51:46,100 --> 00:51:51,270 So these numbers too were -- their purpose was to track 815 00:51:51,270 --> 00:51:54,360 the growth of solution. 816 00:51:54,360 --> 00:51:57,730 And what I want to say is that very often, I mean, 817 00:51:57,730 --> 00:52:03,460 that these numbers and the eigenvalues of the matrices, 818 00:52:03,460 --> 00:52:07,760 they're very closely related ideas. 819 00:52:07,760 --> 00:52:09,960 In fact, they'll be the same idea 820 00:52:09,960 --> 00:52:14,050 if the eigenvector is a pure exponential. 821 00:52:14,050 --> 00:52:17,950 So that's really the link between eigenvalues 822 00:52:17,950 --> 00:52:25,500 of matrices and growth factors G that depended on frequency k. 823 00:52:25,500 --> 00:52:32,890 If the eigenvector was a pure e to the i*k*x on the mesh, 824 00:52:32,890 --> 00:52:39,140 then that growth factor is identical to the eigenvalue 825 00:52:39,140 --> 00:52:43,290 that we get in one step of the difference equation. 826 00:52:43,290 --> 00:52:49,280 So the idea was that the eigenvectors, the details 827 00:52:49,280 --> 00:52:51,340 of the eigenvectors, will be affected 828 00:52:51,340 --> 00:52:53,070 by boundary conditions. 829 00:52:53,070 --> 00:52:54,870 So von Neumann thought, all right, 830 00:52:54,870 --> 00:53:00,590 I'll cut through that Gordian knot and just pretend I have 831 00:53:00,590 --> 00:53:02,790 a pure e to the i*k*x. 832 00:53:02,790 --> 00:53:06,440 I look inside the region, not at the boundary where 833 00:53:06,440 --> 00:53:08,930 we may be doing something different, 834 00:53:08,930 --> 00:53:12,440 and we get a G, a growth factor G. 835 00:53:12,440 --> 00:53:16,980 And again, if the eigenfunction was really an eigenvector, then 836 00:53:16,980 --> 00:53:19,540 that G would be the same as lambda. 837 00:53:19,540 --> 00:53:21,220 OK, we're past time. 838 00:53:21,220 --> 00:53:24,760 I always appreciate that you are patient about that. 839 00:53:24,760 --> 00:53:29,640 So I'll collect those that are complete. 840 00:53:29,640 --> 00:53:32,940 And if you wanted to speak to Mr. Cho 841 00:53:32,940 --> 00:53:37,870 or have a little more time to do more, Monday is good, too. 842 00:53:37,870 --> 00:53:39,460 OK, thanks.