1 00:00:00,000 --> 00:00:01,960 The following content is provided 2 00:00:01,960 --> 00:00:06,310 by MIT OpenCourseWare under a Creative Commons license. 3 00:00:06,310 --> 00:00:08,690 Additional information about our license 4 00:00:08,690 --> 00:00:10,490 and MIT OpenCourseWare in general 5 00:00:10,490 --> 00:00:11,800 is available at ocw.mit.edu. 6 00:00:15,122 --> 00:00:16,580 PROFESSOR: All right, so I'll start 7 00:00:16,580 --> 00:00:23,270 with this review board, which shows the conservation law. 8 00:00:23,270 --> 00:00:27,140 Scalar, one unknown, one dimension, 9 00:00:27,140 --> 00:00:34,310 x, one space variable, and the flux function is f of u. 10 00:00:34,310 --> 00:00:37,100 That's the differential form. 11 00:00:37,100 --> 00:00:42,900 We have seen that at some points that differential form 12 00:00:42,900 --> 00:00:48,400 produces a discontinuous solution, a discontinuity. 13 00:00:48,400 --> 00:00:51,650 In that case we are better off with 14 00:00:51,650 --> 00:00:54,720 the integrated integral form. 15 00:00:54,720 --> 00:01:00,610 So I just integrate that one from a point x_left to a point 16 00:01:00,610 --> 00:01:04,900 x_right, from an a to b, and that gives me-- I 17 00:01:04,900 --> 00:01:07,620 have then the integral. 18 00:01:07,620 --> 00:01:11,070 So if u was density, for example, 19 00:01:11,070 --> 00:01:13,630 the integral of density will be the mass. 20 00:01:13,630 --> 00:01:17,530 That's the quantity that I'm seeing here, 21 00:01:17,530 --> 00:01:20,210 so this is the rate of change of mass 22 00:01:20,210 --> 00:01:29,310 in the interval comes from flux going out minus flux coming in. 23 00:01:29,310 --> 00:01:34,530 So that's the conservation law in its pure statement. 24 00:01:34,530 --> 00:01:40,880 And we learned that solution for this-- 25 00:01:40,880 --> 00:01:43,570 while we're looking at the differential problem, 26 00:01:43,570 --> 00:01:46,170 we can follow these characteristics 27 00:01:46,170 --> 00:01:49,490 from the initial time zero. 28 00:01:49,490 --> 00:01:55,960 So at time zero we look at a point x and the value 29 00:01:55,960 --> 00:02:01,090 at that point just travels-- is sort of safely 30 00:02:01,090 --> 00:02:07,510 carried up the characteristic line, which is that equals x_0. 31 00:02:07,510 --> 00:02:14,970 And the problem arose when characteristics either collided 32 00:02:14,970 --> 00:02:17,690 into each other or separated. 33 00:02:21,340 --> 00:02:23,480 Those are the two difficulties and they 34 00:02:23,480 --> 00:02:28,480 have different resolutions. 35 00:02:28,480 --> 00:02:32,370 So this is what we met in the traffic example 36 00:02:32,370 --> 00:02:35,580 when a light goes red. 37 00:02:35,580 --> 00:02:41,160 When a light goes red cars are coming up to the light, 38 00:02:41,160 --> 00:02:46,950 in front of the light the density is low. 39 00:02:46,950 --> 00:02:49,410 The cars have to brake quickly. 40 00:02:49,410 --> 00:02:52,600 So something happens-- that's a question mark there. 41 00:02:52,600 --> 00:02:55,060 What do you do when characteristics collide? 42 00:02:55,060 --> 00:02:59,430 When the value here cannot be both the initial value 43 00:02:59,430 --> 00:03:02,380 that's coming up that characteristic and the initial 44 00:03:02,380 --> 00:03:05,720 value that's coming up there because they're different? 45 00:03:05,720 --> 00:03:08,180 So we need a rule for what to do there. 46 00:03:08,180 --> 00:03:10,780 And we also need a rule for what to do 47 00:03:10,780 --> 00:03:13,600 in case of an empty space. 48 00:03:13,600 --> 00:03:16,920 When no characteristics-- here, two characteristics 49 00:03:16,920 --> 00:03:20,160 are hitting, here none are touching there. 50 00:03:20,160 --> 00:03:23,280 Characteristics are separating in this case and the question 51 00:03:23,280 --> 00:03:25,540 is, what solution do you fill in the middle 52 00:03:25,540 --> 00:03:27,580 when you haven't got a characteristic? 53 00:03:27,580 --> 00:03:32,030 It's not coming from the start. 54 00:03:32,030 --> 00:03:35,590 OK, what solution to choose? 55 00:03:35,590 --> 00:03:38,210 The answers will be different. 56 00:03:38,210 --> 00:03:44,490 And let me take the first one, the collision question, first. 57 00:03:44,490 --> 00:03:52,850 So what will happen is a shock will form. 58 00:03:52,850 --> 00:03:59,330 There'll be a shock line so that these characteristics carry 59 00:03:59,330 --> 00:04:04,500 information into the shock, by the way, that's a big deal-- 60 00:04:04,500 --> 00:04:07,710 up until they hit the shock, and these characteristics 61 00:04:07,710 --> 00:04:08,970 carry their information. 62 00:04:08,970 --> 00:04:14,110 And across the shock we have different values, 63 00:04:14,110 --> 00:04:16,960 u on the left and u on the right. 64 00:04:16,960 --> 00:04:23,490 So that's one possible form of the solution 65 00:04:23,490 --> 00:04:26,660 because that will solve the equation in this region 66 00:04:26,660 --> 00:04:28,760 where the characteristics from the right 67 00:04:28,760 --> 00:04:30,920 are following the formula. 68 00:04:30,920 --> 00:04:35,060 It'll solve it from the left, so everything comes down 69 00:04:35,060 --> 00:04:39,650 to the jump condition at the shock. 70 00:04:39,650 --> 00:04:41,260 What's the shock speed? 71 00:04:41,260 --> 00:04:44,230 Where does this shock go? 72 00:04:44,230 --> 00:04:47,500 It's to find that shock in the xt-plane 73 00:04:47,500 --> 00:04:52,970 that I have to use this integral form. 74 00:04:52,970 --> 00:04:54,730 I must use the integral form because I've 75 00:04:54,730 --> 00:04:57,880 got a discontinuity and the differential form just 76 00:04:57,880 --> 00:04:59,280 doesn't make sense. 77 00:04:59,280 --> 00:05:04,790 OK, I wrote down the key idea here or the conclusion, but not 78 00:05:04,790 --> 00:05:06,360 the reasoning. 79 00:05:06,360 --> 00:05:08,460 I'm looking for this jump condition 80 00:05:08,460 --> 00:05:13,370 at the shock, which is going to locate the shock. 81 00:05:13,370 --> 00:05:16,230 And the key quantity that I have to find 82 00:05:16,230 --> 00:05:19,040 is the shock speed, which will tell me 83 00:05:19,040 --> 00:05:24,520 where that shock is going, how it's moving up in time. 84 00:05:24,520 --> 00:05:30,010 So time is going that way, x is going this way in this picture. 85 00:05:30,010 --> 00:05:33,400 So this has a speed, s of t. 86 00:05:36,790 --> 00:05:38,330 And here's the conclusion. 87 00:05:38,330 --> 00:05:41,270 There's the jump condition. 88 00:05:41,270 --> 00:05:43,730 So what does it mean? 89 00:05:43,730 --> 00:05:48,190 This notation, this bracket notation is the jump. 90 00:05:51,640 --> 00:05:53,350 The shock equation is beautiful. 91 00:05:53,350 --> 00:05:57,230 It says that s, the shock speed, times the jump 92 00:05:57,230 --> 00:06:00,510 in u is the jump in f of u. 93 00:06:00,510 --> 00:06:03,640 And let me take this example that 94 00:06:03,640 --> 00:06:07,950 is sort of the Burgers' equation example. 95 00:06:07,950 --> 00:06:11,360 So this is the example that comes from Burgers' equation. 96 00:06:18,260 --> 00:06:20,730 That's the example where the flux 97 00:06:20,730 --> 00:06:27,420 is this nice parabolic function u squared over 2. 98 00:06:27,420 --> 00:06:32,520 Then the change in flux, the jump in the flux function 99 00:06:32,520 --> 00:06:36,630 is u right squared minus u left squared over 2 100 00:06:36,630 --> 00:06:41,330 and I divide by the jump in u, which is u_right minus u_left 101 00:06:41,330 --> 00:06:44,600 and that cancels beautifully and leaves me 102 00:06:44,600 --> 00:06:48,680 with the 1/2 of u_right plus u_left. 103 00:06:48,680 --> 00:06:54,270 It tells me that the shock speed is exactly halfway between. 104 00:06:54,270 --> 00:06:58,820 It's exactly the average, in this Burger's example, 105 00:06:58,820 --> 00:07:02,420 of the characteristic speed from the right 106 00:07:02,420 --> 00:07:04,720 and the characteristics speed from the left. 107 00:07:04,720 --> 00:07:07,770 And notice, it's in between. 108 00:07:07,770 --> 00:07:09,320 That's crucial. 109 00:07:09,320 --> 00:07:12,130 That tells us that this figure is good. 110 00:07:12,130 --> 00:07:17,010 It tells us that characteristics are entering the shock 111 00:07:17,010 --> 00:07:18,440 from both sides. 112 00:07:18,440 --> 00:07:20,550 Characteristic lines are entering the shock. 113 00:07:20,550 --> 00:07:25,460 The shock is somehow, in this sort of nonlinear case, 114 00:07:25,460 --> 00:07:31,070 it's somehow kind of cutting out bits of the solution, 115 00:07:31,070 --> 00:07:33,590 reducing the energy, reducing the entropy, 116 00:07:33,590 --> 00:07:37,830 reducing the total variation as we'll see it in a minute. 117 00:07:37,830 --> 00:07:41,740 So that's the answer here and I guess I have to justify, 118 00:07:41,740 --> 00:07:44,880 where did this come from? 119 00:07:44,880 --> 00:07:49,600 Well, it came from this integral form. 120 00:07:52,210 --> 00:07:55,750 Suppose I integrate from a little bit 121 00:07:55,750 --> 00:07:58,117 to the left of the shock to a little bit 122 00:07:58,117 --> 00:07:59,200 to the right of the shock. 123 00:07:59,200 --> 00:08:02,630 Of course, that's where I want to integrate in order 124 00:08:02,630 --> 00:08:04,540 to capture the shock. 125 00:08:04,540 --> 00:08:07,700 OK, so if I just take that equation 126 00:08:07,700 --> 00:08:12,510 and write it over here, and I'll imagine-- because I'm not 127 00:08:12,510 --> 00:08:15,410 going to be-- x_left and x_right are going to be pretty near 128 00:08:15,410 --> 00:08:19,590 the shock-- u up to the shock will pretty much 129 00:08:19,590 --> 00:08:21,630 have a constant value u_left. 130 00:08:21,630 --> 00:08:26,960 So approximately, that is the shock position. 131 00:08:26,960 --> 00:08:30,600 Let me say the shock is at-- capital 132 00:08:30,600 --> 00:08:32,100 X is the position of the shock. 133 00:08:32,100 --> 00:08:42,460 So I have x minus x_left times u on the left of the shock. 134 00:08:45,790 --> 00:08:47,570 I'm doing that integral. 135 00:08:47,570 --> 00:08:51,450 Oh, I need d by dt of all this. 136 00:08:51,450 --> 00:08:52,200 So d by dt. 137 00:08:55,000 --> 00:08:57,540 It's d by dt of the integral of u dx. 138 00:08:57,540 --> 00:09:01,520 Let me just put here what I'm computing. 139 00:09:01,520 --> 00:09:05,840 I'm computing that from x_left to x_right. 140 00:09:05,840 --> 00:09:08,560 Well, I'm going up to x, up to capital 141 00:09:08,560 --> 00:09:10,880 X where the shock is, and then onward. 142 00:09:10,880 --> 00:09:15,260 So here I went up to x using the value u left. 143 00:09:15,260 --> 00:09:18,750 And now will be the onward part, which 144 00:09:18,750 --> 00:09:25,050 will be x_right minus x, the half of the interval that's 145 00:09:25,050 --> 00:09:27,250 on the right of the shock, times u right. 146 00:09:30,160 --> 00:09:38,430 And then copying this term is exactly our jump in f of u. 147 00:09:38,430 --> 00:09:46,530 Plus the jump in f of u is zero. 148 00:09:49,610 --> 00:09:53,020 Those words, total variation are a reminder to myself 149 00:09:53,020 --> 00:09:54,270 for what comes next. 150 00:09:54,270 --> 00:09:58,330 OK, so is that OK? 151 00:09:58,330 --> 00:10:03,440 That will be approximately zero because u wasn't exactly 152 00:10:03,440 --> 00:10:06,700 constant on the left and right. 153 00:10:06,700 --> 00:10:09,340 It didn't stay at the value u_L but we're 154 00:10:09,340 --> 00:10:11,680 squeezing the interval down where 155 00:10:11,680 --> 00:10:13,560 it's essentially constant. 156 00:10:13,560 --> 00:10:20,850 OK, so now, just look and see what we've got. 157 00:10:20,850 --> 00:10:30,990 dx/dt is the shock speed, so that's s, multiplying u_L, 158 00:10:30,990 --> 00:10:34,680 and here I have-- the derivative of that is zero, 159 00:10:34,680 --> 00:10:37,230 the derivative of x_R is zero. 160 00:10:37,230 --> 00:10:39,750 Here I have a minus, so it's minus u_R. 161 00:10:42,540 --> 00:10:50,320 dx/dt times that number plus the jump in f of u, left to right, 162 00:10:50,320 --> 00:10:52,250 is 0. 163 00:10:52,250 --> 00:10:54,560 And that's our equation. 164 00:10:57,570 --> 00:11:00,790 If I put this on the opposite side 165 00:11:00,790 --> 00:11:02,710 it becomes u_right minus u_left. 166 00:11:02,710 --> 00:11:04,330 Shall I do that? 167 00:11:04,330 --> 00:11:09,360 So I'll keep this on this side, got a good sign, 168 00:11:09,360 --> 00:11:11,820 but I'll move this onto the other side 169 00:11:11,820 --> 00:11:16,030 so it will be the shock speed, d capital X dt times 170 00:11:16,030 --> 00:11:22,110 u_R minus u_L, which is the jump in u. 171 00:11:22,110 --> 00:11:23,790 So that's great. 172 00:11:23,790 --> 00:11:28,420 We now have a way to deal with shocks that collide. 173 00:11:31,740 --> 00:11:40,690 So that was a derivation of this jump condition, which first 174 00:11:40,690 --> 00:11:45,490 had a name, two authors Rankine, Hugoniot in gas dynamics, 175 00:11:45,490 --> 00:11:48,860 but then it was seen as the right condition 176 00:11:48,860 --> 00:11:52,420 in all conservation laws. 177 00:11:52,420 --> 00:11:53,170 OK. 178 00:11:53,170 --> 00:11:56,630 And now I'm ready to comment on the second issue. 179 00:11:56,630 --> 00:12:00,350 What do you do when characteristics separate? 180 00:12:00,350 --> 00:12:03,220 Well, let me say what you don't do. 181 00:12:03,220 --> 00:12:08,560 You could, without wrecking the equation, 182 00:12:08,560 --> 00:12:12,440 you could try to put a shock in there 183 00:12:12,440 --> 00:12:15,680 and give this the value u on the left, 184 00:12:15,680 --> 00:12:18,250 even though a characteristic isn't 185 00:12:18,250 --> 00:12:22,450 telling us what's going on in this space or this one. 186 00:12:22,450 --> 00:12:24,810 You could try to put a shock in there. 187 00:12:24,810 --> 00:12:26,900 You could maybe satisfy the jump condition 188 00:12:26,900 --> 00:12:31,740 if you put the right shock in, but the characteristics would 189 00:12:31,740 --> 00:12:36,150 be coming out and that's wrong. 190 00:12:36,150 --> 00:12:40,400 And it's called the entropy condition, which 191 00:12:40,400 --> 00:12:42,740 says that this can't happen. 192 00:12:42,740 --> 00:12:46,390 That the derivative of entropy has a certain sign 193 00:12:46,390 --> 00:12:48,860 and it rules this out. 194 00:12:48,860 --> 00:12:52,230 It rules out this possibility of characteristics 195 00:12:52,230 --> 00:12:55,330 emerging from the shock and we need a totally different 196 00:12:55,330 --> 00:12:56,170 solution. 197 00:12:56,170 --> 00:12:59,050 And it turns out to be a fan in there. 198 00:12:59,050 --> 00:13:04,710 So we put in a fan of characteristics 199 00:13:04,710 --> 00:13:09,500 coming from here and that gives a solution that actually 200 00:13:09,500 --> 00:13:14,490 will depend on x over t. 201 00:13:14,490 --> 00:13:23,870 This solution will be a simple-- a simple ratio x over t 202 00:13:23,870 --> 00:13:24,820 enters that. 203 00:13:28,500 --> 00:13:30,910 In other words, that gives us a smooth solution, 204 00:13:30,910 --> 00:13:34,410 and I guess, thinking about the traffic example, that's 205 00:13:34,410 --> 00:13:35,640 exactly what happens. 206 00:13:35,640 --> 00:13:37,650 So this is the shock part. 207 00:13:37,650 --> 00:13:42,850 Now for the fan part, what happens 208 00:13:42,850 --> 00:13:48,740 when we have a bunch of cars sitting there, ready to go? 209 00:13:48,740 --> 00:13:56,430 The light goes green, the first cars start moving, 210 00:13:56,430 --> 00:14:02,150 cars behind them see the car in front moving, they start moving 211 00:14:02,150 --> 00:14:11,140 and so the initial profile would be, in this case, 212 00:14:11,140 --> 00:14:20,180 for cars, would be high density and then zero density. 213 00:14:20,180 --> 00:14:28,800 So this would be maximum density sitting at the light and zero 214 00:14:28,800 --> 00:14:33,180 density in front of the light and then the light goes green 215 00:14:33,180 --> 00:14:35,180 now and what happens? 216 00:14:35,180 --> 00:14:38,920 Well, as we all know those first cars pull ahead 217 00:14:38,920 --> 00:14:47,670 and the profile is these cars haven't yet started to move. 218 00:14:47,670 --> 00:14:50,340 These have started to move and there's a profile 219 00:14:50,340 --> 00:14:56,160 there and then zero density. 220 00:14:56,160 --> 00:15:00,290 So this represents the first car, the lead car. 221 00:15:00,290 --> 00:15:08,900 In front of that car is zero density and that's the fan. 222 00:15:08,900 --> 00:15:14,600 That's the x over t expression. 223 00:15:14,600 --> 00:15:19,880 So a simple linear expression, which will flatten out; 224 00:15:19,880 --> 00:15:24,190 as time goes on, these cars will be going faster and faster 225 00:15:24,190 --> 00:15:28,440 and then eventually up to full speed. 226 00:15:28,440 --> 00:15:30,860 OK, so that's the fan. 227 00:15:34,950 --> 00:15:41,060 Maybe now is my chance to speak about this concept 228 00:15:41,060 --> 00:15:44,350 of the variation in a function. 229 00:15:44,350 --> 00:15:49,790 The idea will be now, how do we get the finite difference 230 00:15:49,790 --> 00:15:54,370 method to do this automatically? 231 00:15:54,370 --> 00:16:00,180 Well, OK, one way to do it is put in a little viscosity. 232 00:16:03,540 --> 00:16:10,180 This rule, the rule that picks out this shock in one case, 233 00:16:10,180 --> 00:16:17,810 of colliding characteristics, and this fan 234 00:16:17,810 --> 00:16:22,250 in the case of separating characteristics, one way 235 00:16:22,250 --> 00:16:25,520 to get to those two solution would 236 00:16:25,520 --> 00:16:32,560 be to put in a little bit of viscosity 237 00:16:32,560 --> 00:16:35,930 and then let epsilon go to zero. 238 00:16:35,930 --> 00:16:42,070 And the result would be these two solutions. 239 00:16:42,070 --> 00:16:45,950 And in a way, that's what finite differences is going to do. 240 00:16:45,950 --> 00:16:47,760 Finite difference or finite volume 241 00:16:47,760 --> 00:16:52,530 is going to be a numerical viscosity 242 00:16:52,530 --> 00:16:58,400 and the epsilon in the Burgers' equation 243 00:16:58,400 --> 00:17:05,844 would become delta t or a delta x in the numerical method, 244 00:17:05,844 --> 00:17:07,260 in the finite difference equation. 245 00:17:10,810 --> 00:17:12,600 And then what's the concept? 246 00:17:12,600 --> 00:17:15,740 So we have to find finite differences. 247 00:17:15,740 --> 00:17:25,100 And we want those-- the point is finite differences 248 00:17:25,100 --> 00:17:27,300 sort of takes more care. 249 00:17:27,300 --> 00:17:31,310 We can't just create something that's 250 00:17:31,310 --> 00:17:35,620 consistent in the smooth case. 251 00:17:35,620 --> 00:17:40,210 Up till now we've just had no hesitation. 252 00:17:40,210 --> 00:17:43,680 We replaced derivatives by differences 253 00:17:43,680 --> 00:17:46,270 without too many worries, really. 254 00:17:49,040 --> 00:17:51,500 And that was fine as long as the solution is smooth. 255 00:17:51,500 --> 00:17:55,630 But now if we're going to replace derivatives 256 00:17:55,630 --> 00:17:59,370 by differences, we have to do it in such a way 257 00:17:59,370 --> 00:18:05,400 that this entropy condition, these solutions 258 00:18:05,400 --> 00:18:07,800 emerge automatically. 259 00:18:07,800 --> 00:18:15,420 And really, in a way that-- our finite difference 260 00:18:15,420 --> 00:18:18,780 or our finite volume methods should somehow 261 00:18:18,780 --> 00:18:24,100 conserve mass the same way the true equation does. 262 00:18:24,100 --> 00:18:27,020 So those are two different points. 263 00:18:27,020 --> 00:18:31,720 First point is this question of total variation. 264 00:18:31,720 --> 00:18:33,610 What's the total variation in a function? 265 00:18:36,240 --> 00:18:41,220 Well, let me draw a function. 266 00:18:41,220 --> 00:18:44,190 The total variation is-- well, let me see. 267 00:18:44,190 --> 00:18:50,480 Total variation in a function u is going to be, 268 00:18:50,480 --> 00:18:57,130 by definition, the integral of the absolute value of du/dx. 269 00:19:03,180 --> 00:19:09,060 OK now we have to understand this concept. 270 00:19:09,060 --> 00:19:21,350 It's beautiful that the analysis has led to a quantity as simple 271 00:19:21,350 --> 00:19:25,750 as that, as a guide. 272 00:19:25,750 --> 00:19:31,250 So when we create finite difference methods we will ask, 273 00:19:31,250 --> 00:19:36,260 are they total variation diminishing? 274 00:19:36,260 --> 00:19:42,060 So TVD will be total variation diminishing methods. 275 00:19:42,060 --> 00:19:46,600 Methods where this quantity does not increase in time. 276 00:19:49,180 --> 00:19:52,050 Let me just understand this quantity first. 277 00:19:52,050 --> 00:19:53,600 The absolute value is crucial. 278 00:19:53,600 --> 00:19:57,550 If the absolute value was not there, 279 00:19:57,550 --> 00:20:01,770 then I would just have the integral of the derivative, 280 00:20:01,770 --> 00:20:06,250 so it would be u at that point minus u at that. 281 00:20:06,250 --> 00:20:11,110 But here I'm taking absolute values, so for a while 282 00:20:11,110 --> 00:20:14,220 the derivative is positive, up to there. 283 00:20:14,220 --> 00:20:16,770 So the integral gives me this value minus this. 284 00:20:19,840 --> 00:20:24,160 Then, in this region, the derivative is negative, 285 00:20:24,160 --> 00:20:26,560 but I'm taking absolute value. 286 00:20:26,560 --> 00:20:30,450 Do you see that that integral will split 287 00:20:30,450 --> 00:20:32,180 into four separate integrals? 288 00:20:32,180 --> 00:20:37,170 The total variation will be this difference plus this difference 289 00:20:37,170 --> 00:20:39,350 plus this plus this. 290 00:20:39,350 --> 00:20:42,470 I'll have four-- in this example-- 291 00:20:42,470 --> 00:20:46,960 have four positive contributions to the variation. 292 00:20:51,350 --> 00:20:54,100 And that's the key quantity there. 293 00:20:54,100 --> 00:20:56,470 And what a shock will do-- can I maybe 294 00:20:56,470 --> 00:21:03,630 mention in this figure what a shock does, can do, would 295 00:21:03,630 --> 00:21:10,970 be to somehow cut out part of the variation. 296 00:21:10,970 --> 00:21:14,940 When a shock occurs it might happen that these things are 297 00:21:14,940 --> 00:21:17,410 going into the shock, these are going 298 00:21:17,410 --> 00:21:23,310 in from that side, and the shock kind of jumps 299 00:21:23,310 --> 00:21:28,710 from this value, u_left, to this value, u_right, 300 00:21:28,710 --> 00:21:32,670 following the shock condition when it emerges. 301 00:21:32,670 --> 00:21:37,030 And this additional variation is wiped out. 302 00:21:37,030 --> 00:21:43,020 You see that the variation for the solution that 303 00:21:43,020 --> 00:21:46,810 contains a shock is still-- this part is still there, 304 00:21:46,810 --> 00:21:50,150 but this part is all downhill, so it's just this. 305 00:21:52,710 --> 00:21:54,460 It's that minus that. 306 00:21:54,460 --> 00:21:58,350 And the middle two parts have been cut out. 307 00:21:58,350 --> 00:22:01,840 OK, so that's the total variation in a function. 308 00:22:01,840 --> 00:22:06,560 Let me just take the most important example 309 00:22:06,560 --> 00:22:09,274 that's not total variation diminishing, which 310 00:22:09,274 --> 00:22:09,940 is Lax-Wendroff. 311 00:22:12,720 --> 00:22:16,250 The one-- even in the linear case. 312 00:22:21,330 --> 00:22:27,400 The one experiment we've done with that profile, where, 313 00:22:27,400 --> 00:22:34,220 you remember, we had a wall of water and the exact solution 314 00:22:34,220 --> 00:22:39,820 brought the wall to the left with velocity c, with speed c, 315 00:22:39,820 --> 00:22:44,000 but Lax-Wendroff, which did quite well 316 00:22:44,000 --> 00:22:53,410 staying near the discontinuity, near the jump, 317 00:22:53,410 --> 00:23:00,970 had a little oscillation behind. 318 00:23:00,970 --> 00:23:05,590 That was essentially our one objection to Lax-Wendroff, 319 00:23:05,590 --> 00:23:09,900 was that oscillation showing up. 320 00:23:09,900 --> 00:23:14,970 And so that's telling us that the variation-- see, 321 00:23:14,970 --> 00:23:20,550 the variation in the true solution is just that step. 322 00:23:20,550 --> 00:23:24,520 The variation in the Lax-Wendroff answer 323 00:23:24,520 --> 00:23:28,540 is this bigger step plus this one 324 00:23:28,540 --> 00:23:31,750 plus this one plus this, larger than 1. 325 00:23:36,010 --> 00:23:45,000 Roughly the match is that monotone methods-- so here's 326 00:23:45,000 --> 00:23:46,750 the general picture. 327 00:23:46,750 --> 00:23:57,730 We can get monotone, or I'll say, all positive coefficients, 328 00:23:57,730 --> 00:24:02,330 and those will be total variation diminishing. 329 00:24:02,330 --> 00:24:04,460 Great. 330 00:24:04,460 --> 00:24:07,990 Why don't I draw-- recall a monotone example 331 00:24:07,990 --> 00:24:12,650 was Lax-Friedrichs or upwind. 332 00:24:12,650 --> 00:24:17,320 Upwind will be our model, so let me draw the upwind one. 333 00:24:17,320 --> 00:24:24,830 The upwind one didn't increase the variation, it stayed at 1. 334 00:24:24,830 --> 00:24:31,630 We had no oscillations, but it was low accuracy. 335 00:24:31,630 --> 00:24:32,860 So first order. 336 00:24:38,440 --> 00:24:45,355 We can't do better than first-order accuracy, the order 337 00:24:45,355 --> 00:24:50,860 of accuracy I'm measuring in the smooth part of the solution 338 00:24:50,860 --> 00:24:54,950 because that's where I can take Taylor series and match terms 339 00:24:54,950 --> 00:24:58,730 and see which is the first one that doesn't match. 340 00:24:58,730 --> 00:25:05,950 And Lax-Friedrichs and upwind was first order. 341 00:25:05,950 --> 00:25:22,540 Then our problem is that if I allow a negative coefficient 342 00:25:22,540 --> 00:25:32,080 then variation increases, oscillations will come in, 343 00:25:32,080 --> 00:25:37,730 but I can get the accuracy up to second order or higher. 344 00:25:41,480 --> 00:25:46,130 And Lax-Wendroff is the example, the key example 345 00:25:46,130 --> 00:25:47,960 there that I'll continue with. 346 00:25:47,960 --> 00:25:53,780 OK, that's the idea of total variation. 347 00:25:53,780 --> 00:25:59,650 It just gives us a quantity to work 348 00:25:59,650 --> 00:26:07,040 with when we see oscillation. 349 00:26:07,040 --> 00:26:10,770 And of course, I'm interested in oscillations, 350 00:26:10,770 --> 00:26:14,920 so this is a linear case with a discontinuity. 351 00:26:14,920 --> 00:26:18,540 I don't think it's quite fair to call it a shock. 352 00:26:18,540 --> 00:26:21,500 At least, the shocks that we meet 353 00:26:21,500 --> 00:26:25,230 in conservation laws in the nonlinear problems, 354 00:26:25,230 --> 00:26:28,130 they're created-- the ones we saw here 355 00:26:28,130 --> 00:26:34,970 are created for nonlinear reasons, colliding 356 00:26:34,970 --> 00:26:38,960 characteristics, whereas the linear case has characteristics 357 00:26:38,960 --> 00:26:39,830 all parallel. 358 00:26:42,930 --> 00:26:46,240 So here, we're in the lecture, we're 359 00:26:46,240 --> 00:26:52,460 really into the second lecture into the scientific computing 360 00:26:52,460 --> 00:26:54,010 part of our problem. 361 00:26:54,010 --> 00:26:58,930 So that's the analysis part, the PDE part. 362 00:26:58,930 --> 00:27:08,380 Now comes the experience of quite a lot of people working 363 00:27:08,380 --> 00:27:17,640 hard to identify a way to get second-order accuracy where 364 00:27:17,640 --> 00:27:25,220 it's smooth and monotonicity or something like it 365 00:27:25,220 --> 00:27:28,090 where the shocks are. 366 00:27:28,090 --> 00:27:32,640 So that's essentially the goal. 367 00:27:32,640 --> 00:27:35,850 Our method is going to turn out to be nonlinear. 368 00:27:35,850 --> 00:27:42,510 It's going to have some-- it'll be called a flux limiter. 369 00:27:42,510 --> 00:27:47,680 The key quantity will be the flux. 370 00:27:47,680 --> 00:27:56,940 And for a smooth function when u_right and u_left are 371 00:27:56,940 --> 00:28:02,000 very close, the flux is small, we don't have a problem. 372 00:28:02,000 --> 00:28:06,250 But at a shock where there's a significant jump in u 373 00:28:06,250 --> 00:28:14,250 and especially in f of u, the flux limiter 374 00:28:14,250 --> 00:28:20,160 is going to come in, into the discrete method. 375 00:28:20,160 --> 00:28:23,400 Let me begin discrete methods with a new word, 376 00:28:23,400 --> 00:28:24,520 finite volumes. 377 00:28:27,650 --> 00:28:34,320 And with a new interpretation of capital U. 378 00:28:34,320 --> 00:28:40,950 So the discrete value, capital U_(j, n), which, 379 00:28:40,950 --> 00:28:44,300 up to now we've been thinking of as that's an approximation 380 00:28:44,300 --> 00:28:48,860 to the true solution, at j delta x, n delta t. 381 00:28:48,860 --> 00:28:52,290 Now, because we're dealing with discontinuities, 382 00:28:52,290 --> 00:28:54,060 we better integrate. 383 00:28:54,060 --> 00:28:58,230 We integrate over a delta x interval. 384 00:28:58,230 --> 00:29:01,130 We take the average over a little delta x interval 385 00:29:01,130 --> 00:29:04,120 instead of hitting just the point. 386 00:29:04,120 --> 00:29:08,270 So the delta x interval from j minus 1/2 to j plus 1/2, 387 00:29:08,270 --> 00:29:13,340 we take the average of u at the time n delta t 388 00:29:13,340 --> 00:29:23,670 and that's what finite volume-- and let's think of as capital 389 00:29:23,670 --> 00:29:27,360 U, the discrete approximation. 390 00:29:27,360 --> 00:29:31,940 The appropriateness of that is in the integral form 391 00:29:31,940 --> 00:29:34,260 of the conservation law. 392 00:29:34,260 --> 00:29:36,330 So the integral form of a conservation law 393 00:29:36,330 --> 00:29:49,670 leads us to that natural discretization, d by dt of U, 394 00:29:49,670 --> 00:29:52,060 of the average. 395 00:29:52,060 --> 00:29:55,750 You're keeping your eye on the integral form 396 00:29:55,750 --> 00:29:56,960 of the conservation law. 397 00:29:56,960 --> 00:29:59,720 I'll just come back to it again. 398 00:29:59,720 --> 00:30:04,280 The derivative of the mass, the integral 399 00:30:04,280 --> 00:30:10,430 of u-- so in the discrete case, the difference of capital U-- 400 00:30:10,430 --> 00:30:15,560 plus the integrated term of flux at one end 401 00:30:15,560 --> 00:30:17,530 of the interval minus a flux at the other end. 402 00:30:17,530 --> 00:30:22,410 And notice that we are getting for the flux, half-- 403 00:30:22,410 --> 00:30:32,320 so a staggered grid is presenting itself naturally. 404 00:30:32,320 --> 00:30:35,780 And so that f_(j+1/2), the fluxes, 405 00:30:35,780 --> 00:30:41,130 which is the critical thing, is on the j plus 1/2 grid. 406 00:30:41,130 --> 00:30:45,650 All right, now this could be exactly defined 407 00:30:45,650 --> 00:30:52,200 as the flux at that grid point. 408 00:30:52,200 --> 00:30:55,990 Capital U squared over 2 in the Burger equation, 409 00:30:55,990 --> 00:30:58,550 at that grid point. 410 00:30:58,550 --> 00:31:04,560 That would give us a totally straightforward, first-order 411 00:31:04,560 --> 00:31:05,800 numerical method. 412 00:31:05,800 --> 00:31:08,710 First order, I see first order here in time 413 00:31:08,710 --> 00:31:11,550 and it might be better than that in space 414 00:31:11,550 --> 00:31:13,250 because in some way that's centered. 415 00:31:15,970 --> 00:31:19,450 But if we want higher accuracy, if we 416 00:31:19,450 --> 00:31:21,700 want to include Lax-Wendroff, for example, 417 00:31:21,700 --> 00:31:26,530 that includes several values of u, 418 00:31:26,530 --> 00:31:31,510 we go to a more general flux function. 419 00:31:34,670 --> 00:31:39,610 The whole job here is create a smart flux function. 420 00:31:39,610 --> 00:31:48,790 So a flux function depends not only on u in the center, u_j, 421 00:31:48,790 --> 00:31:52,100 at time n-- this is all at time n-- 422 00:31:52,100 --> 00:31:56,350 it can also depend on values to the left of it 423 00:31:56,350 --> 00:31:57,400 and values to the right. 424 00:32:00,350 --> 00:32:02,970 It may use those values, it may not, 425 00:32:02,970 --> 00:32:05,870 but we allow that possibility then, 426 00:32:05,870 --> 00:32:08,900 that the flux-- so this is the numerical method 427 00:32:08,900 --> 00:32:10,440 that you aim for. 428 00:32:10,440 --> 00:32:16,210 Is to choose a flux function and of course, 429 00:32:16,210 --> 00:32:19,620 this should be consistent with the original problem 430 00:32:19,620 --> 00:32:26,180 and that simply says that if all the u's were the same, 431 00:32:26,180 --> 00:32:32,159 that the numerical flux function would give you the flux 432 00:32:32,159 --> 00:32:33,450 function in the given equation. 433 00:32:36,860 --> 00:32:39,750 OK, so this is our general pattern. 434 00:32:39,750 --> 00:32:43,780 And now, the question is really, maybe the first step 435 00:32:43,780 --> 00:33:00,900 is write a-- I would like to write upwind and Lax-Wendroff-- 436 00:33:00,900 --> 00:33:02,960 I'd like to find their flux functions. 437 00:33:02,960 --> 00:33:06,710 I guess I'm saying that when we look 438 00:33:06,710 --> 00:33:14,630 at an equation in this format, we 439 00:33:14,630 --> 00:33:17,150 can take methods that we have already used 440 00:33:17,150 --> 00:33:22,000 and see, OK, even in the linear case, the one-way wave 441 00:33:22,000 --> 00:33:26,850 equation, what are we picking for the flux function? 442 00:33:26,850 --> 00:33:32,170 OK, now the one distinction is: the one-way wave equation, 443 00:33:32,170 --> 00:33:38,150 we had a constant, and upwind, for example, 444 00:33:38,150 --> 00:33:44,070 went in the direction decided by that speed, c. 445 00:33:44,070 --> 00:33:49,350 If the speed reversed sign so that the waves were going 446 00:33:49,350 --> 00:33:52,390 the other way, the wind is coming the other way, 447 00:33:52,390 --> 00:33:55,360 the upwind direction is the other way, 448 00:33:55,360 --> 00:33:57,870 so our method would change. 449 00:34:02,770 --> 00:34:05,450 Let me just write down the flux function 450 00:34:05,450 --> 00:34:12,210 that turns out to come for upwind 451 00:34:12,210 --> 00:34:14,290 and then for Lax-Wendroff. 452 00:34:14,290 --> 00:34:18,120 OK, so for upwind the flux turns out to be-- 453 00:34:18,120 --> 00:34:24,730 F turns out to be 1/2 of-- because it's at this, 454 00:34:24,730 --> 00:34:26,720 this is f_(j+1/2), really. 455 00:34:33,460 --> 00:34:34,510 It's at the center point. 456 00:34:34,510 --> 00:34:44,000 The flux function will be-- I'll put in the a-- a delta t-- 457 00:34:44,000 --> 00:34:52,620 sorry, a over 2 U_(j+1) plus U_j, now, plus-- 458 00:34:52,620 --> 00:35:05,290 or rather minus actually-- comes absolute value of a over 2 459 00:35:05,290 --> 00:35:08,710 times u_(j+1) minus u_j. 460 00:35:12,510 --> 00:35:19,660 If we're patient we can track down the two cases, a positive 461 00:35:19,660 --> 00:35:28,770 and a negative, and recognize that this form switches 462 00:35:28,770 --> 00:35:35,160 direction correctly on the difference method. 463 00:35:35,160 --> 00:35:37,230 So it keeps us upwind, in other words. 464 00:35:37,230 --> 00:35:39,750 If the wind changes direction our flux 465 00:35:39,750 --> 00:35:44,390 changes appropriately to follow that direction. 466 00:35:44,390 --> 00:35:48,720 So that's an OK flux function. 467 00:35:48,720 --> 00:35:57,480 It's a monotone one and it's first order. 468 00:35:57,480 --> 00:36:00,190 Now let me put in Lax-Wendroff. 469 00:36:00,190 --> 00:36:02,420 What would Lax-Wendroff be? 470 00:36:02,420 --> 00:36:08,790 Would be the same, the same upwind flux, if I call that UP 471 00:36:08,790 --> 00:36:12,310 for upwind and then a correction, 472 00:36:12,310 --> 00:36:19,670 a Lax-Wendroff term that gets the second-order accuracy. 473 00:36:19,670 --> 00:36:26,730 OK, and that turns out to be-- it turns out that I multiply 474 00:36:26,730 --> 00:36:32,780 this term by that ratio r. 475 00:36:32,780 --> 00:36:34,340 By the absolute value of r. 476 00:36:34,340 --> 00:36:39,700 So it turns out that it's-- the Lax-Wendroff thing, 477 00:36:39,700 --> 00:36:42,240 to get that extra accuracy, you might remember, 478 00:36:42,240 --> 00:36:45,200 had an extra delta t over delta x. 479 00:36:45,200 --> 00:36:51,150 And so it's this same thing times a over 2. 480 00:36:51,150 --> 00:36:55,730 In the nonlinear case it'll be f of u_(j+1). 481 00:36:59,050 --> 00:37:03,200 f at u_(j+1) minus f at u_j. 482 00:37:10,090 --> 00:37:16,640 Next time, I'll come back to this point 483 00:37:16,640 --> 00:37:24,040 and pin down this more carefully. 484 00:37:24,040 --> 00:37:26,450 All I want to do in completing today 485 00:37:26,450 --> 00:37:35,210 is introduce this idea of a flux limiter. 486 00:37:35,210 --> 00:37:42,920 So what people say when they look at Lax-Wendroff 487 00:37:42,920 --> 00:37:49,010 is that this term that's been added on 488 00:37:49,010 --> 00:37:52,280 is fine in the case of a smooth solution, 489 00:37:52,280 --> 00:37:59,320 but in that case of a jump, like our picture 490 00:37:59,320 --> 00:38:07,030 here, that term is increasing variation. 491 00:38:07,030 --> 00:38:09,650 It's an anti-diffusion term. 492 00:38:09,650 --> 00:38:14,390 Instead of smoothing things out the way F upwind does, 493 00:38:14,390 --> 00:38:21,830 that extra term has the effect of increasing the jump. 494 00:38:21,830 --> 00:38:23,740 It's anti-diffusive. 495 00:38:23,740 --> 00:38:31,500 And therefore, when it's a significant term, 496 00:38:31,500 --> 00:38:39,120 the flux limiter cuts back on it. 497 00:38:39,120 --> 00:38:41,070 The result will be a numerical method, 498 00:38:41,070 --> 00:38:43,580 which is not linear anymore even if we apply it 499 00:38:43,580 --> 00:38:45,300 to a linear problem. 500 00:38:45,300 --> 00:38:49,190 The result won't be strictly linear 501 00:38:49,190 --> 00:38:54,400 and it will produce a profile that 502 00:38:54,400 --> 00:38:56,960 doesn't have that oscillation. 503 00:38:56,960 --> 00:39:00,630 So it has what we want. 504 00:39:00,630 --> 00:39:05,410 It gives us good answers when the solution is smooth, 505 00:39:05,410 --> 00:39:08,820 it stays near the shock the way Lax-Wendroff does, 506 00:39:08,820 --> 00:39:11,490 but it doesn't oscillate. 507 00:39:11,490 --> 00:39:15,550 OK, so what's the flux limiter? 508 00:39:15,550 --> 00:39:20,000 The flux limiter is, I multiply in here by some function, 509 00:39:20,000 --> 00:39:26,590 the limiter function, phi_(j+1) and I multiply in here 510 00:39:26,590 --> 00:39:28,230 by phi_j. 511 00:39:28,230 --> 00:39:33,050 In other words, there's an extra factor, 512 00:39:33,050 --> 00:39:36,500 this flux limiting factor. 513 00:39:36,500 --> 00:39:40,090 So the question is, how to construct that? 514 00:39:40,090 --> 00:39:43,700 What formula do we use for the flux limiter? 515 00:39:43,700 --> 00:39:49,160 So the flux limiter, after lot of experiments 516 00:39:49,160 --> 00:40:00,540 is chosen to depend on-- so we want it to be 1. 517 00:40:00,540 --> 00:40:02,980 So what do we want? 518 00:40:02,980 --> 00:40:10,380 We want it to be 1 or close to 1 at smooth part 519 00:40:10,380 --> 00:40:17,480 and we want it to be greater than 1 at a discontinuity. 520 00:40:20,080 --> 00:40:23,030 OK, so our question is, how do you 521 00:40:23,030 --> 00:40:25,420 identify that discontinuity? 522 00:40:25,420 --> 00:40:29,850 Well, the convenient way has turned out 523 00:40:29,850 --> 00:40:35,980 to be look at the ratio-- look at differences. 524 00:40:35,980 --> 00:40:40,000 Look at the ratio of-- so what comes in here 525 00:40:40,000 --> 00:40:46,160 will be the ratio r, say, j plus 1/2. 526 00:40:46,160 --> 00:40:53,610 That'll be the ratio of U, say j plus 2 minus U_(j+1) 1 527 00:40:53,610 --> 00:40:58,090 to U_(j+1) minus U_j. 528 00:40:58,090 --> 00:41:00,910 All at time n. 529 00:41:00,910 --> 00:41:04,250 That ratio tells us-- that ratio is close to 1, 530 00:41:04,250 --> 00:41:16,870 if-- so the ratio r is close to 1 for smooth and our function, 531 00:41:16,870 --> 00:41:23,620 phi of r then, we want our function phi, when the ratio is 532 00:41:23,620 --> 00:41:26,630 1, to give the value 1. 533 00:41:26,630 --> 00:41:29,270 So if I try to draw a graph of-- let 534 00:41:29,270 --> 00:41:43,960 me draw a graph of people's flux functions-- of the phi of r. 535 00:41:43,960 --> 00:41:47,180 So here's a graph of phi of r. 536 00:41:47,180 --> 00:41:50,230 First on the graph I'll put upwind and Lax-Wendroff 537 00:41:50,230 --> 00:41:51,890 without anything. 538 00:41:51,890 --> 00:41:59,640 So Lax-Wendroff would take-- phi Lax-Wendroff is identically 1. 539 00:41:59,640 --> 00:42:02,660 This is the ratio r and this is phi. 540 00:42:06,490 --> 00:42:15,250 So Lax-Wendroff, before fixing, had no flux limiter function 541 00:42:15,250 --> 00:42:17,100 and it just had 1. 542 00:42:17,100 --> 00:42:22,680 And actually, if I can keep the same formula, 543 00:42:22,680 --> 00:42:27,400 phi_upwind would just be zero. 544 00:42:27,400 --> 00:42:34,280 Then it'll turn out that-- let me draw these ratios. 545 00:42:34,280 --> 00:42:45,750 There is a region here and here where the flux limiter 546 00:42:45,750 --> 00:42:47,390 is doing the right thing. 547 00:42:47,390 --> 00:42:53,710 So in that region it's acting the right way. 548 00:42:53,710 --> 00:42:57,350 It's limiting the anti-diffusive flux. 549 00:42:57,350 --> 00:42:59,670 It's controlling the variation. 550 00:42:59,670 --> 00:43:01,710 It's controlling the oscillation. 551 00:43:01,710 --> 00:43:11,020 So various people have proposed-- let's see. 552 00:43:11,020 --> 00:43:12,760 No, I haven't got this right. 553 00:43:16,180 --> 00:43:20,330 I think it's-- because I know that part is wrong 554 00:43:20,330 --> 00:43:21,290 for Lax-Wendroff. 555 00:43:21,290 --> 00:43:22,270 So maybe it's there. 556 00:43:25,200 --> 00:43:28,790 No, I've forgotten-- sorry that I'm 557 00:43:28,790 --> 00:43:41,210 not-- I've lost this because I know that one possibility is 558 00:43:41,210 --> 00:43:44,850 this as the flux limiter. 559 00:43:44,850 --> 00:43:48,850 One possibility would be this one 560 00:43:48,850 --> 00:43:51,600 and then there's another possibility that 561 00:43:51,600 --> 00:43:54,440 is OK that goes through here. 562 00:43:59,950 --> 00:44:06,770 Let me stop for today with this idea of a flux limiter 563 00:44:06,770 --> 00:44:12,380 and you see that it's a task to fix the method 564 00:44:12,380 --> 00:44:15,650 to do the right thing at smooth parts 565 00:44:15,650 --> 00:44:18,090 and also the right thing at discontinuities. 566 00:44:18,090 --> 00:44:20,130 OK, I'll pick up on that next time. 567 00:44:20,130 --> 00:44:21,380 Thanks.