1 00:00:00,000 --> 00:00:02,320 The following content is provided 2 00:00:02,320 --> 00:00:04,910 by MIT OpenCourseWare under a Creative Commons license. 3 00:00:04,910 --> 00:00:09,510 Additional information about our license and MIT OpenCourseWare, 4 00:00:09,510 --> 00:00:11,770 in general, is available at ocw.mit.edu. 5 00:00:16,950 --> 00:00:19,020 PROFESSOR: OK. 6 00:00:19,020 --> 00:00:27,170 Last time, I started to talk about the heat equation, also 7 00:00:27,170 --> 00:00:29,830 known as the diffusion equation. 8 00:00:29,830 --> 00:00:36,790 And what we did then was to solve it. 9 00:00:36,790 --> 00:00:39,540 Give an analytical solution. 10 00:00:39,540 --> 00:00:44,180 Both for any initial value u of x and 0, 11 00:00:44,180 --> 00:00:49,540 and then, specifically for the point source. 12 00:00:49,540 --> 00:00:52,290 So then, we're talking about the impulse response, 13 00:00:52,290 --> 00:00:55,160 or the fundamental solution, when 14 00:00:55,160 --> 00:00:56,960 we start from a delta function. 15 00:00:56,960 --> 00:01:01,090 And, it was that Gaussian, bell-shaped curve. 16 00:01:01,090 --> 00:01:04,710 Fantastic fundamental solution. 17 00:01:04,710 --> 00:01:09,270 But the formulas you get involving 18 00:01:09,270 --> 00:01:12,700 infinite integrals are not the greatest for actually getting 19 00:01:12,700 --> 00:01:17,560 numbers for the whole temperature 20 00:01:17,560 --> 00:01:21,930 history, u of t and x. 21 00:01:21,930 --> 00:01:25,120 So we have to go to finite differences. 22 00:01:25,120 --> 00:01:28,600 And, of course, we need finite differences also 23 00:01:28,600 --> 00:01:38,500 in more general cases, where non-linear terms could appear. 24 00:01:38,500 --> 00:01:45,260 Or even just variation in x and t in the linear coefficient, 25 00:01:45,260 --> 00:01:48,360 which I had just set to be 1. 26 00:01:48,360 --> 00:01:50,520 So, this is the first job. 27 00:01:50,520 --> 00:01:55,020 Difference methods for the heat equation. 28 00:01:55,020 --> 00:02:02,590 And you'll see that we get pushed toward implicit methods. 29 00:02:02,590 --> 00:02:06,400 Because explicit method will require delta t 30 00:02:06,400 --> 00:02:11,130 to be that very small size delta x squared, 31 00:02:11,130 --> 00:02:15,130 and that's pretty slow-going. 32 00:02:15,130 --> 00:02:22,040 And of course, what I'm saying applies equally to -- 33 00:02:22,040 --> 00:02:28,950 we might be in 2D or in 3D diffusion of pollution, 34 00:02:28,950 --> 00:02:33,850 for example, in environmental engineering. 35 00:02:33,850 --> 00:02:36,620 Our fundamental solution from the delta function 36 00:02:36,620 --> 00:02:42,420 would easily extend to two dimensions or three. 37 00:02:42,420 --> 00:02:44,820 And, our difference methods will extend, 38 00:02:44,820 --> 00:02:47,650 but we'll see what comes up. 39 00:02:47,650 --> 00:02:53,540 And then, let me point to this second problem which 40 00:02:53,540 --> 00:02:55,360 combines the two. 41 00:02:55,360 --> 00:02:58,670 We already know about one-way waves, 42 00:02:58,670 --> 00:03:01,300 coming from a convection term. 43 00:03:01,300 --> 00:03:07,880 Now we're getting the basic ideas for the diffusion term 44 00:03:07,880 --> 00:03:14,330 and lot of very important, real problems combine the two. 45 00:03:14,330 --> 00:03:20,480 And then you have a situation that I'm just 46 00:03:20,480 --> 00:03:21,490 sort of pointing to. 47 00:03:21,490 --> 00:03:26,020 Here, that coefficient, as before, 48 00:03:26,020 --> 00:03:29,400 has the units of velocity. 49 00:03:29,400 --> 00:03:30,660 Distance over time. 50 00:03:30,660 --> 00:03:32,600 You see, if I bring time up here, 51 00:03:32,600 --> 00:03:34,880 I have time over distance there. 52 00:03:34,880 --> 00:03:38,600 So I need distance over time from c. 53 00:03:38,600 --> 00:03:42,520 It's what we've seen always. 54 00:03:42,520 --> 00:03:47,950 And, then the c delta t over delta x, the quantity r, 55 00:03:47,950 --> 00:03:49,050 was dimensionless. 56 00:03:49,050 --> 00:03:54,670 So, I'm taking three minutes to speak about dimensions. 57 00:03:54,670 --> 00:03:57,890 Dimensional analysis is a very simple but useful thing 58 00:03:57,890 --> 00:03:59,990 to do at the beginning of a problem. 59 00:03:59,990 --> 00:04:05,290 So here we had r equals c delta t over delta x, which 60 00:04:05,290 --> 00:04:07,000 is dimensionless. 61 00:04:07,000 --> 00:04:12,490 And, typically for explicit methods, 62 00:04:12,490 --> 00:04:16,950 it had to be bounded by 1, for example. 63 00:04:16,950 --> 00:04:23,570 Now, for the diffusion term d. 64 00:04:23,570 --> 00:04:25,870 What are the units of d? 65 00:04:25,870 --> 00:04:29,550 Well, here, now I have distance squared in the denominator, 66 00:04:29,550 --> 00:04:35,070 so d needs to be a distance squared over time. 67 00:04:35,070 --> 00:04:43,870 So it has different dimensions, and the quantity capital 68 00:04:43,870 --> 00:04:47,050 R that we'll meet, parallel to this, 69 00:04:47,050 --> 00:04:53,890 will be d delta t over delta x squared. 70 00:04:53,890 --> 00:04:58,790 And, again, it's dimensionless. 71 00:04:58,790 --> 00:05:02,820 Because if d has dimensions of distance squared over time, 72 00:05:02,820 --> 00:05:06,200 that reverses it to give dimensionless quantity, 73 00:05:06,200 --> 00:05:09,680 and we'll see that explicit methods require 74 00:05:09,680 --> 00:05:16,660 a bound on capital R, and that's what so forces us 75 00:05:16,660 --> 00:05:18,711 into that small time step. 76 00:05:18,711 --> 00:05:19,210 OK. 77 00:05:19,210 --> 00:05:23,900 So you'll see capital R playing a similar role to little r. 78 00:05:23,900 --> 00:05:26,120 And actually, while we're talking 79 00:05:26,120 --> 00:05:32,510 about dimensional analysis, back to the differential equation, 80 00:05:32,510 --> 00:05:35,900 suppose we want, as we certainly do, 81 00:05:35,900 --> 00:05:40,670 to compare the importance of the convection 82 00:05:40,670 --> 00:05:45,250 term and the diffusion term. 83 00:05:45,250 --> 00:05:55,130 Well, that comparison is somehow the ratio of c to d, So, 84 00:05:55,130 --> 00:05:58,540 somehow, it's natural to think of the ratio of c 85 00:05:58,540 --> 00:06:06,440 to d as telling us a very important fact. 86 00:06:06,440 --> 00:06:10,570 Is our problem more like a wave problem? 87 00:06:10,570 --> 00:06:20,220 Do we expect to be on the unit circle, or close, 88 00:06:20,220 --> 00:06:25,620 with a term that in Fourier space is pure imaginary, 89 00:06:25,620 --> 00:06:27,680 from the i*k? 90 00:06:27,680 --> 00:06:35,730 Or, if d is big, we have then a small ratio, 91 00:06:35,730 --> 00:06:38,870 and that tells us that the diffusion term is important, 92 00:06:38,870 --> 00:06:41,650 that we've got lots of dissipation, 93 00:06:41,650 --> 00:06:46,740 that the stability, at least in the differential equation, 94 00:06:46,740 --> 00:06:52,640 is much greater, because this is producing a minus k squared 95 00:06:52,640 --> 00:06:57,720 in Fourier space, and a strongly negative term. 96 00:06:57,720 --> 00:07:01,260 So it's that ratio c to d but of course, my point 97 00:07:01,260 --> 00:07:06,500 here in the dimensional stuff is that c to d 98 00:07:06,500 --> 00:07:12,220 isn't quite OK as it is, because they have different units. 99 00:07:12,220 --> 00:07:14,410 And we need one more distance. 100 00:07:14,410 --> 00:07:19,400 So we need the ratio of c -- where do I need an L? 101 00:07:19,400 --> 00:07:21,390 I need another L up there. 102 00:07:21,390 --> 00:07:30,210 I need another distance there to get something that 103 00:07:30,210 --> 00:07:33,620 is altogether dimensionless. 104 00:07:36,150 --> 00:07:39,980 And so that's the critical number. 105 00:07:39,980 --> 00:07:42,930 And, it's a known as the Peclet Number. 106 00:07:42,930 --> 00:07:46,910 Pe, I'll call it. 107 00:07:46,910 --> 00:07:49,160 What is L? 108 00:07:49,160 --> 00:07:53,740 We're talking here only about the differential equation. l is 109 00:07:53,740 --> 00:07:55,580 characteristic length. 110 00:07:55,580 --> 00:08:01,110 If our problem is set on a finite interval, 111 00:08:01,110 --> 00:08:03,430 then it's probably the length of the interval. 112 00:08:03,430 --> 00:08:06,540 It's a characteristic length of the problem. 113 00:08:06,540 --> 00:08:16,340 And, then this, quantity c*L over d gives us a good measure 114 00:08:16,340 --> 00:08:18,800 of the importance of these two. 115 00:08:18,800 --> 00:08:23,580 And, we'll see in the finite difference case, 116 00:08:23,580 --> 00:08:26,870 there will be cell Peclet number, 117 00:08:26,870 --> 00:08:28,360 which is entirely different. 118 00:08:28,360 --> 00:08:33,740 A Peclet number for the little finite different cell, 119 00:08:33,740 --> 00:08:39,440 where this L is changed to like a delta x. 120 00:08:39,440 --> 00:08:41,600 OK. 121 00:08:41,600 --> 00:08:47,640 I think we've got it there. 122 00:08:47,640 --> 00:08:52,470 Have I got the numerator and denominator correct there, 123 00:08:52,470 --> 00:08:54,760 or should they be reversed? 124 00:08:54,760 --> 00:08:56,500 Is it time over distance? 125 00:09:05,210 --> 00:09:06,410 It is, isn't it? 126 00:09:06,410 --> 00:09:07,510 Time over distance. 127 00:09:07,510 --> 00:09:10,500 That shows that I'm not a dimensional analyst. 128 00:09:10,500 --> 00:09:12,740 Let me fix it. 129 00:09:12,740 --> 00:09:16,110 And, then that will mean fixing this. 130 00:09:16,110 --> 00:09:22,310 So, let me just fix that. 131 00:09:22,310 --> 00:09:25,220 It doesn't sound right as I said it, but I read it wrong. 132 00:09:25,220 --> 00:09:25,870 OK. 133 00:09:25,870 --> 00:09:28,770 So c is a time over distance. 134 00:09:28,770 --> 00:09:30,700 Thanks. 135 00:09:30,700 --> 00:09:33,400 Is that right? 136 00:09:33,400 --> 00:09:35,400 It was right the first time? 137 00:09:35,400 --> 00:09:38,800 In that case, we'll move on. 138 00:09:38,800 --> 00:09:40,270 OK. 139 00:09:40,270 --> 00:09:41,240 All right. 140 00:09:41,240 --> 00:09:42,080 OK. 141 00:09:42,080 --> 00:09:43,070 Good. 142 00:09:43,070 --> 00:09:45,190 Thanks. 143 00:09:45,190 --> 00:09:45,900 Yeah. 144 00:09:45,900 --> 00:09:47,220 You're right. 145 00:09:47,220 --> 00:09:49,800 OK. 146 00:09:49,800 --> 00:09:51,210 Sorry that's on the tapes. 147 00:09:51,210 --> 00:09:51,710 OK. 148 00:09:54,260 --> 00:09:58,890 So now let's begin with a natural explicit method 149 00:09:58,890 --> 00:10:00,620 for the heat equation. 150 00:10:00,620 --> 00:10:02,660 What everybody would think of. 151 00:10:02,660 --> 00:10:08,270 So, this will be explicit for the heat equations, 152 00:10:08,270 --> 00:10:11,080 and I'm in 1D. 153 00:10:11,080 --> 00:10:13,670 And, what do you think of? 154 00:10:13,670 --> 00:10:17,090 Naturally, forward difference in time. 155 00:10:17,090 --> 00:10:19,120 So, can I write it this way? 156 00:10:19,120 --> 00:10:23,580 Time difference -- I'll use capital U, as always, 157 00:10:23,580 --> 00:10:28,650 for the finite difference solution -- divided by delta x, 158 00:10:28,650 --> 00:10:35,510 and I'm doing the heat equation with c equal to 1. 159 00:10:35,510 --> 00:10:38,270 Just to keep things simple, let's just 160 00:10:38,270 --> 00:10:43,470 do the heat equation normalized. 161 00:10:43,470 --> 00:10:47,090 So, here I'm going to take a second difference 162 00:10:47,090 --> 00:10:51,950 in the x direction of u divided by delta x squared. 163 00:10:54,630 --> 00:10:55,130 OK. 164 00:10:58,490 --> 00:10:59,780 What's the accuracy? 165 00:11:03,000 --> 00:11:06,750 Well, we know the accuracy of each term, 166 00:11:06,750 --> 00:11:11,290 the centered difference, the error is the discretization 167 00:11:11,290 --> 00:11:20,040 error, the local error, the truncation error, 168 00:11:20,040 --> 00:11:25,620 it's also called, will be proportional to -- oh, 169 00:11:25,620 --> 00:11:30,420 I'm sorry, that's delta t of course. 170 00:11:30,420 --> 00:11:32,360 OK. 171 00:11:32,360 --> 00:11:37,010 So, the discretization error in a forward difference, 172 00:11:37,010 --> 00:11:50,400 when I divide by delta t -- so it's first order in t, 173 00:11:50,400 --> 00:11:55,040 and centered difference because of the symmetry is second order 174 00:11:55,040 --> 00:11:55,540 in x. 175 00:11:58,470 --> 00:11:58,970 OK. 176 00:12:01,580 --> 00:12:03,640 And this would be another problem 177 00:12:03,640 --> 00:12:09,830 where I could have begun with a semi-discretization. 178 00:12:09,830 --> 00:12:13,950 I could have begun with keeping this 179 00:12:13,950 --> 00:12:19,450 as a derivative, du/dt, equal this centered difference, 180 00:12:19,450 --> 00:12:25,210 and the discussion would have been entirely parallel to what 181 00:12:25,210 --> 00:12:27,700 we'll have here. 182 00:12:27,700 --> 00:12:32,560 So, I've jumped right to the finite difference in x and t. 183 00:12:32,560 --> 00:12:35,110 So, in other words, what's happening here? 184 00:12:35,110 --> 00:12:40,730 U_(n+1), U at that center point, n plus 1, of course, 185 00:12:40,730 --> 00:12:46,730 our little molecule is computing this backward difference, 186 00:12:46,730 --> 00:12:51,490 and this -- or rather forward in time, 187 00:12:51,490 --> 00:12:55,310 so this is the time direction, n plus 1 -- 188 00:12:55,310 --> 00:12:57,550 and this centered difference in space. 189 00:12:57,550 --> 00:13:03,580 So if I just write out what that does, it's u_(j, 190 00:13:03,580 --> 00:13:08,950 n) coming from the time difference. 191 00:13:08,950 --> 00:13:13,330 Then, I'll multiply by delta t, and I have a delta x squared, 192 00:13:13,330 --> 00:13:18,390 so that ratio is what I'm calling capital R. Capital R, 193 00:13:18,390 --> 00:13:28,100 to distinguish the diffusion ratio from the wave ratio. 194 00:13:28,100 --> 00:13:29,170 Times what? 195 00:13:29,170 --> 00:13:36,570 So it's just going to be U_(j+1) at time n minus 2*U_j at time 196 00:13:36,570 --> 00:13:41,480 n, plus U_(j-1) at time n. 197 00:13:41,480 --> 00:13:45,430 We're good at writing that expression down. 198 00:13:45,430 --> 00:13:50,440 We can follow e to the i*k*x, of course. 199 00:13:50,440 --> 00:13:52,920 That's always the idea. 200 00:13:52,920 --> 00:14:01,550 Start with U_(j, 0) as e to the i*k*j delta x, 201 00:14:01,550 --> 00:14:03,990 and follow it up. 202 00:14:03,990 --> 00:14:14,200 We're looking for a growth factor G, 203 00:14:14,200 --> 00:14:16,340 the growth factor in a single step, 204 00:14:16,340 --> 00:14:20,750 and, of course, the point is, it depends on k. 205 00:14:24,800 --> 00:14:30,460 At at every step, the exponential is always with us. 206 00:14:30,460 --> 00:14:35,260 It's the factor that multiplies that exponential that we want, 207 00:14:35,260 --> 00:14:37,980 and you see what it is, don't you? 208 00:14:37,980 --> 00:14:41,620 Can you see what G is here? 209 00:14:41,620 --> 00:14:47,650 So, if I put in that, and then ultimately factor out 210 00:14:47,650 --> 00:14:51,290 the exponential, it'll be multiplied by a 1, 211 00:14:51,290 --> 00:14:56,980 from that part, and an R from this, and now in parentheses, 212 00:14:56,980 --> 00:15:01,995 we have the now very familiar quantity that comes from e 213 00:15:01,995 --> 00:15:09,270 to the i*k delta x, e to the minus i*k delta x. 214 00:15:09,270 --> 00:15:14,650 Those combine to give 2 cosine k delta x, 215 00:15:14,650 --> 00:15:18,590 and the center one is minus 2. 216 00:15:18,590 --> 00:15:20,060 OK. 217 00:15:20,060 --> 00:15:22,690 That's our growth factor. 218 00:15:22,690 --> 00:15:26,180 And the question is, does it stay 219 00:15:26,180 --> 00:15:32,910 below 1, which is the stability requirement, for all k delta x. 220 00:15:32,910 --> 00:15:37,860 And, so you can see the dangerous frequency as, 221 00:15:37,860 --> 00:15:44,030 often, the dangerous frequency is when k delta x is pi, 222 00:15:44,030 --> 00:15:47,820 and that cosine becomes minus 1. 223 00:15:47,820 --> 00:15:59,980 And, then we have -- So, stability at k delta x equal pi 224 00:15:59,980 --> 00:16:06,820 needs -- so this becomes 1 and this is negative 2, 225 00:16:06,820 --> 00:16:07,890 and that's negative 2. 226 00:16:07,890 --> 00:16:13,840 We have 1 minus 4R, and that quantity 227 00:16:13,840 --> 00:16:18,960 has to be, in absolute value, not exceeding 1. 228 00:16:18,960 --> 00:16:23,800 Well, there's no danger that it's going to go above 1. 229 00:16:23,800 --> 00:16:27,370 The danger is, will it go below minus 1, 230 00:16:27,370 --> 00:16:30,070 and it's not allowed to. 231 00:16:30,070 --> 00:16:34,670 So we need it to be greater than or equal to minus 1. 232 00:16:34,670 --> 00:16:37,910 And, now, if I put the minus 1 on this side and the 4R 233 00:16:37,910 --> 00:16:42,890 on that side, that means 2 is greater or equal to 4R. 234 00:16:42,890 --> 00:16:49,690 And, I get the famous condition R less or equal to half. 235 00:16:49,690 --> 00:16:52,300 That's the stability requirement. 236 00:16:55,090 --> 00:16:58,880 And if there was a coefficient d in the heat equation, 237 00:16:58,880 --> 00:17:03,620 the coefficient would be there in R. OK. 238 00:17:03,620 --> 00:17:05,440 So, that's the stability condition, 239 00:17:05,440 --> 00:17:12,710 and that's the thing we would like to find a way around, 240 00:17:12,710 --> 00:17:17,120 because, it says that this delta t -- 241 00:17:17,120 --> 00:17:20,390 I'll just put that in, less or equal a half, 242 00:17:20,390 --> 00:17:25,120 and that means a very small delta t. 243 00:17:25,120 --> 00:17:26,910 OK. 244 00:17:26,910 --> 00:17:30,220 What's the way around it? 245 00:17:30,220 --> 00:17:33,220 To go to an explicit method. 246 00:17:33,220 --> 00:17:38,660 So, let me go first to a completely explicit method. 247 00:17:38,660 --> 00:17:39,680 Fully explicit. 248 00:17:42,390 --> 00:17:43,390 I'm sorry. 249 00:17:43,390 --> 00:17:46,140 To go to an implicit method. 250 00:17:46,140 --> 00:17:54,780 Now, implicit method -- now I'm going to do -- 251 00:17:54,780 --> 00:17:59,090 switch to an implicit method by computing this difference 252 00:17:59,090 --> 00:18:01,120 at time n plus 1. 253 00:18:01,120 --> 00:18:03,100 So this is now n plus 1. 254 00:18:05,890 --> 00:18:09,410 And, what's the change in G? 255 00:18:09,410 --> 00:18:10,861 A big change. 256 00:18:10,861 --> 00:18:11,360 OK. 257 00:18:11,360 --> 00:18:19,960 So now I want to recompute G, for the growth 258 00:18:19,960 --> 00:18:26,290 factor over a single step for the exponential, when 259 00:18:26,290 --> 00:18:29,805 the right-hand side, the second difference 260 00:18:29,805 --> 00:18:31,930 is taken at the new time. 261 00:18:31,930 --> 00:18:38,130 So this has to move over to the other side of the equation. 262 00:18:38,130 --> 00:18:40,700 That makes the stability better, but, of course, 263 00:18:40,700 --> 00:18:46,100 it's going to make the computation more difficult. 264 00:18:46,100 --> 00:18:48,500 Let me just think about the stability for now. 265 00:18:48,500 --> 00:18:51,950 So, instead of this board, this is the explicit case. 266 00:18:57,136 --> 00:18:57,781 Yeah. 267 00:18:57,781 --> 00:18:58,656 AUDIENCE: [INAUDIBLE] 268 00:19:05,930 --> 00:19:08,820 PROFESSOR: So it has to be, let's draw a little picture 269 00:19:08,820 --> 00:19:10,440 here. 270 00:19:10,440 --> 00:19:11,890 So, there's 0. 271 00:19:11,890 --> 00:19:12,840 There's 1. 272 00:19:12,840 --> 00:19:15,760 And here's minus 1. 273 00:19:15,760 --> 00:19:21,640 And this quantity is starting at 1, and subtracting 4R. 274 00:19:21,640 --> 00:19:24,280 So, the question is, does it stop here? 275 00:19:24,280 --> 00:19:30,810 If it stops there, then in this picture, 276 00:19:30,810 --> 00:19:32,060 these would both be negative. 277 00:19:32,060 --> 00:19:34,240 This would be about minus a half, 278 00:19:34,240 --> 00:19:36,800 but that is greater than minus 1. 279 00:19:36,800 --> 00:19:41,000 So this would be the OK case, when this is true. 280 00:19:41,000 --> 00:19:44,820 And that's what required R less or equal to half. 281 00:19:44,820 --> 00:19:50,470 So in this case, yeah. 282 00:19:50,470 --> 00:19:56,780 In that case, R is less than -- but if R equal to half -- 283 00:19:56,780 --> 00:19:59,340 everybody sees this -- if R equaled a half, 284 00:19:59,340 --> 00:20:01,060 this would be 1 minus 2. 285 00:20:01,060 --> 00:20:03,320 This would bring us all the way here. 286 00:20:03,320 --> 00:20:07,260 And if R is bigger than a half, then that 287 00:20:07,260 --> 00:20:12,450 would reverse that inequality, and we would be unstable. 288 00:20:12,450 --> 00:20:12,950 Yeah. 289 00:20:12,950 --> 00:20:15,130 I think it's OK. 290 00:20:15,130 --> 00:20:15,630 Thanks. 291 00:20:19,040 --> 00:20:22,760 It's the difference between looking at the magnitude -- 292 00:20:22,760 --> 00:20:27,830 if I put magnitude, then I would need less or equal or just 293 00:20:27,830 --> 00:20:30,840 realizing that it's just taking it as a number. 294 00:20:30,840 --> 00:20:32,501 1 minus 4. 295 00:20:32,501 --> 00:20:33,000 OK. 296 00:20:33,000 --> 00:20:34,600 I think that's all right. 297 00:20:34,600 --> 00:20:37,440 Now the implicit case. 298 00:20:37,440 --> 00:20:39,870 What happened to G? 299 00:20:39,870 --> 00:20:42,490 What happened to G? 300 00:20:42,490 --> 00:20:43,820 Let me put it here. 301 00:20:43,820 --> 00:20:44,400 Implicit. 302 00:20:50,560 --> 00:20:53,240 So, what's changed? 303 00:20:53,240 --> 00:20:57,630 This term, which came from the second difference, 304 00:20:57,630 --> 00:21:02,350 is moved over to the n plus 1 side of the equation. 305 00:21:02,350 --> 00:21:12,830 So, we have -- in Fourier space, we would have 1 from the u_(j, 306 00:21:12,830 --> 00:21:21,240 n+1), and then minus 2R, and times our expression 2 cosine k 307 00:21:21,240 --> 00:21:31,730 delta x minus 2 -- that's what multiplies e to the i -- 308 00:21:31,730 --> 00:21:39,900 do you see that that's showing up on the new side 309 00:21:39,900 --> 00:21:42,900 of the equation, multiplying e to the i*k delta x, 310 00:21:42,900 --> 00:21:57,280 and on the old side, the old time step, we only have the 1. 311 00:21:57,280 --> 00:21:58,960 In other words, the growth factor 312 00:21:58,960 --> 00:22:02,150 is, this goes down into the denominator, 313 00:22:02,150 --> 00:22:04,810 and we can see what's going on. 314 00:22:04,810 --> 00:22:09,520 So the growth factor is 1 over this expression. 315 00:22:09,520 --> 00:22:18,910 1 minus 2R times 2 cosine k delta x minus 2. 316 00:22:18,910 --> 00:22:20,780 We keep running into that quantity 317 00:22:20,780 --> 00:22:22,670 from the second difference. 318 00:22:22,670 --> 00:22:23,310 OK. 319 00:22:23,310 --> 00:22:27,180 So if I look at that, what's going on? 320 00:22:27,180 --> 00:22:29,005 The key point of this quantity is 321 00:22:29,005 --> 00:22:35,390 that it's always negative, or zero, but never positive. 322 00:22:35,390 --> 00:22:42,150 Because this term is never as large as this one. 323 00:22:42,150 --> 00:22:43,240 So it's negative. 324 00:22:43,240 --> 00:22:45,160 And now, I have another minus sign here, 325 00:22:45,160 --> 00:22:48,380 so I have 1 divided by 1 plus something. 326 00:22:48,380 --> 00:22:56,830 And always, for all k and all ratios R, it's OK. 327 00:22:56,830 --> 00:22:58,060 Right? 328 00:22:58,060 --> 00:23:00,460 1 over 1 plus something positive. 329 00:23:00,460 --> 00:23:02,390 And you see that everything's real. 330 00:23:02,390 --> 00:23:05,610 That's because we're doing the diffusion equation, where 331 00:23:05,610 --> 00:23:09,970 the wave equation took us, we had some imaginary part. 332 00:23:15,130 --> 00:23:17,050 Absolutely stable. 333 00:23:17,050 --> 00:23:18,490 Stable for all delta t. 334 00:23:24,320 --> 00:23:26,740 I'll say even for large delta t. 335 00:23:31,120 --> 00:23:32,090 That's good. 336 00:23:32,090 --> 00:23:40,970 Of course, we're not going to take delta t equal to 1000. 337 00:23:40,970 --> 00:23:45,640 Even if stability allowed us, the accuracy would be hopeless. 338 00:23:45,640 --> 00:23:47,120 It'd be terrible. 339 00:23:47,120 --> 00:23:54,430 I can't take a giant time step and expect to follow the true 340 00:23:54,430 --> 00:23:58,330 exponentials for e to the i*k*x. 341 00:23:58,330 --> 00:24:05,870 But, I can take a larger step than this very small delta, 342 00:24:05,870 --> 00:24:08,120 one half delta x squared step. 343 00:24:08,120 --> 00:24:10,060 So this good. 344 00:24:10,060 --> 00:24:11,430 What's the price that we pay? 345 00:24:14,320 --> 00:24:19,660 The price is that at every time step, we have to solve. 346 00:24:19,660 --> 00:24:25,030 This is now on the other side, so this brings a whole matrix 347 00:24:25,030 --> 00:24:26,720 over. 348 00:24:26,720 --> 00:24:33,670 So, this implicit case -- when I bring this second difference 349 00:24:33,670 --> 00:24:41,930 matrix over -- I've used capital K, 350 00:24:41,930 --> 00:24:46,200 because the second difference matrix there is probably 351 00:24:46,200 --> 00:24:52,500 the most interesting difference matrix, I'll say. 352 00:24:52,500 --> 00:24:56,350 One of the most interesting matrices, period. 353 00:24:56,350 --> 00:25:00,320 It has the minus 2's on the diagonal and 1's 354 00:25:00,320 --> 00:25:03,430 above the diagonal. 355 00:25:03,430 --> 00:25:06,950 So, actually, with signs exchanged, 356 00:25:06,950 --> 00:25:10,570 I called it K. So can I just remember 357 00:25:10,570 --> 00:25:15,210 the meaning of this matrix K. I'm 358 00:25:15,210 --> 00:25:19,750 going to stay with that letter because it's 359 00:25:19,750 --> 00:25:24,220 standard in finite elements for the stiffness matrix. 360 00:25:24,220 --> 00:25:27,480 So, K, by definition, is the matrix with the 2's 361 00:25:27,480 --> 00:25:32,940 on the diagonal, minus 1's above the diagonal, and minus 1's 362 00:25:32,940 --> 00:25:35,170 below the diagonal. 363 00:25:35,170 --> 00:25:37,340 And, what size is it? 364 00:25:40,100 --> 00:25:46,130 On the whole line, where we're working for simplicity, 365 00:25:46,130 --> 00:25:51,850 I guess K is infinite in both directions. 366 00:25:51,850 --> 00:25:54,000 2's and minus 1's. 367 00:25:54,000 --> 00:25:58,910 On an interval where we would really compute, 368 00:25:58,910 --> 00:26:01,240 K would be a finite matrix. 369 00:26:01,240 --> 00:26:04,395 But, the main point is that this matrix K is minus 370 00:26:04,395 --> 00:26:06,480 the second difference matrix. 371 00:26:09,460 --> 00:26:14,900 So that when I bring it over to this side, which 372 00:26:14,900 --> 00:26:20,750 changes that plus to a minus, I get this equation to solve. 373 00:26:20,750 --> 00:26:31,100 The identity, coming from the 1, plus R, the ratio, 374 00:26:31,100 --> 00:26:34,800 times this matrix K, multiplies this vector 375 00:26:34,800 --> 00:26:38,900 of all the values at the new time step, 376 00:26:38,900 --> 00:26:45,590 and equals, in this case, all the values at the old time 377 00:26:45,590 --> 00:26:46,850 step. 378 00:26:46,850 --> 00:26:51,370 Because, at the old time step, I just have the identity. 379 00:26:51,370 --> 00:26:52,860 Are you OK with this? 380 00:26:52,860 --> 00:26:58,010 This is the problem that we actually have to solve. 381 00:26:58,010 --> 00:27:00,430 I mean, the main point about it is 382 00:27:00,430 --> 00:27:03,910 that we have a linear system. 383 00:27:03,910 --> 00:27:07,440 The u's at the new time step are coupled. 384 00:27:07,440 --> 00:27:09,430 Coupled by the second difference. 385 00:27:09,430 --> 00:27:12,830 Coupled by the matrix K, so we have to deal with it. 386 00:27:12,830 --> 00:27:17,200 Well actually, in one dimension, no problem. 387 00:27:17,200 --> 00:27:22,980 Because in one dimension, K is just a tri-diagonal matrix, 388 00:27:22,980 --> 00:27:29,360 and we can solve tri-diagonal systems 389 00:27:29,360 --> 00:27:31,930 almost as fast as we could work with the identity matrix. 390 00:27:34,870 --> 00:27:40,820 And actually, what we're seeing in this ratio 391 00:27:40,820 --> 00:27:46,850 are exactly the eigenvalues of the inverse matrix. 392 00:27:46,850 --> 00:27:51,910 So, if I use a formula, I would say U_(n+1) is the inverse 393 00:27:51,910 --> 00:27:58,910 matrix times U_n, and the growth factor is just the eigenvalues. 394 00:27:58,910 --> 00:28:02,490 The growth factor G is simply the eigenvalues 395 00:28:02,490 --> 00:28:05,070 of this inverse matrix. 396 00:28:05,070 --> 00:28:12,800 And, notice again, that all the eigenvalues of K 397 00:28:12,800 --> 00:28:15,350 are greater or equal zero. 398 00:28:15,350 --> 00:28:17,840 It was a positive definite matrix, 399 00:28:17,840 --> 00:28:20,210 because I switched sign. 400 00:28:20,210 --> 00:28:22,940 And, the effect of the identity is 401 00:28:22,940 --> 00:28:25,760 to push it a little more positive. 402 00:28:25,760 --> 00:28:26,730 Here. 403 00:28:26,730 --> 00:28:31,550 It's making -- this is a positive number there. 404 00:28:31,550 --> 00:28:34,460 That's the eigenvalues of R*K, and then, 405 00:28:34,460 --> 00:28:37,100 I making it a little more positive with this 1, 406 00:28:37,100 --> 00:28:49,910 so we have a problem that is safely invertible, stable, 407 00:28:49,910 --> 00:28:52,330 and in one dimension, quite simple. 408 00:28:52,330 --> 00:29:00,900 So in one dimension we're entirely ready to go implicit. 409 00:29:00,900 --> 00:29:04,570 Let me just, while I'm thinking about it, 410 00:29:04,570 --> 00:29:09,630 recall that we had some other matrix in the finite matrix 411 00:29:09,630 --> 00:29:19,270 case, we had some other matrices in the very first section 412 00:29:19,270 --> 00:29:25,350 of the original notes, which had different boundary conditions. 413 00:29:25,350 --> 00:29:32,110 In other words, if we're on a finite interval, 414 00:29:32,110 --> 00:29:37,560 then this matrix K corresponds to zero boundary conditions. 415 00:29:37,560 --> 00:29:39,790 Zero temperature at the end. 416 00:29:39,790 --> 00:29:43,970 But we might have other boundary conditions, 417 00:29:43,970 --> 00:29:46,310 like the derivative of the temperature 418 00:29:46,310 --> 00:29:49,220 is zero at one end, or the other end. 419 00:29:49,220 --> 00:29:54,450 That would change this matrix in the first row and the last row 420 00:29:54,450 --> 00:29:55,880 to another nice matrix. 421 00:29:55,880 --> 00:29:59,520 So, we might have not K but one of the other tri-diagonal 422 00:29:59,520 --> 00:30:02,630 matrices. 423 00:30:02,630 --> 00:30:06,660 But, here's the real point. 424 00:30:06,660 --> 00:30:08,780 What happens in two dimensions? 425 00:30:08,780 --> 00:30:18,250 Suppose I include now a u_yy term in the heat equation, 426 00:30:18,250 --> 00:30:21,470 or three dimensions with a u_zz term. 427 00:30:21,470 --> 00:30:25,280 What's different? 428 00:30:25,280 --> 00:30:29,260 Well, we can follow through the stability 429 00:30:29,260 --> 00:30:33,840 test for the explicit method, and it'll all 430 00:30:33,840 --> 00:30:37,640 have a term from delta x squared, 431 00:30:37,640 --> 00:30:40,290 and a term from delta y squared, and a delta z squared. 432 00:30:40,290 --> 00:30:43,940 It'll be even slightly tighter on delta t. 433 00:30:43,940 --> 00:30:48,740 Or I can follow through the implicit method, where 434 00:30:48,740 --> 00:30:51,900 now this second difference in the x direction 435 00:30:51,900 --> 00:30:54,940 appears also with a second difference in the y direction 436 00:30:54,940 --> 00:30:58,370 and a second different in the z direction. 437 00:30:58,370 --> 00:31:02,050 So what's happening? 438 00:31:02,050 --> 00:31:03,990 From the matrix theory point of view, 439 00:31:03,990 --> 00:31:08,920 I still have nice matrices with 2's and minus 1's or maybe with 440 00:31:08,920 --> 00:31:12,830 -- I guess what we would have -- can we just remember what we 441 00:31:12,830 --> 00:31:16,950 would have in the two-dimensional case. 442 00:31:16,950 --> 00:31:19,150 In the two-dimensional case, instead of 443 00:31:19,150 --> 00:31:24,700 this second difference with a minus 1, a 2, and a minus 1 444 00:31:24,700 --> 00:31:29,150 in the x direction, we also have a second difference 445 00:31:29,150 --> 00:31:34,260 coming from the u_yy in the y direction, so another minus 1, 446 00:31:34,260 --> 00:31:38,480 and that 2 moves up to a 4. 447 00:31:38,480 --> 00:31:43,840 So our matrix in the 2D, two space dimension case, 448 00:31:43,840 --> 00:31:47,520 will have a 4's on the diagonal. 449 00:31:47,520 --> 00:31:50,000 So it grows from this one. 450 00:31:50,000 --> 00:31:55,230 It's called the tensor product coming out of this. 451 00:31:55,230 --> 00:31:59,910 Anyway, it produces this famous five-point molecule 452 00:31:59,910 --> 00:32:04,220 with 4's on the diagonal, and four minus 1's 453 00:32:04,220 --> 00:32:06,960 going down other diagonals. 454 00:32:06,960 --> 00:32:12,930 And the problem is that those diagonal are not -- 455 00:32:12,930 --> 00:32:15,430 two of the diagonals might be right next to the main 456 00:32:15,430 --> 00:32:19,570 diagonal, but the other two will be farther away. 457 00:32:19,570 --> 00:32:22,240 Because we can't number the nodes to keep all 458 00:32:22,240 --> 00:32:28,000 these four numbered adjacently at all the mesh points. 459 00:32:28,000 --> 00:32:30,510 In other words, we have a serious matrix problem. 460 00:32:30,510 --> 00:32:32,490 That's my comment. 461 00:32:32,490 --> 00:32:38,670 We have a matrix K in 2D. 462 00:32:38,670 --> 00:32:41,820 So what I'm speaking about is this problem. 463 00:32:41,820 --> 00:32:49,980 Identity plus R*K(2D) u at the new time equal u at the old 464 00:32:49,980 --> 00:32:50,810 time. 465 00:32:50,810 --> 00:32:54,370 And K(2D) is by no means so easy to deal with as K(1D). 466 00:32:58,910 --> 00:33:02,770 It doesn't have a nice narrow band like this. 467 00:33:02,770 --> 00:33:05,010 The band stretches further, and we'll 468 00:33:05,010 --> 00:33:11,120 see in spades, in then next part of the course which 469 00:33:11,120 --> 00:33:17,410 is about solving large linear systems. 470 00:33:17,410 --> 00:33:19,150 So we have a large system because we've 471 00:33:19,150 --> 00:33:21,560 got lots of mesh points. 472 00:33:21,560 --> 00:33:25,840 We have matrix that has a wider band, 473 00:33:25,840 --> 00:33:28,850 and so you need new ideas in solving it. 474 00:33:28,850 --> 00:33:36,160 So this problem is what will be a model -- is the connection, 475 00:33:36,160 --> 00:33:40,000 really, between this section of the course, 476 00:33:40,000 --> 00:33:42,380 difference methods for initial-value problems, 477 00:33:42,380 --> 00:33:47,050 and the next section, which is solving large linear systems. 478 00:33:47,050 --> 00:33:53,570 Because this will be a typical model large linear system. 479 00:33:53,570 --> 00:33:59,500 So that's a comment that for implicit equations, 480 00:33:59,500 --> 00:34:04,130 the price was low in 1D, but the price is not so low in 2D. 481 00:34:04,130 --> 00:34:05,040 OK. 482 00:34:05,040 --> 00:34:10,210 Now, have we got any other candidates 483 00:34:10,210 --> 00:34:14,440 between explicit and implicit? 484 00:34:14,440 --> 00:34:20,810 Well, we did, way back, for ordinary differential 485 00:34:20,810 --> 00:34:22,510 equations, we had a trapezoidal rule. 486 00:34:22,510 --> 00:34:25,170 And what did that do? 487 00:34:25,170 --> 00:34:29,530 That increased the accuracy by centering it 488 00:34:29,530 --> 00:34:31,190 at time n plus a half. 489 00:34:31,190 --> 00:34:32,810 Do you remember that? 490 00:34:32,810 --> 00:34:35,160 It was half on the right-hand side; 491 00:34:35,160 --> 00:34:38,830 we used half the explicit difference 492 00:34:38,830 --> 00:34:40,590 and half of the implicit difference, 493 00:34:40,590 --> 00:34:44,020 so natural to do it again. 494 00:34:44,020 --> 00:34:47,140 It worked then, and it'll work now. 495 00:34:47,140 --> 00:34:49,570 So this will be the trapezoidal method. 496 00:34:49,570 --> 00:34:57,310 But most people in the diffusion application 497 00:34:57,310 --> 00:35:03,320 name it after these authors, Crank and Nicolson. 498 00:35:03,320 --> 00:35:06,560 So it's the Crank-Nicolson method. 499 00:35:06,560 --> 00:35:12,990 It's half the explicit part plus half the implicit part. 500 00:35:12,990 --> 00:35:17,390 And the order of accuracy is increased 501 00:35:17,390 --> 00:35:19,420 by centering it this way. 502 00:35:19,420 --> 00:35:22,520 That's the point, of course. 503 00:35:22,520 --> 00:35:25,320 So I guess the point will be we don't 504 00:35:25,320 --> 00:35:33,560 have any extra work compared to a fully implicit one, 505 00:35:33,560 --> 00:35:37,120 and we get an extra order of accuracy, 506 00:35:37,120 --> 00:35:41,350 by centering it at a half. 507 00:35:41,350 --> 00:35:55,020 So u_(j, n+1) minus u_(j, n) over delta t is R, 508 00:35:55,020 --> 00:35:59,320 that ratio -- oh sorry. 509 00:35:59,320 --> 00:36:07,220 We've got delta t there -- so will be the second difference 510 00:36:07,220 --> 00:36:24,210 at the old time times a half -- it would just take half of that 511 00:36:24,210 --> 00:36:29,040 explicit part plus half of the implicit part -- 512 00:36:29,040 --> 00:36:34,830 this is at time n plus one, over delta x squared. 513 00:36:34,830 --> 00:36:36,050 You would think of that. 514 00:36:36,050 --> 00:36:37,730 Everybody would. 515 00:36:37,730 --> 00:36:39,830 It's just a natural idea. 516 00:36:39,830 --> 00:36:43,820 It increases the accuracy. 517 00:36:43,820 --> 00:36:45,100 It's still implicit. 518 00:36:45,100 --> 00:36:50,050 We still have to bring this term over onto this side when u_(j, 519 00:36:50,050 --> 00:36:53,850 n) moves over onto that side. 520 00:36:53,850 --> 00:36:56,470 And we can look at the growth factor. 521 00:36:56,470 --> 00:36:58,870 So that's what we have to do. 522 00:36:58,870 --> 00:37:04,140 So bring that term over, and bring delta t up, 523 00:37:04,140 --> 00:37:08,750 so it's delta t divided by delta x squared and R. OK. 524 00:37:08,750 --> 00:37:10,790 So I've got this term. 525 00:37:10,790 --> 00:37:12,060 The 1. 526 00:37:12,060 --> 00:37:13,840 And when I bring that term over, it's 527 00:37:13,840 --> 00:37:17,340 a half of what we brought over before, 528 00:37:17,340 --> 00:37:22,790 which was the R, 2 cosine k delta x minus the 2. 529 00:37:25,400 --> 00:37:29,750 That multiplies G. And on the right side, 530 00:37:29,750 --> 00:37:35,450 we have the 1 from the u_(j, n) plus the half, 531 00:37:35,450 --> 00:37:41,080 times the R times 2 cosine k delta x minus 2. 532 00:37:47,450 --> 00:37:50,720 I'm certainly getting to the point where I write these 533 00:37:50,720 --> 00:37:57,260 expressions down faster without explicitly writing e to the i*k 534 00:37:57,260 --> 00:38:00,650 delta x and then canceling it. 535 00:38:00,650 --> 00:38:05,650 So, G is a ratio again. 536 00:38:05,650 --> 00:38:14,480 In the explicit case G was all in the numerator. 537 00:38:14,480 --> 00:38:16,490 There was no denominator. 538 00:38:16,490 --> 00:38:24,480 In the fully implicit case, g was all denominator, 539 00:38:24,480 --> 00:38:27,460 with just a one in the numerator. 540 00:38:27,460 --> 00:38:30,250 Now we have G as the ratio. 541 00:38:30,250 --> 00:38:32,660 So I'm going to bring this down below. 542 00:38:39,000 --> 00:38:40,850 And, G is a ratio. 543 00:38:40,850 --> 00:38:43,360 So it's a fraction because it's got 544 00:38:43,360 --> 00:38:48,810 an explicit part in the numerator and an implicit part 545 00:38:48,810 --> 00:38:49,690 in the denominator. 546 00:38:49,690 --> 00:38:51,400 And, of course, the main question 547 00:38:51,400 --> 00:38:54,090 is, when is G less or equal to one? 548 00:38:59,570 --> 00:39:00,640 Well. 549 00:39:00,640 --> 00:39:07,260 Actually it's cool, because G is simply -- look, 550 00:39:07,260 --> 00:39:11,250 the numerator is one, and all this stuff is negative. 551 00:39:11,250 --> 00:39:15,070 So this is one minus something. 552 00:39:15,070 --> 00:39:18,630 And, what's going down here? 553 00:39:18,630 --> 00:39:20,020 The signs just reverse. 554 00:39:20,020 --> 00:39:22,270 It's one plus something. 555 00:39:22,270 --> 00:39:25,710 And this is a positive quantity, and this 556 00:39:25,710 --> 00:39:28,450 is the same positive quantity. 557 00:39:28,450 --> 00:39:31,590 So, this is one minus -- do you see this -- 558 00:39:31,590 --> 00:39:34,955 because this is negative, I'll write it as 1 minus something 559 00:39:34,955 --> 00:39:36,660 -- something positive. 560 00:39:36,660 --> 00:39:40,700 And here, this is minus and this is negative, 561 00:39:40,700 --> 00:39:43,250 so this is 1 plus the same thing. 562 00:39:43,250 --> 00:39:44,180 So what's the answer? 563 00:39:51,830 --> 00:39:53,710 [UNINTELLIGIBLE PHRASE] 564 00:39:53,710 --> 00:39:57,390 One minus a positive number over 1 plus that positive number 565 00:39:57,390 --> 00:39:58,390 can't go wrong. 566 00:40:02,410 --> 00:40:09,230 So, for all R. OK. 567 00:40:09,230 --> 00:40:17,440 In a way, that has to be more attractive than fully implicit, 568 00:40:17,440 --> 00:40:19,340 because it didn't take more work, 569 00:40:19,340 --> 00:40:23,400 it wasn't any less stable, and it was more accurate. 570 00:40:23,400 --> 00:40:28,740 So Crank-Nicolson is sort of a natural idea 571 00:40:28,740 --> 00:40:33,420 at this level of analysis. 572 00:40:36,080 --> 00:40:36,770 OK. 573 00:40:36,770 --> 00:40:39,970 So those are the three methods to think about. 574 00:40:39,970 --> 00:40:49,720 Of course, actual computations will bring in new questions. 575 00:40:49,720 --> 00:40:51,906 Boundary conditions, non-linearities, 576 00:40:51,906 --> 00:40:52,530 and everything. 577 00:40:52,530 --> 00:40:55,280 I'm just staying at this simple level, 578 00:40:55,280 --> 00:41:00,850 where I'm comparing the first ideas that would occur to us. 579 00:41:00,850 --> 00:41:01,680 OK. 580 00:41:01,680 --> 00:41:05,250 We may want higher accuracy, all sorts of things. 581 00:41:05,250 --> 00:41:12,780 But, we're seeing the contrast with the wave equation case. 582 00:41:12,780 --> 00:41:19,060 But, here we move more toward the implicit side, 583 00:41:19,060 --> 00:41:25,930 and the reward is unconditional stability, no condition on r, 584 00:41:25,930 --> 00:41:29,100 but the price is more computation. 585 00:41:29,100 --> 00:41:31,080 OK. 586 00:41:31,080 --> 00:41:35,770 Well, at that same level of thinking about the most 587 00:41:35,770 --> 00:41:41,040 natural methods that occur to you, 588 00:41:41,040 --> 00:41:45,860 I want to speak now about the big, big issue 589 00:41:45,860 --> 00:41:47,430 of convection-diffusion. 590 00:41:47,430 --> 00:41:52,370 Suppose both types of problem are with us. 591 00:41:52,370 --> 00:41:55,320 So convention-diffusion. 592 00:41:55,320 --> 00:41:56,340 What would you do then? 593 00:42:01,520 --> 00:42:11,110 Now, I have both convection and diffusion in the equation u_t 594 00:42:11,110 --> 00:42:14,650 equals c*u_x plus d*u_xx. 595 00:42:18,650 --> 00:42:22,920 And I guess I want to say that this is something 596 00:42:22,920 --> 00:42:25,700 that the world has not come to a decision on. 597 00:42:25,700 --> 00:42:29,720 The computational world is still debating 598 00:42:29,720 --> 00:42:32,290 what to do with convection-diffusion. 599 00:42:32,290 --> 00:42:43,740 So it's a very appropriate equation model to experiment 600 00:42:43,740 --> 00:42:44,240 with. 601 00:42:47,240 --> 00:42:53,500 Maybe I'll ask you, as the next, kind of informal homework, 602 00:42:53,500 --> 00:42:58,660 to begin to experiment with that. 603 00:42:58,660 --> 00:43:02,050 And when I say experiment, you can think immediately 604 00:43:02,050 --> 00:43:07,700 of the methods you might think of. 605 00:43:07,700 --> 00:43:09,850 The first methods to try. 606 00:43:09,850 --> 00:43:12,290 So the first method of course, will be explicit. 607 00:43:16,580 --> 00:43:19,850 Fully explicit. 608 00:43:19,850 --> 00:43:22,110 What do you expect to have happen there? 609 00:43:22,110 --> 00:43:25,690 Let's just see how the two different terms 610 00:43:25,690 --> 00:43:33,100 each give their little push towards the stability 611 00:43:33,100 --> 00:43:34,460 requirement. 612 00:43:34,460 --> 00:43:44,250 So explicit is going to be the time difference, forward 613 00:43:44,250 --> 00:43:49,970 time difference equals c times-- what shall I choose first 614 00:43:49,970 --> 00:43:52,690 for the space difference? 615 00:43:55,230 --> 00:43:57,030 Upwind, maybe? 616 00:43:57,030 --> 00:43:59,330 Upwind in space? 617 00:43:59,330 --> 00:44:04,540 So, this is going to be explicit and upwind. 618 00:44:04,540 --> 00:44:09,330 Of course, everybody knows that refers to the windy term 619 00:44:09,330 --> 00:44:13,441 in the equation, the wave term. 620 00:44:13,441 --> 00:44:13,940 OK. 621 00:44:13,940 --> 00:44:17,510 So, of course, I'll get a forward difference, then, 622 00:44:17,510 --> 00:44:23,150 in the x direction of u at time n divided delta x, 623 00:44:23,150 --> 00:44:27,080 and then I'm going to do the natural one. 624 00:44:27,080 --> 00:44:35,120 A centered second difference at time n for the diffusion term. 625 00:44:35,120 --> 00:44:35,620 OK. 626 00:44:39,660 --> 00:44:43,040 We'll see what comes out of that. 627 00:44:43,040 --> 00:44:45,510 It's going to involve not only the ratio, 628 00:44:45,510 --> 00:44:48,900 it will involve both ratios, naturally. 629 00:44:48,900 --> 00:44:52,890 Little r, as soon as I multiply up by delta x, 630 00:44:52,890 --> 00:44:57,830 I have a little r there, my wave equation ratio. 631 00:44:57,830 --> 00:45:00,520 And I have a d delta t over delta x squared, 632 00:45:00,520 --> 00:45:06,710 which I'm calling capital R. That's the diffusion ratio. 633 00:45:06,710 --> 00:45:11,510 And multiplied by all sorts of stuff. 634 00:45:11,510 --> 00:45:15,330 So actually, why don't I write down what the equation really 635 00:45:15,330 --> 00:45:16,000 is then? 636 00:45:16,000 --> 00:45:21,830 Explicitly, it's u_(j, n+1) is u_(j, n) -- 637 00:45:21,830 --> 00:45:27,220 so maybe I should collect all the terms here. 638 00:45:27,220 --> 00:45:33,360 Can I collect the terms that multiply u_(j+1, n)? 639 00:45:38,080 --> 00:45:40,780 So what's going to multiply u_(j+1, n)? 640 00:45:40,780 --> 00:45:43,870 Here I have a c delta t over delta x, 641 00:45:43,870 --> 00:45:46,690 I have an r, from this term. 642 00:45:46,690 --> 00:45:49,860 And here I have a d. 643 00:45:49,860 --> 00:45:50,360 Right? 644 00:45:53,720 --> 00:45:54,760 No. 645 00:45:54,760 --> 00:45:55,260 Sorry. 646 00:45:55,260 --> 00:46:00,100 I have a capital R. Isn't that right? 647 00:46:00,100 --> 00:46:03,620 When I multiply through by delta t, I have d delta t over delta 648 00:46:03,620 --> 00:46:04,840 x squared. 649 00:46:04,840 --> 00:46:08,660 I have a capital R that's going to multiply the -- 650 00:46:08,660 --> 00:46:11,690 does that look right? 651 00:46:14,580 --> 00:46:18,900 And then what multiplies the center guy, the u_(j, n)? 652 00:46:23,980 --> 00:46:29,090 Well that's where the time difference comes over. 653 00:46:29,090 --> 00:46:31,480 There's the usual 1. 654 00:46:31,480 --> 00:46:36,760 Then from this difference, this upwind difference, 655 00:46:36,760 --> 00:46:38,810 there's a minus at the center term, 656 00:46:38,810 --> 00:46:41,110 because it's a forward difference. 657 00:46:41,110 --> 00:46:48,160 So that will be a minus r, and from this centered difference, 658 00:46:48,160 --> 00:46:51,270 well remember, that that's a second difference, 659 00:46:51,270 --> 00:46:57,660 so there's a coefficient 2 in the middle term, or minus 2, 660 00:46:57,660 --> 00:47:05,170 and when I bring the delta t up, it's minus 2 capital R. 661 00:47:05,170 --> 00:47:10,150 And then, the third one, whatever this last coefficient 662 00:47:10,150 --> 00:47:16,220 is, is multiplying u_(j, n-1), and now, where does that come? 663 00:47:16,220 --> 00:47:18,010 It doesn't come from here, it doesn't 664 00:47:18,010 --> 00:47:19,530 come from the upwind difference, I 665 00:47:19,530 --> 00:47:22,490 guess it only comes from there. 666 00:47:22,490 --> 00:47:30,220 And it's probably just R. And I guess 667 00:47:30,220 --> 00:47:35,730 I would check any formula like that, 668 00:47:35,730 --> 00:47:40,000 to be sure that the coefficients add to 1. 669 00:47:43,670 --> 00:47:45,700 And they do. 670 00:47:45,700 --> 00:47:46,990 And why should they? 671 00:47:46,990 --> 00:47:50,910 Why should those coefficients add to 1? 672 00:47:50,910 --> 00:47:57,940 So that's my little test that I've got it right. 673 00:47:57,940 --> 00:48:04,945 They should add to 1, because a constant initial -- if u_(j, 674 00:48:04,945 --> 00:48:09,490 0) is a constant, if I'm starting from a constant 675 00:48:09,490 --> 00:48:15,330 function, its x derivatives would all be zero. 676 00:48:15,330 --> 00:48:18,757 If I had u equal constant, then the x derivatives 677 00:48:18,757 --> 00:48:20,340 would be zero, so the time derivatives 678 00:48:20,340 --> 00:48:24,380 would be zero, so that constant would stay forever. 679 00:48:24,380 --> 00:48:27,790 And so I would want to get the same constant out 680 00:48:27,790 --> 00:48:31,040 from constants going in. 681 00:48:31,040 --> 00:48:35,480 And so the coefficient should add to one. 682 00:48:35,480 --> 00:48:36,190 Well. 683 00:48:40,830 --> 00:48:45,130 I can see a situation when I know stability will be OK. 684 00:48:48,420 --> 00:48:53,520 Stability can't fail if those three coefficients are all 685 00:48:53,520 --> 00:48:55,380 positive. 686 00:48:55,380 --> 00:48:59,610 If those three coefficients are all positive, 687 00:48:59,610 --> 00:49:04,760 then my growth factor can't be larger than 1. 688 00:49:04,760 --> 00:49:07,030 Those coefficients, they add to 1. 689 00:49:07,030 --> 00:49:11,600 So the growth factor is this times -- 690 00:49:11,600 --> 00:49:12,990 should I write it down? 691 00:49:12,990 --> 00:49:15,950 The growth factor will be this -- 692 00:49:15,950 --> 00:49:17,450 so this will be the growth factor -- 693 00:49:17,450 --> 00:49:22,760 G will be that number times e to the i*k delta x, 694 00:49:22,760 --> 00:49:25,530 from shifting over one. 695 00:49:25,530 --> 00:49:31,780 It'll be this number times 1, and it'll be this number times 696 00:49:31,780 --> 00:49:37,970 e to the minus i*k delta x, and my point is just that 697 00:49:37,970 --> 00:49:42,210 the easiest estimate would say, well this has magnitude 1, 698 00:49:42,210 --> 00:49:45,840 that has magnitude 1, that has magnitude 1. 699 00:49:45,840 --> 00:49:49,560 If I just add up magnitudes, I'll 700 00:49:49,560 --> 00:49:53,570 get this one, and this one and this one. 701 00:49:53,570 --> 00:49:57,520 And if they're all positive, they add to 1, I'm safe. 702 00:49:57,520 --> 00:50:00,130 So what's the stability condition 703 00:50:00,130 --> 00:50:04,420 that this simple argument requires? 704 00:50:04,420 --> 00:50:05,770 Well. 705 00:50:05,770 --> 00:50:06,900 That's positive. 706 00:50:06,900 --> 00:50:07,930 That's positive. 707 00:50:07,930 --> 00:50:09,440 The danger is this one. 708 00:50:13,040 --> 00:50:16,340 So I guess I would like that to be positive. 709 00:50:16,340 --> 00:50:23,250 So I would want -- I'm certainly stable if this is positive, 710 00:50:23,250 --> 00:50:30,971 which I might as well write as r plus 2 big R being not larger 711 00:50:30,971 --> 00:50:31,470 than 1. 712 00:50:37,590 --> 00:50:41,050 Because if that's true, I'm taking a combination 713 00:50:41,050 --> 00:50:45,030 with positive coefficients adding to 1, I can't go wrong. 714 00:50:45,030 --> 00:50:45,530 OK. 715 00:50:45,530 --> 00:50:51,290 And in fact, that's pretty familiar. 716 00:50:55,450 --> 00:50:58,950 It combines the two conditions we know. 717 00:50:58,950 --> 00:51:01,200 The condition that r should be less or equal 718 00:51:01,200 --> 00:51:03,840 1, small r, the wave Courant-Friedrichs-Lewy 719 00:51:03,840 --> 00:51:06,300 condition. 720 00:51:06,300 --> 00:51:10,260 And the condition that capital R should 721 00:51:10,260 --> 00:51:14,010 be less or equal to half the diffusion condition. 722 00:51:14,010 --> 00:51:19,440 Now they're both -- now I'm forcing them even more 723 00:51:19,440 --> 00:51:23,251 by saying that the sum has to be less or equal 1. 724 00:51:23,251 --> 00:51:23,750 OK. 725 00:51:26,720 --> 00:51:35,600 So let me leave the explicit upwind at this point, 726 00:51:35,600 --> 00:51:40,270 and ask you what other method might you want to try? 727 00:51:40,270 --> 00:51:43,100 Well, it would be natural to think 728 00:51:43,100 --> 00:51:51,330 of explicit centered, where you would center this 729 00:51:51,330 --> 00:51:54,120 and we could follow that one through. 730 00:51:56,680 --> 00:51:58,750 I guess I'll do that next time very quickly, 731 00:51:58,750 --> 00:52:01,340 and the notes will do it. 732 00:52:01,340 --> 00:52:05,650 And I'm not sure which of those two ways is better. 733 00:52:05,650 --> 00:52:11,080 Or, you're going to say go implicit. 734 00:52:11,080 --> 00:52:14,520 And now I'm going to ask, all right, but what part of it 735 00:52:14,520 --> 00:52:15,890 do you want me to make implicit? 736 00:52:19,720 --> 00:52:22,000 Maybe all I need to make implicit 737 00:52:22,000 --> 00:52:25,380 is the heat equation part. 738 00:52:28,860 --> 00:52:35,310 And I can keep this explicit perhaps. 739 00:52:35,310 --> 00:52:39,810 And I could even think of splitting the step. 740 00:52:39,810 --> 00:52:41,370 This is something we haven't done. 741 00:52:41,370 --> 00:52:46,900 I could think of taking a wave step and then a diffusion step. 742 00:52:49,490 --> 00:52:50,320 Whatever. 743 00:52:50,320 --> 00:52:53,500 Putting them together into one step, where the diffusion 744 00:52:53,500 --> 00:53:02,210 step is changed to implicit, because of the great stability 745 00:53:02,210 --> 00:53:05,540 that that provides. 746 00:53:05,540 --> 00:53:06,620 OK. 747 00:53:06,620 --> 00:53:10,370 I think you're seeing the reality 748 00:53:10,370 --> 00:53:22,190 of scientific computation even at this fundamental model 749 00:53:22,190 --> 00:53:30,130 level of a model with two quite different physical effects 750 00:53:30,130 --> 00:53:37,530 of motion with finite speed, energy preserving, magnitude 751 00:53:37,530 --> 00:53:43,140 of G tending to be right at 1; and diffusion 752 00:53:43,140 --> 00:53:52,510 at infinite speed, giving us a growth factor that 753 00:53:52,510 --> 00:54:00,605 tends to be real and below 1, but the price 754 00:54:00,605 --> 00:54:05,020 is it pushes us toward implicit methods. 755 00:54:05,020 --> 00:54:10,490 I think that's really -- this one choice that I've worked 756 00:54:10,490 --> 00:54:14,370 through, and the others that I've mention, 757 00:54:14,370 --> 00:54:16,850 present a real problem. 758 00:54:16,850 --> 00:54:22,980 I'll say another word about this topic next time 759 00:54:22,980 --> 00:54:27,390 and then we'll pretty soon move to non-linear 760 00:54:27,390 --> 00:54:35,980 equations and more realistic, more difficult applications. 761 00:54:35,980 --> 00:54:36,670 OK. 762 00:54:36,670 --> 00:54:37,920 Thanks.