WEBVTT 00:00:00.000 --> 00:00:02.320 The following content is provided 00:00:02.320 --> 00:00:04.910 by MIT OpenCourseWare under a Creative Commons license. 00:00:04.910 --> 00:00:09.510 Additional information about our license and MIT OpenCourseWare, 00:00:09.510 --> 00:00:11.770 in general, is available at ocw.mit.edu. 00:00:16.950 --> 00:00:19.020 PROFESSOR: OK. 00:00:19.020 --> 00:00:27.170 Last time, I started to talk about the heat equation, also 00:00:27.170 --> 00:00:29.830 known as the diffusion equation. 00:00:29.830 --> 00:00:36.790 And what we did then was to solve it. 00:00:36.790 --> 00:00:39.540 Give an analytical solution. 00:00:39.540 --> 00:00:44.180 Both for any initial value u of x and 0, 00:00:44.180 --> 00:00:49.540 and then, specifically for the point source. 00:00:49.540 --> 00:00:52.290 So then, we're talking about the impulse response, 00:00:52.290 --> 00:00:55.160 or the fundamental solution, when 00:00:55.160 --> 00:00:56.960 we start from a delta function. 00:00:56.960 --> 00:01:01.090 And, it was that Gaussian, bell-shaped curve. 00:01:01.090 --> 00:01:04.710 Fantastic fundamental solution. 00:01:04.710 --> 00:01:09.270 But the formulas you get involving 00:01:09.270 --> 00:01:12.700 infinite integrals are not the greatest for actually getting 00:01:12.700 --> 00:01:17.560 numbers for the whole temperature 00:01:17.560 --> 00:01:21.930 history, u of t and x. 00:01:21.930 --> 00:01:25.120 So we have to go to finite differences. 00:01:25.120 --> 00:01:28.600 And, of course, we need finite differences also 00:01:28.600 --> 00:01:38.500 in more general cases, where non-linear terms could appear. 00:01:38.500 --> 00:01:45.260 Or even just variation in x and t in the linear coefficient, 00:01:45.260 --> 00:01:48.360 which I had just set to be 1. 00:01:48.360 --> 00:01:50.520 So, this is the first job. 00:01:50.520 --> 00:01:55.020 Difference methods for the heat equation. 00:01:55.020 --> 00:02:02.590 And you'll see that we get pushed toward implicit methods. 00:02:02.590 --> 00:02:06.400 Because explicit method will require delta t 00:02:06.400 --> 00:02:11.130 to be that very small size delta x squared, 00:02:11.130 --> 00:02:15.130 and that's pretty slow-going. 00:02:15.130 --> 00:02:22.040 And of course, what I'm saying applies equally to -- 00:02:22.040 --> 00:02:28.950 we might be in 2D or in 3D diffusion of pollution, 00:02:28.950 --> 00:02:33.850 for example, in environmental engineering. 00:02:33.850 --> 00:02:36.620 Our fundamental solution from the delta function 00:02:36.620 --> 00:02:42.420 would easily extend to two dimensions or three. 00:02:42.420 --> 00:02:44.820 And, our difference methods will extend, 00:02:44.820 --> 00:02:47.650 but we'll see what comes up. 00:02:47.650 --> 00:02:53.540 And then, let me point to this second problem which 00:02:53.540 --> 00:02:55.360 combines the two. 00:02:55.360 --> 00:02:58.670 We already know about one-way waves, 00:02:58.670 --> 00:03:01.300 coming from a convection term. 00:03:01.300 --> 00:03:07.880 Now we're getting the basic ideas for the diffusion term 00:03:07.880 --> 00:03:14.330 and lot of very important, real problems combine the two. 00:03:14.330 --> 00:03:20.480 And then you have a situation that I'm just 00:03:20.480 --> 00:03:21.490 sort of pointing to. 00:03:21.490 --> 00:03:26.020 Here, that coefficient, as before, 00:03:26.020 --> 00:03:29.400 has the units of velocity. 00:03:29.400 --> 00:03:30.660 Distance over time. 00:03:30.660 --> 00:03:32.600 You see, if I bring time up here, 00:03:32.600 --> 00:03:34.880 I have time over distance there. 00:03:34.880 --> 00:03:38.600 So I need distance over time from c. 00:03:38.600 --> 00:03:42.520 It's what we've seen always. 00:03:42.520 --> 00:03:47.950 And, then the c delta t over delta x, the quantity r, 00:03:47.950 --> 00:03:49.050 was dimensionless. 00:03:49.050 --> 00:03:54.670 So, I'm taking three minutes to speak about dimensions. 00:03:54.670 --> 00:03:57.890 Dimensional analysis is a very simple but useful thing 00:03:57.890 --> 00:03:59.990 to do at the beginning of a problem. 00:03:59.990 --> 00:04:05.290 So here we had r equals c delta t over delta x, which 00:04:05.290 --> 00:04:07.000 is dimensionless. 00:04:07.000 --> 00:04:12.490 And, typically for explicit methods, 00:04:12.490 --> 00:04:16.950 it had to be bounded by 1, for example. 00:04:16.950 --> 00:04:23.570 Now, for the diffusion term d. 00:04:23.570 --> 00:04:25.870 What are the units of d? 00:04:25.870 --> 00:04:29.550 Well, here, now I have distance squared in the denominator, 00:04:29.550 --> 00:04:35.070 so d needs to be a distance squared over time. 00:04:35.070 --> 00:04:43.870 So it has different dimensions, and the quantity capital 00:04:43.870 --> 00:04:47.050 R that we'll meet, parallel to this, 00:04:47.050 --> 00:04:53.890 will be d delta t over delta x squared. 00:04:53.890 --> 00:04:58.790 And, again, it's dimensionless. 00:04:58.790 --> 00:05:02.820 Because if d has dimensions of distance squared over time, 00:05:02.820 --> 00:05:06.200 that reverses it to give dimensionless quantity, 00:05:06.200 --> 00:05:09.680 and we'll see that explicit methods require 00:05:09.680 --> 00:05:16.660 a bound on capital R, and that's what so forces us 00:05:16.660 --> 00:05:18.711 into that small time step. 00:05:18.711 --> 00:05:19.210 OK. 00:05:19.210 --> 00:05:23.900 So you'll see capital R playing a similar role to little r. 00:05:23.900 --> 00:05:26.120 And actually, while we're talking 00:05:26.120 --> 00:05:32.510 about dimensional analysis, back to the differential equation, 00:05:32.510 --> 00:05:35.900 suppose we want, as we certainly do, 00:05:35.900 --> 00:05:40.670 to compare the importance of the convection 00:05:40.670 --> 00:05:45.250 term and the diffusion term. 00:05:45.250 --> 00:05:55.130 Well, that comparison is somehow the ratio of c to d, So, 00:05:55.130 --> 00:05:58.540 somehow, it's natural to think of the ratio of c 00:05:58.540 --> 00:06:06.440 to d as telling us a very important fact. 00:06:06.440 --> 00:06:10.570 Is our problem more like a wave problem? 00:06:10.570 --> 00:06:20.220 Do we expect to be on the unit circle, or close, 00:06:20.220 --> 00:06:25.620 with a term that in Fourier space is pure imaginary, 00:06:25.620 --> 00:06:27.680 from the i*k? 00:06:27.680 --> 00:06:35.730 Or, if d is big, we have then a small ratio, 00:06:35.730 --> 00:06:38.870 and that tells us that the diffusion term is important, 00:06:38.870 --> 00:06:41.650 that we've got lots of dissipation, 00:06:41.650 --> 00:06:46.740 that the stability, at least in the differential equation, 00:06:46.740 --> 00:06:52.640 is much greater, because this is producing a minus k squared 00:06:52.640 --> 00:06:57.720 in Fourier space, and a strongly negative term. 00:06:57.720 --> 00:07:01.260 So it's that ratio c to d but of course, my point 00:07:01.260 --> 00:07:06.500 here in the dimensional stuff is that c to d 00:07:06.500 --> 00:07:12.220 isn't quite OK as it is, because they have different units. 00:07:12.220 --> 00:07:14.410 And we need one more distance. 00:07:14.410 --> 00:07:19.400 So we need the ratio of c -- where do I need an L? 00:07:19.400 --> 00:07:21.390 I need another L up there. 00:07:21.390 --> 00:07:30.210 I need another distance there to get something that 00:07:30.210 --> 00:07:33.620 is altogether dimensionless. 00:07:36.150 --> 00:07:39.980 And so that's the critical number. 00:07:39.980 --> 00:07:42.930 And, it's a known as the Peclet Number. 00:07:42.930 --> 00:07:46.910 Pe, I'll call it. 00:07:46.910 --> 00:07:49.160 What is L? 00:07:49.160 --> 00:07:53.740 We're talking here only about the differential equation. l is 00:07:53.740 --> 00:07:55.580 characteristic length. 00:07:55.580 --> 00:08:01.110 If our problem is set on a finite interval, 00:08:01.110 --> 00:08:03.430 then it's probably the length of the interval. 00:08:03.430 --> 00:08:06.540 It's a characteristic length of the problem. 00:08:06.540 --> 00:08:16.340 And, then this, quantity c*L over d gives us a good measure 00:08:16.340 --> 00:08:18.800 of the importance of these two. 00:08:18.800 --> 00:08:23.580 And, we'll see in the finite difference case, 00:08:23.580 --> 00:08:26.870 there will be cell Peclet number, 00:08:26.870 --> 00:08:28.360 which is entirely different. 00:08:28.360 --> 00:08:33.740 A Peclet number for the little finite different cell, 00:08:33.740 --> 00:08:39.440 where this L is changed to like a delta x. 00:08:39.440 --> 00:08:41.600 OK. 00:08:41.600 --> 00:08:47.640 I think we've got it there. 00:08:47.640 --> 00:08:52.470 Have I got the numerator and denominator correct there, 00:08:52.470 --> 00:08:54.760 or should they be reversed? 00:08:54.760 --> 00:08:56.500 Is it time over distance? 00:09:05.210 --> 00:09:06.410 It is, isn't it? 00:09:06.410 --> 00:09:07.510 Time over distance. 00:09:07.510 --> 00:09:10.500 That shows that I'm not a dimensional analyst. 00:09:10.500 --> 00:09:12.740 Let me fix it. 00:09:12.740 --> 00:09:16.110 And, then that will mean fixing this. 00:09:16.110 --> 00:09:22.310 So, let me just fix that. 00:09:22.310 --> 00:09:25.220 It doesn't sound right as I said it, but I read it wrong. 00:09:25.220 --> 00:09:25.870 OK. 00:09:25.870 --> 00:09:28.770 So c is a time over distance. 00:09:28.770 --> 00:09:30.700 Thanks. 00:09:30.700 --> 00:09:33.400 Is that right? 00:09:33.400 --> 00:09:35.400 It was right the first time? 00:09:35.400 --> 00:09:38.800 In that case, we'll move on. 00:09:38.800 --> 00:09:40.270 OK. 00:09:40.270 --> 00:09:41.240 All right. 00:09:41.240 --> 00:09:42.080 OK. 00:09:42.080 --> 00:09:43.070 Good. 00:09:43.070 --> 00:09:45.190 Thanks. 00:09:45.190 --> 00:09:45.900 Yeah. 00:09:45.900 --> 00:09:47.220 You're right. 00:09:47.220 --> 00:09:49.800 OK. 00:09:49.800 --> 00:09:51.210 Sorry that's on the tapes. 00:09:51.210 --> 00:09:51.710 OK. 00:09:54.260 --> 00:09:58.890 So now let's begin with a natural explicit method 00:09:58.890 --> 00:10:00.620 for the heat equation. 00:10:00.620 --> 00:10:02.660 What everybody would think of. 00:10:02.660 --> 00:10:08.270 So, this will be explicit for the heat equations, 00:10:08.270 --> 00:10:11.080 and I'm in 1D. 00:10:11.080 --> 00:10:13.670 And, what do you think of? 00:10:13.670 --> 00:10:17.090 Naturally, forward difference in time. 00:10:17.090 --> 00:10:19.120 So, can I write it this way? 00:10:19.120 --> 00:10:23.580 Time difference -- I'll use capital U, as always, 00:10:23.580 --> 00:10:28.650 for the finite difference solution -- divided by delta x, 00:10:28.650 --> 00:10:35.510 and I'm doing the heat equation with c equal to 1. 00:10:35.510 --> 00:10:38.270 Just to keep things simple, let's just 00:10:38.270 --> 00:10:43.470 do the heat equation normalized. 00:10:43.470 --> 00:10:47.090 So, here I'm going to take a second difference 00:10:47.090 --> 00:10:51.950 in the x direction of u divided by delta x squared. 00:10:54.630 --> 00:10:55.130 OK. 00:10:58.490 --> 00:10:59.780 What's the accuracy? 00:11:03.000 --> 00:11:06.750 Well, we know the accuracy of each term, 00:11:06.750 --> 00:11:11.290 the centered difference, the error is the discretization 00:11:11.290 --> 00:11:20.040 error, the local error, the truncation error, 00:11:20.040 --> 00:11:25.620 it's also called, will be proportional to -- oh, 00:11:25.620 --> 00:11:30.420 I'm sorry, that's delta t of course. 00:11:30.420 --> 00:11:32.360 OK. 00:11:32.360 --> 00:11:37.010 So, the discretization error in a forward difference, 00:11:37.010 --> 00:11:50.400 when I divide by delta t -- so it's first order in t, 00:11:50.400 --> 00:11:55.040 and centered difference because of the symmetry is second order 00:11:55.040 --> 00:11:55.540 in x. 00:11:58.470 --> 00:11:58.970 OK. 00:12:01.580 --> 00:12:03.640 And this would be another problem 00:12:03.640 --> 00:12:09.830 where I could have begun with a semi-discretization. 00:12:09.830 --> 00:12:13.950 I could have begun with keeping this 00:12:13.950 --> 00:12:19.450 as a derivative, du/dt, equal this centered difference, 00:12:19.450 --> 00:12:25.210 and the discussion would have been entirely parallel to what 00:12:25.210 --> 00:12:27.700 we'll have here. 00:12:27.700 --> 00:12:32.560 So, I've jumped right to the finite difference in x and t. 00:12:32.560 --> 00:12:35.110 So, in other words, what's happening here? 00:12:35.110 --> 00:12:40.730 U_(n+1), U at that center point, n plus 1, of course, 00:12:40.730 --> 00:12:46.730 our little molecule is computing this backward difference, 00:12:46.730 --> 00:12:51.490 and this -- or rather forward in time, 00:12:51.490 --> 00:12:55.310 so this is the time direction, n plus 1 -- 00:12:55.310 --> 00:12:57.550 and this centered difference in space. 00:12:57.550 --> 00:13:03.580 So if I just write out what that does, it's u_(j, 00:13:03.580 --> 00:13:08.950 n) coming from the time difference. 00:13:08.950 --> 00:13:13.330 Then, I'll multiply by delta t, and I have a delta x squared, 00:13:13.330 --> 00:13:18.390 so that ratio is what I'm calling capital R. Capital R, 00:13:18.390 --> 00:13:28.100 to distinguish the diffusion ratio from the wave ratio. 00:13:28.100 --> 00:13:29.170 Times what? 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 00:13:36.570 --> 00:13:41.480 n, plus U_(j-1) at time n. 00:13:41.480 --> 00:13:45.430 We're good at writing that expression down. 00:13:45.430 --> 00:13:50.440 We can follow e to the i*k*x, of course. 00:13:50.440 --> 00:13:52.920 That's always the idea. 00:13:52.920 --> 00:14:01.550 Start with U_(j, 0) as e to the i*k*j delta x, 00:14:01.550 --> 00:14:03.990 and follow it up. 00:14:03.990 --> 00:14:14.200 We're looking for a growth factor G, 00:14:14.200 --> 00:14:16.340 the growth factor in a single step, 00:14:16.340 --> 00:14:20.750 and, of course, the point is, it depends on k. 00:14:24.800 --> 00:14:30.460 At at every step, the exponential is always with us. 00:14:30.460 --> 00:14:35.260 It's the factor that multiplies that exponential that we want, 00:14:35.260 --> 00:14:37.980 and you see what it is, don't you? 00:14:37.980 --> 00:14:41.620 Can you see what G is here? 00:14:41.620 --> 00:14:47.650 So, if I put in that, and then ultimately factor out 00:14:47.650 --> 00:14:51.290 the exponential, it'll be multiplied by a 1, 00:14:51.290 --> 00:14:56.980 from that part, and an R from this, and now in parentheses, 00:14:56.980 --> 00:15:01.995 we have the now very familiar quantity that comes from e 00:15:01.995 --> 00:15:09.270 to the i*k delta x, e to the minus i*k delta x. 00:15:09.270 --> 00:15:14.650 Those combine to give 2 cosine k delta x, 00:15:14.650 --> 00:15:18.590 and the center one is minus 2. 00:15:18.590 --> 00:15:20.060 OK. 00:15:20.060 --> 00:15:22.690 That's our growth factor. 00:15:22.690 --> 00:15:26.180 And the question is, does it stay 00:15:26.180 --> 00:15:32.910 below 1, which is the stability requirement, for all k delta x. 00:15:32.910 --> 00:15:37.860 And, so you can see the dangerous frequency as, 00:15:37.860 --> 00:15:44.030 often, the dangerous frequency is when k delta x is pi, 00:15:44.030 --> 00:15:47.820 and that cosine becomes minus 1. 00:15:47.820 --> 00:15:59.980 And, then we have -- So, stability at k delta x equal pi 00:15:59.980 --> 00:16:06.820 needs -- so this becomes 1 and this is negative 2, 00:16:06.820 --> 00:16:07.890 and that's negative 2. 00:16:07.890 --> 00:16:13.840 We have 1 minus 4R, and that quantity 00:16:13.840 --> 00:16:18.960 has to be, in absolute value, not exceeding 1. 00:16:18.960 --> 00:16:23.800 Well, there's no danger that it's going to go above 1. 00:16:23.800 --> 00:16:27.370 The danger is, will it go below minus 1, 00:16:27.370 --> 00:16:30.070 and it's not allowed to. 00:16:30.070 --> 00:16:34.670 So we need it to be greater than or equal to minus 1. 00:16:34.670 --> 00:16:37.910 And, now, if I put the minus 1 on this side and the 4R 00:16:37.910 --> 00:16:42.890 on that side, that means 2 is greater or equal to 4R. 00:16:42.890 --> 00:16:49.690 And, I get the famous condition R less or equal to half. 00:16:49.690 --> 00:16:52.300 That's the stability requirement. 00:16:55.090 --> 00:16:58.880 And if there was a coefficient d in the heat equation, 00:16:58.880 --> 00:17:03.620 the coefficient would be there in R. OK. 00:17:03.620 --> 00:17:05.440 So, that's the stability condition, 00:17:05.440 --> 00:17:12.710 and that's the thing we would like to find a way around, 00:17:12.710 --> 00:17:17.120 because, it says that this delta t -- 00:17:17.120 --> 00:17:20.390 I'll just put that in, less or equal a half, 00:17:20.390 --> 00:17:25.120 and that means a very small delta t. 00:17:25.120 --> 00:17:26.910 OK. 00:17:26.910 --> 00:17:30.220 What's the way around it? 00:17:30.220 --> 00:17:33.220 To go to an explicit method. 00:17:33.220 --> 00:17:38.660 So, let me go first to a completely explicit method. 00:17:38.660 --> 00:17:39.680 Fully explicit. 00:17:42.390 --> 00:17:43.390 I'm sorry. 00:17:43.390 --> 00:17:46.140 To go to an implicit method. 00:17:46.140 --> 00:17:54.780 Now, implicit method -- now I'm going to do -- 00:17:54.780 --> 00:17:59.090 switch to an implicit method by computing this difference 00:17:59.090 --> 00:18:01.120 at time n plus 1. 00:18:01.120 --> 00:18:03.100 So this is now n plus 1. 00:18:05.890 --> 00:18:09.410 And, what's the change in G? 00:18:09.410 --> 00:18:10.861 A big change. 00:18:10.861 --> 00:18:11.360 OK. 00:18:11.360 --> 00:18:19.960 So now I want to recompute G, for the growth 00:18:19.960 --> 00:18:26.290 factor over a single step for the exponential, when 00:18:26.290 --> 00:18:29.805 the right-hand side, the second difference 00:18:29.805 --> 00:18:31.930 is taken at the new time. 00:18:31.930 --> 00:18:38.130 So this has to move over to the other side of the equation. 00:18:38.130 --> 00:18:40.700 That makes the stability better, but, of course, 00:18:40.700 --> 00:18:46.100 it's going to make the computation more difficult. 00:18:46.100 --> 00:18:48.500 Let me just think about the stability for now. 00:18:48.500 --> 00:18:51.950 So, instead of this board, this is the explicit case. 00:18:57.136 --> 00:18:57.781 Yeah. 00:18:57.781 --> 00:18:58.656 AUDIENCE: [INAUDIBLE] 00:19:05.930 --> 00:19:08.820 PROFESSOR: So it has to be, let's draw a little picture 00:19:08.820 --> 00:19:10.440 here. 00:19:10.440 --> 00:19:11.890 So, there's 0. 00:19:11.890 --> 00:19:12.840 There's 1. 00:19:12.840 --> 00:19:15.760 And here's minus 1. 00:19:15.760 --> 00:19:21.640 And this quantity is starting at 1, and subtracting 4R. 00:19:21.640 --> 00:19:24.280 So, the question is, does it stop here? 00:19:24.280 --> 00:19:30.810 If it stops there, then in this picture, 00:19:30.810 --> 00:19:32.060 these would both be negative. 00:19:32.060 --> 00:19:34.240 This would be about minus a half, 00:19:34.240 --> 00:19:36.800 but that is greater than minus 1. 00:19:36.800 --> 00:19:41.000 So this would be the OK case, when this is true. 00:19:41.000 --> 00:19:44.820 And that's what required R less or equal to half. 00:19:44.820 --> 00:19:50.470 So in this case, yeah. 00:19:50.470 --> 00:19:56.780 In that case, R is less than -- but if R equal to half -- 00:19:56.780 --> 00:19:59.340 everybody sees this -- if R equaled a half, 00:19:59.340 --> 00:20:01.060 this would be 1 minus 2. 00:20:01.060 --> 00:20:03.320 This would bring us all the way here. 00:20:03.320 --> 00:20:07.260 And if R is bigger than a half, then that 00:20:07.260 --> 00:20:12.450 would reverse that inequality, and we would be unstable. 00:20:12.450 --> 00:20:12.950 Yeah. 00:20:12.950 --> 00:20:15.130 I think it's OK. 00:20:15.130 --> 00:20:15.630 Thanks. 00:20:19.040 --> 00:20:22.760 It's the difference between looking at the magnitude -- 00:20:22.760 --> 00:20:27.830 if I put magnitude, then I would need less or equal or just 00:20:27.830 --> 00:20:30.840 realizing that it's just taking it as a number. 00:20:30.840 --> 00:20:32.501 1 minus 4. 00:20:32.501 --> 00:20:33.000 OK. 00:20:33.000 --> 00:20:34.600 I think that's all right. 00:20:34.600 --> 00:20:37.440 Now the implicit case. 00:20:37.440 --> 00:20:39.870 What happened to G? 00:20:39.870 --> 00:20:42.490 What happened to G? 00:20:42.490 --> 00:20:43.820 Let me put it here. 00:20:43.820 --> 00:20:44.400 Implicit. 00:20:50.560 --> 00:20:53.240 So, what's changed? 00:20:53.240 --> 00:20:57.630 This term, which came from the second difference, 00:20:57.630 --> 00:21:02.350 is moved over to the n plus 1 side of the equation. 00:21:02.350 --> 00:21:12.830 So, we have -- in Fourier space, we would have 1 from the u_(j, 00:21:12.830 --> 00:21:21.240 n+1), and then minus 2R, and times our expression 2 cosine k 00:21:21.240 --> 00:21:31.730 delta x minus 2 -- that's what multiplies e to the i -- 00:21:31.730 --> 00:21:39.900 do you see that that's showing up on the new side 00:21:39.900 --> 00:21:42.900 of the equation, multiplying e to the i*k delta x, 00:21:42.900 --> 00:21:57.280 and on the old side, the old time step, we only have the 1. 00:21:57.280 --> 00:21:58.960 In other words, the growth factor 00:21:58.960 --> 00:22:02.150 is, this goes down into the denominator, 00:22:02.150 --> 00:22:04.810 and we can see what's going on. 00:22:04.810 --> 00:22:09.520 So the growth factor is 1 over this expression. 00:22:09.520 --> 00:22:18.910 1 minus 2R times 2 cosine k delta x minus 2. 00:22:18.910 --> 00:22:20.780 We keep running into that quantity 00:22:20.780 --> 00:22:22.670 from the second difference. 00:22:22.670 --> 00:22:23.310 OK. 00:22:23.310 --> 00:22:27.180 So if I look at that, what's going on? 00:22:27.180 --> 00:22:29.005 The key point of this quantity is 00:22:29.005 --> 00:22:35.390 that it's always negative, or zero, but never positive. 00:22:35.390 --> 00:22:42.150 Because this term is never as large as this one. 00:22:42.150 --> 00:22:43.240 So it's negative. 00:22:43.240 --> 00:22:45.160 And now, I have another minus sign here, 00:22:45.160 --> 00:22:48.380 so I have 1 divided by 1 plus something. 00:22:48.380 --> 00:22:56.830 And always, for all k and all ratios R, it's OK. 00:22:56.830 --> 00:22:58.060 Right? 00:22:58.060 --> 00:23:00.460 1 over 1 plus something positive. 00:23:00.460 --> 00:23:02.390 And you see that everything's real. 00:23:02.390 --> 00:23:05.610 That's because we're doing the diffusion equation, where 00:23:05.610 --> 00:23:09.970 the wave equation took us, we had some imaginary part. 00:23:15.130 --> 00:23:17.050 Absolutely stable. 00:23:17.050 --> 00:23:18.490 Stable for all delta t. 00:23:24.320 --> 00:23:26.740 I'll say even for large delta t. 00:23:31.120 --> 00:23:32.090 That's good. 00:23:32.090 --> 00:23:40.970 Of course, we're not going to take delta t equal to 1000. 00:23:40.970 --> 00:23:45.640 Even if stability allowed us, the accuracy would be hopeless. 00:23:45.640 --> 00:23:47.120 It'd be terrible. 00:23:47.120 --> 00:23:54.430 I can't take a giant time step and expect to follow the true 00:23:54.430 --> 00:23:58.330 exponentials for e to the i*k*x. 00:23:58.330 --> 00:24:05.870 But, I can take a larger step than this very small delta, 00:24:05.870 --> 00:24:08.120 one half delta x squared step. 00:24:08.120 --> 00:24:10.060 So this good. 00:24:10.060 --> 00:24:11.430 What's the price that we pay? 00:24:14.320 --> 00:24:19.660 The price is that at every time step, we have to solve. 00:24:19.660 --> 00:24:25.030 This is now on the other side, so this brings a whole matrix 00:24:25.030 --> 00:24:26.720 over. 00:24:26.720 --> 00:24:33.670 So, this implicit case -- when I bring this second difference 00:24:33.670 --> 00:24:41.930 matrix over -- I've used capital K, 00:24:41.930 --> 00:24:46.200 because the second difference matrix there is probably 00:24:46.200 --> 00:24:52.500 the most interesting difference matrix, I'll say. 00:24:52.500 --> 00:24:56.350 One of the most interesting matrices, period. 00:24:56.350 --> 00:25:00.320 It has the minus 2's on the diagonal and 1's 00:25:00.320 --> 00:25:03.430 above the diagonal. 00:25:03.430 --> 00:25:06.950 So, actually, with signs exchanged, 00:25:06.950 --> 00:25:10.570 I called it K. So can I just remember 00:25:10.570 --> 00:25:15.210 the meaning of this matrix K. I'm 00:25:15.210 --> 00:25:19.750 going to stay with that letter because it's 00:25:19.750 --> 00:25:24.220 standard in finite elements for the stiffness matrix. 00:25:24.220 --> 00:25:27.480 So, K, by definition, is the matrix with the 2's 00:25:27.480 --> 00:25:32.940 on the diagonal, minus 1's above the diagonal, and minus 1's 00:25:32.940 --> 00:25:35.170 below the diagonal. 00:25:35.170 --> 00:25:37.340 And, what size is it? 00:25:40.100 --> 00:25:46.130 On the whole line, where we're working for simplicity, 00:25:46.130 --> 00:25:51.850 I guess K is infinite in both directions. 00:25:51.850 --> 00:25:54.000 2's and minus 1's. 00:25:54.000 --> 00:25:58.910 On an interval where we would really compute, 00:25:58.910 --> 00:26:01.240 K would be a finite matrix. 00:26:01.240 --> 00:26:04.395 But, the main point is that this matrix K is minus 00:26:04.395 --> 00:26:06.480 the second difference matrix. 00:26:09.460 --> 00:26:14.900 So that when I bring it over to this side, which 00:26:14.900 --> 00:26:20.750 changes that plus to a minus, I get this equation to solve. 00:26:20.750 --> 00:26:31.100 The identity, coming from the 1, plus R, the ratio, 00:26:31.100 --> 00:26:34.800 times this matrix K, multiplies this vector 00:26:34.800 --> 00:26:38.900 of all the values at the new time step, 00:26:38.900 --> 00:26:45.590 and equals, in this case, all the values at the old time 00:26:45.590 --> 00:26:46.850 step. 00:26:46.850 --> 00:26:51.370 Because, at the old time step, I just have the identity. 00:26:51.370 --> 00:26:52.860 Are you OK with this? 00:26:52.860 --> 00:26:58.010 This is the problem that we actually have to solve. 00:26:58.010 --> 00:27:00.430 I mean, the main point about it is 00:27:00.430 --> 00:27:03.910 that we have a linear system. 00:27:03.910 --> 00:27:07.440 The u's at the new time step are coupled. 00:27:07.440 --> 00:27:09.430 Coupled by the second difference. 00:27:09.430 --> 00:27:12.830 Coupled by the matrix K, so we have to deal with it. 00:27:12.830 --> 00:27:17.200 Well actually, in one dimension, no problem. 00:27:17.200 --> 00:27:22.980 Because in one dimension, K is just a tri-diagonal matrix, 00:27:22.980 --> 00:27:29.360 and we can solve tri-diagonal systems 00:27:29.360 --> 00:27:31.930 almost as fast as we could work with the identity matrix. 00:27:34.870 --> 00:27:40.820 And actually, what we're seeing in this ratio 00:27:40.820 --> 00:27:46.850 are exactly the eigenvalues of the inverse matrix. 00:27:46.850 --> 00:27:51.910 So, if I use a formula, I would say U_(n+1) is the inverse 00:27:51.910 --> 00:27:58.910 matrix times U_n, and the growth factor is just the eigenvalues. 00:27:58.910 --> 00:28:02.490 The growth factor G is simply the eigenvalues 00:28:02.490 --> 00:28:05.070 of this inverse matrix. 00:28:05.070 --> 00:28:12.800 And, notice again, that all the eigenvalues of K 00:28:12.800 --> 00:28:15.350 are greater or equal zero. 00:28:15.350 --> 00:28:17.840 It was a positive definite matrix, 00:28:17.840 --> 00:28:20.210 because I switched sign. 00:28:20.210 --> 00:28:22.940 And, the effect of the identity is 00:28:22.940 --> 00:28:25.760 to push it a little more positive. 00:28:25.760 --> 00:28:26.730 Here. 00:28:26.730 --> 00:28:31.550 It's making -- this is a positive number there. 00:28:31.550 --> 00:28:34.460 That's the eigenvalues of R*K, and then, 00:28:34.460 --> 00:28:37.100 I making it a little more positive with this 1, 00:28:37.100 --> 00:28:49.910 so we have a problem that is safely invertible, stable, 00:28:49.910 --> 00:28:52.330 and in one dimension, quite simple. 00:28:52.330 --> 00:29:00.900 So in one dimension we're entirely ready to go implicit. 00:29:00.900 --> 00:29:04.570 Let me just, while I'm thinking about it, 00:29:04.570 --> 00:29:09.630 recall that we had some other matrix in the finite matrix 00:29:09.630 --> 00:29:19.270 case, we had some other matrices in the very first section 00:29:19.270 --> 00:29:25.350 of the original notes, which had different boundary conditions. 00:29:25.350 --> 00:29:32.110 In other words, if we're on a finite interval, 00:29:32.110 --> 00:29:37.560 then this matrix K corresponds to zero boundary conditions. 00:29:37.560 --> 00:29:39.790 Zero temperature at the end. 00:29:39.790 --> 00:29:43.970 But we might have other boundary conditions, 00:29:43.970 --> 00:29:46.310 like the derivative of the temperature 00:29:46.310 --> 00:29:49.220 is zero at one end, or the other end. 00:29:49.220 --> 00:29:54.450 That would change this matrix in the first row and the last row 00:29:54.450 --> 00:29:55.880 to another nice matrix. 00:29:55.880 --> 00:29:59.520 So, we might have not K but one of the other tri-diagonal 00:29:59.520 --> 00:30:02.630 matrices. 00:30:02.630 --> 00:30:06.660 But, here's the real point. 00:30:06.660 --> 00:30:08.780 What happens in two dimensions? 00:30:08.780 --> 00:30:18.250 Suppose I include now a u_yy term in the heat equation, 00:30:18.250 --> 00:30:21.470 or three dimensions with a u_zz term. 00:30:21.470 --> 00:30:25.280 What's different? 00:30:25.280 --> 00:30:29.260 Well, we can follow through the stability 00:30:29.260 --> 00:30:33.840 test for the explicit method, and it'll all 00:30:33.840 --> 00:30:37.640 have a term from delta x squared, 00:30:37.640 --> 00:30:40.290 and a term from delta y squared, and a delta z squared. 00:30:40.290 --> 00:30:43.940 It'll be even slightly tighter on delta t. 00:30:43.940 --> 00:30:48.740 Or I can follow through the implicit method, where 00:30:48.740 --> 00:30:51.900 now this second difference in the x direction 00:30:51.900 --> 00:30:54.940 appears also with a second difference in the y direction 00:30:54.940 --> 00:30:58.370 and a second different in the z direction. 00:30:58.370 --> 00:31:02.050 So what's happening? 00:31:02.050 --> 00:31:03.990 From the matrix theory point of view, 00:31:03.990 --> 00:31:08.920 I still have nice matrices with 2's and minus 1's or maybe with 00:31:08.920 --> 00:31:12.830 -- I guess what we would have -- can we just remember what we 00:31:12.830 --> 00:31:16.950 would have in the two-dimensional case. 00:31:16.950 --> 00:31:19.150 In the two-dimensional case, instead of 00:31:19.150 --> 00:31:24.700 this second difference with a minus 1, a 2, and a minus 1 00:31:24.700 --> 00:31:29.150 in the x direction, we also have a second difference 00:31:29.150 --> 00:31:34.260 coming from the u_yy in the y direction, so another minus 1, 00:31:34.260 --> 00:31:38.480 and that 2 moves up to a 4. 00:31:38.480 --> 00:31:43.840 So our matrix in the 2D, two space dimension case, 00:31:43.840 --> 00:31:47.520 will have a 4's on the diagonal. 00:31:47.520 --> 00:31:50.000 So it grows from this one. 00:31:50.000 --> 00:31:55.230 It's called the tensor product coming out of this. 00:31:55.230 --> 00:31:59.910 Anyway, it produces this famous five-point molecule 00:31:59.910 --> 00:32:04.220 with 4's on the diagonal, and four minus 1's 00:32:04.220 --> 00:32:06.960 going down other diagonals. 00:32:06.960 --> 00:32:12.930 And the problem is that those diagonal are not -- 00:32:12.930 --> 00:32:15.430 two of the diagonals might be right next to the main 00:32:15.430 --> 00:32:19.570 diagonal, but the other two will be farther away. 00:32:19.570 --> 00:32:22.240 Because we can't number the nodes to keep all 00:32:22.240 --> 00:32:28.000 these four numbered adjacently at all the mesh points. 00:32:28.000 --> 00:32:30.510 In other words, we have a serious matrix problem. 00:32:30.510 --> 00:32:32.490 That's my comment. 00:32:32.490 --> 00:32:38.670 We have a matrix K in 2D. 00:32:38.670 --> 00:32:41.820 So what I'm speaking about is this problem. 00:32:41.820 --> 00:32:49.980 Identity plus R*K(2D) u at the new time equal u at the old 00:32:49.980 --> 00:32:50.810 time. 00:32:50.810 --> 00:32:54.370 And K(2D) is by no means so easy to deal with as K(1D). 00:32:58.910 --> 00:33:02.770 It doesn't have a nice narrow band like this. 00:33:02.770 --> 00:33:05.010 The band stretches further, and we'll 00:33:05.010 --> 00:33:11.120 see in spades, in then next part of the course which 00:33:11.120 --> 00:33:17.410 is about solving large linear systems. 00:33:17.410 --> 00:33:19.150 So we have a large system because we've 00:33:19.150 --> 00:33:21.560 got lots of mesh points. 00:33:21.560 --> 00:33:25.840 We have matrix that has a wider band, 00:33:25.840 --> 00:33:28.850 and so you need new ideas in solving it. 00:33:28.850 --> 00:33:36.160 So this problem is what will be a model -- is the connection, 00:33:36.160 --> 00:33:40.000 really, between this section of the course, 00:33:40.000 --> 00:33:42.380 difference methods for initial-value problems, 00:33:42.380 --> 00:33:47.050 and the next section, which is solving large linear systems. 00:33:47.050 --> 00:33:53.570 Because this will be a typical model large linear system. 00:33:53.570 --> 00:33:59.500 So that's a comment that for implicit equations, 00:33:59.500 --> 00:34:04.130 the price was low in 1D, but the price is not so low in 2D. 00:34:04.130 --> 00:34:05.040 OK. 00:34:05.040 --> 00:34:10.210 Now, have we got any other candidates 00:34:10.210 --> 00:34:14.440 between explicit and implicit? 00:34:14.440 --> 00:34:20.810 Well, we did, way back, for ordinary differential 00:34:20.810 --> 00:34:22.510 equations, we had a trapezoidal rule. 00:34:22.510 --> 00:34:25.170 And what did that do? 00:34:25.170 --> 00:34:29.530 That increased the accuracy by centering it 00:34:29.530 --> 00:34:31.190 at time n plus a half. 00:34:31.190 --> 00:34:32.810 Do you remember that? 00:34:32.810 --> 00:34:35.160 It was half on the right-hand side; 00:34:35.160 --> 00:34:38.830 we used half the explicit difference 00:34:38.830 --> 00:34:40.590 and half of the implicit difference, 00:34:40.590 --> 00:34:44.020 so natural to do it again. 00:34:44.020 --> 00:34:47.140 It worked then, and it'll work now. 00:34:47.140 --> 00:34:49.570 So this will be the trapezoidal method. 00:34:49.570 --> 00:34:57.310 But most people in the diffusion application 00:34:57.310 --> 00:35:03.320 name it after these authors, Crank and Nicolson. 00:35:03.320 --> 00:35:06.560 So it's the Crank-Nicolson method. 00:35:06.560 --> 00:35:12.990 It's half the explicit part plus half the implicit part. 00:35:12.990 --> 00:35:17.390 And the order of accuracy is increased 00:35:17.390 --> 00:35:19.420 by centering it this way. 00:35:19.420 --> 00:35:22.520 That's the point, of course. 00:35:22.520 --> 00:35:25.320 So I guess the point will be we don't 00:35:25.320 --> 00:35:33.560 have any extra work compared to a fully implicit one, 00:35:33.560 --> 00:35:37.120 and we get an extra order of accuracy, 00:35:37.120 --> 00:35:41.350 by centering it at a half. 00:35:41.350 --> 00:35:55.020 So u_(j, n+1) minus u_(j, n) over delta t is R, 00:35:55.020 --> 00:35:59.320 that ratio -- oh sorry. 00:35:59.320 --> 00:36:07.220 We've got delta t there -- so will be the second difference 00:36:07.220 --> 00:36:24.210 at the old time times a half -- it would just take half of that 00:36:24.210 --> 00:36:29.040 explicit part plus half of the implicit part -- 00:36:29.040 --> 00:36:34.830 this is at time n plus one, over delta x squared. 00:36:34.830 --> 00:36:36.050 You would think of that. 00:36:36.050 --> 00:36:37.730 Everybody would. 00:36:37.730 --> 00:36:39.830 It's just a natural idea. 00:36:39.830 --> 00:36:43.820 It increases the accuracy. 00:36:43.820 --> 00:36:45.100 It's still implicit. 00:36:45.100 --> 00:36:50.050 We still have to bring this term over onto this side when u_(j, 00:36:50.050 --> 00:36:53.850 n) moves over onto that side. 00:36:53.850 --> 00:36:56.470 And we can look at the growth factor. 00:36:56.470 --> 00:36:58.870 So that's what we have to do. 00:36:58.870 --> 00:37:04.140 So bring that term over, and bring delta t up, 00:37:04.140 --> 00:37:08.750 so it's delta t divided by delta x squared and R. OK. 00:37:08.750 --> 00:37:10.790 So I've got this term. 00:37:10.790 --> 00:37:12.060 The 1. 00:37:12.060 --> 00:37:13.840 And when I bring that term over, it's 00:37:13.840 --> 00:37:17.340 a half of what we brought over before, 00:37:17.340 --> 00:37:22.790 which was the R, 2 cosine k delta x minus the 2. 00:37:25.400 --> 00:37:29.750 That multiplies G. And on the right side, 00:37:29.750 --> 00:37:35.450 we have the 1 from the u_(j, n) plus the half, 00:37:35.450 --> 00:37:41.080 times the R times 2 cosine k delta x minus 2. 00:37:47.450 --> 00:37:50.720 I'm certainly getting to the point where I write these 00:37:50.720 --> 00:37:57.260 expressions down faster without explicitly writing e to the i*k 00:37:57.260 --> 00:38:00.650 delta x and then canceling it. 00:38:00.650 --> 00:38:05.650 So, G is a ratio again. 00:38:05.650 --> 00:38:14.480 In the explicit case G was all in the numerator. 00:38:14.480 --> 00:38:16.490 There was no denominator. 00:38:16.490 --> 00:38:24.480 In the fully implicit case, g was all denominator, 00:38:24.480 --> 00:38:27.460 with just a one in the numerator. 00:38:27.460 --> 00:38:30.250 Now we have G as the ratio. 00:38:30.250 --> 00:38:32.660 So I'm going to bring this down below. 00:38:39.000 --> 00:38:40.850 And, G is a ratio. 00:38:40.850 --> 00:38:43.360 So it's a fraction because it's got 00:38:43.360 --> 00:38:48.810 an explicit part in the numerator and an implicit part 00:38:48.810 --> 00:38:49.690 in the denominator. 00:38:49.690 --> 00:38:51.400 And, of course, the main question 00:38:51.400 --> 00:38:54.090 is, when is G less or equal to one? 00:38:59.570 --> 00:39:00.640 Well. 00:39:00.640 --> 00:39:07.260 Actually it's cool, because G is simply -- look, 00:39:07.260 --> 00:39:11.250 the numerator is one, and all this stuff is negative. 00:39:11.250 --> 00:39:15.070 So this is one minus something. 00:39:15.070 --> 00:39:18.630 And, what's going down here? 00:39:18.630 --> 00:39:20.020 The signs just reverse. 00:39:20.020 --> 00:39:22.270 It's one plus something. 00:39:22.270 --> 00:39:25.710 And this is a positive quantity, and this 00:39:25.710 --> 00:39:28.450 is the same positive quantity. 00:39:28.450 --> 00:39:31.590 So, this is one minus -- do you see this -- 00:39:31.590 --> 00:39:34.955 because this is negative, I'll write it as 1 minus something 00:39:34.955 --> 00:39:36.660 -- something positive. 00:39:36.660 --> 00:39:40.700 And here, this is minus and this is negative, 00:39:40.700 --> 00:39:43.250 so this is 1 plus the same thing. 00:39:43.250 --> 00:39:44.180 So what's the answer? 00:39:51.830 --> 00:39:53.710 [UNINTELLIGIBLE PHRASE] 00:39:53.710 --> 00:39:57.390 One minus a positive number over 1 plus that positive number 00:39:57.390 --> 00:39:58.390 can't go wrong. 00:40:02.410 --> 00:40:09.230 So, for all R. OK. 00:40:09.230 --> 00:40:17.440 In a way, that has to be more attractive than fully implicit, 00:40:17.440 --> 00:40:19.340 because it didn't take more work, 00:40:19.340 --> 00:40:23.400 it wasn't any less stable, and it was more accurate. 00:40:23.400 --> 00:40:28.740 So Crank-Nicolson is sort of a natural idea 00:40:28.740 --> 00:40:33.420 at this level of analysis. 00:40:36.080 --> 00:40:36.770 OK. 00:40:36.770 --> 00:40:39.970 So those are the three methods to think about. 00:40:39.970 --> 00:40:49.720 Of course, actual computations will bring in new questions. 00:40:49.720 --> 00:40:51.906 Boundary conditions, non-linearities, 00:40:51.906 --> 00:40:52.530 and everything. 00:40:52.530 --> 00:40:55.280 I'm just staying at this simple level, 00:40:55.280 --> 00:41:00.850 where I'm comparing the first ideas that would occur to us. 00:41:00.850 --> 00:41:01.680 OK. 00:41:01.680 --> 00:41:05.250 We may want higher accuracy, all sorts of things. 00:41:05.250 --> 00:41:12.780 But, we're seeing the contrast with the wave equation case. 00:41:12.780 --> 00:41:19.060 But, here we move more toward the implicit side, 00:41:19.060 --> 00:41:25.930 and the reward is unconditional stability, no condition on r, 00:41:25.930 --> 00:41:29.100 but the price is more computation. 00:41:29.100 --> 00:41:31.080 OK. 00:41:31.080 --> 00:41:35.770 Well, at that same level of thinking about the most 00:41:35.770 --> 00:41:41.040 natural methods that occur to you, 00:41:41.040 --> 00:41:45.860 I want to speak now about the big, big issue 00:41:45.860 --> 00:41:47.430 of convection-diffusion. 00:41:47.430 --> 00:41:52.370 Suppose both types of problem are with us. 00:41:52.370 --> 00:41:55.320 So convention-diffusion. 00:41:55.320 --> 00:41:56.340 What would you do then? 00:42:01.520 --> 00:42:11.110 Now, I have both convection and diffusion in the equation u_t 00:42:11.110 --> 00:42:14.650 equals c*u_x plus d*u_xx. 00:42:18.650 --> 00:42:22.920 And I guess I want to say that this is something 00:42:22.920 --> 00:42:25.700 that the world has not come to a decision on. 00:42:25.700 --> 00:42:29.720 The computational world is still debating 00:42:29.720 --> 00:42:32.290 what to do with convection-diffusion. 00:42:32.290 --> 00:42:43.740 So it's a very appropriate equation model to experiment 00:42:43.740 --> 00:42:44.240 with. 00:42:47.240 --> 00:42:53.500 Maybe I'll ask you, as the next, kind of informal homework, 00:42:53.500 --> 00:42:58.660 to begin to experiment with that. 00:42:58.660 --> 00:43:02.050 And when I say experiment, you can think immediately 00:43:02.050 --> 00:43:07.700 of the methods you might think of. 00:43:07.700 --> 00:43:09.850 The first methods to try. 00:43:09.850 --> 00:43:12.290 So the first method of course, will be explicit. 00:43:16.580 --> 00:43:19.850 Fully explicit. 00:43:19.850 --> 00:43:22.110 What do you expect to have happen there? 00:43:22.110 --> 00:43:25.690 Let's just see how the two different terms 00:43:25.690 --> 00:43:33.100 each give their little push towards the stability 00:43:33.100 --> 00:43:34.460 requirement. 00:43:34.460 --> 00:43:44.250 So explicit is going to be the time difference, forward 00:43:44.250 --> 00:43:49.970 time difference equals c times-- what shall I choose first 00:43:49.970 --> 00:43:52.690 for the space difference? 00:43:55.230 --> 00:43:57.030 Upwind, maybe? 00:43:57.030 --> 00:43:59.330 Upwind in space? 00:43:59.330 --> 00:44:04.540 So, this is going to be explicit and upwind. 00:44:04.540 --> 00:44:09.330 Of course, everybody knows that refers to the windy term 00:44:09.330 --> 00:44:13.441 in the equation, the wave term. 00:44:13.441 --> 00:44:13.940 OK. 00:44:13.940 --> 00:44:17.510 So, of course, I'll get a forward difference, then, 00:44:17.510 --> 00:44:23.150 in the x direction of u at time n divided delta x, 00:44:23.150 --> 00:44:27.080 and then I'm going to do the natural one. 00:44:27.080 --> 00:44:35.120 A centered second difference at time n for the diffusion term. 00:44:35.120 --> 00:44:35.620 OK. 00:44:39.660 --> 00:44:43.040 We'll see what comes out of that. 00:44:43.040 --> 00:44:45.510 It's going to involve not only the ratio, 00:44:45.510 --> 00:44:48.900 it will involve both ratios, naturally. 00:44:48.900 --> 00:44:52.890 Little r, as soon as I multiply up by delta x, 00:44:52.890 --> 00:44:57.830 I have a little r there, my wave equation ratio. 00:44:57.830 --> 00:45:00.520 And I have a d delta t over delta x squared, 00:45:00.520 --> 00:45:06.710 which I'm calling capital R. That's the diffusion ratio. 00:45:06.710 --> 00:45:11.510 And multiplied by all sorts of stuff. 00:45:11.510 --> 00:45:15.330 So actually, why don't I write down what the equation really 00:45:15.330 --> 00:45:16.000 is then? 00:45:16.000 --> 00:45:21.830 Explicitly, it's u_(j, n+1) is u_(j, n) -- 00:45:21.830 --> 00:45:27.220 so maybe I should collect all the terms here. 00:45:27.220 --> 00:45:33.360 Can I collect the terms that multiply u_(j+1, n)? 00:45:38.080 --> 00:45:40.780 So what's going to multiply u_(j+1, n)? 00:45:40.780 --> 00:45:43.870 Here I have a c delta t over delta x, 00:45:43.870 --> 00:45:46.690 I have an r, from this term. 00:45:46.690 --> 00:45:49.860 And here I have a d. 00:45:49.860 --> 00:45:50.360 Right? 00:45:53.720 --> 00:45:54.760 No. 00:45:54.760 --> 00:45:55.260 Sorry. 00:45:55.260 --> 00:46:00.100 I have a capital R. Isn't that right? 00:46:00.100 --> 00:46:03.620 When I multiply through by delta t, I have d delta t over delta 00:46:03.620 --> 00:46:04.840 x squared. 00:46:04.840 --> 00:46:08.660 I have a capital R that's going to multiply the -- 00:46:08.660 --> 00:46:11.690 does that look right? 00:46:14.580 --> 00:46:18.900 And then what multiplies the center guy, the u_(j, n)? 00:46:23.980 --> 00:46:29.090 Well that's where the time difference comes over. 00:46:29.090 --> 00:46:31.480 There's the usual 1. 00:46:31.480 --> 00:46:36.760 Then from this difference, this upwind difference, 00:46:36.760 --> 00:46:38.810 there's a minus at the center term, 00:46:38.810 --> 00:46:41.110 because it's a forward difference. 00:46:41.110 --> 00:46:48.160 So that will be a minus r, and from this centered difference, 00:46:48.160 --> 00:46:51.270 well remember, that that's a second difference, 00:46:51.270 --> 00:46:57.660 so there's a coefficient 2 in the middle term, or minus 2, 00:46:57.660 --> 00:47:05.170 and when I bring the delta t up, it's minus 2 capital R. 00:47:05.170 --> 00:47:10.150 And then, the third one, whatever this last coefficient 00:47:10.150 --> 00:47:16.220 is, is multiplying u_(j, n-1), and now, where does that come? 00:47:16.220 --> 00:47:18.010 It doesn't come from here, it doesn't 00:47:18.010 --> 00:47:19.530 come from the upwind difference, I 00:47:19.530 --> 00:47:22.490 guess it only comes from there. 00:47:22.490 --> 00:47:30.220 And it's probably just R. And I guess 00:47:30.220 --> 00:47:35.730 I would check any formula like that, 00:47:35.730 --> 00:47:40.000 to be sure that the coefficients add to 1. 00:47:43.670 --> 00:47:45.700 And they do. 00:47:45.700 --> 00:47:46.990 And why should they? 00:47:46.990 --> 00:47:50.910 Why should those coefficients add to 1? 00:47:50.910 --> 00:47:57.940 So that's my little test that I've got it right. 00:47:57.940 --> 00:48:04.945 They should add to 1, because a constant initial -- if u_(j, 00:48:04.945 --> 00:48:09.490 0) is a constant, if I'm starting from a constant 00:48:09.490 --> 00:48:15.330 function, its x derivatives would all be zero. 00:48:15.330 --> 00:48:18.757 If I had u equal constant, then the x derivatives 00:48:18.757 --> 00:48:20.340 would be zero, so the time derivatives 00:48:20.340 --> 00:48:24.380 would be zero, so that constant would stay forever. 00:48:24.380 --> 00:48:27.790 And so I would want to get the same constant out 00:48:27.790 --> 00:48:31.040 from constants going in. 00:48:31.040 --> 00:48:35.480 And so the coefficient should add to one. 00:48:35.480 --> 00:48:36.190 Well. 00:48:40.830 --> 00:48:45.130 I can see a situation when I know stability will be OK. 00:48:48.420 --> 00:48:53.520 Stability can't fail if those three coefficients are all 00:48:53.520 --> 00:48:55.380 positive. 00:48:55.380 --> 00:48:59.610 If those three coefficients are all positive, 00:48:59.610 --> 00:49:04.760 then my growth factor can't be larger than 1. 00:49:04.760 --> 00:49:07.030 Those coefficients, they add to 1. 00:49:07.030 --> 00:49:11.600 So the growth factor is this times -- 00:49:11.600 --> 00:49:12.990 should I write it down? 00:49:12.990 --> 00:49:15.950 The growth factor will be this -- 00:49:15.950 --> 00:49:17.450 so this will be the growth factor -- 00:49:17.450 --> 00:49:22.760 G will be that number times e to the i*k delta x, 00:49:22.760 --> 00:49:25.530 from shifting over one. 00:49:25.530 --> 00:49:31.780 It'll be this number times 1, and it'll be this number times 00:49:31.780 --> 00:49:37.970 e to the minus i*k delta x, and my point is just that 00:49:37.970 --> 00:49:42.210 the easiest estimate would say, well this has magnitude 1, 00:49:42.210 --> 00:49:45.840 that has magnitude 1, that has magnitude 1. 00:49:45.840 --> 00:49:49.560 If I just add up magnitudes, I'll 00:49:49.560 --> 00:49:53.570 get this one, and this one and this one. 00:49:53.570 --> 00:49:57.520 And if they're all positive, they add to 1, I'm safe. 00:49:57.520 --> 00:50:00.130 So what's the stability condition 00:50:00.130 --> 00:50:04.420 that this simple argument requires? 00:50:04.420 --> 00:50:05.770 Well. 00:50:05.770 --> 00:50:06.900 That's positive. 00:50:06.900 --> 00:50:07.930 That's positive. 00:50:07.930 --> 00:50:09.440 The danger is this one. 00:50:13.040 --> 00:50:16.340 So I guess I would like that to be positive. 00:50:16.340 --> 00:50:23.250 So I would want -- I'm certainly stable if this is positive, 00:50:23.250 --> 00:50:30.971 which I might as well write as r plus 2 big R being not larger 00:50:30.971 --> 00:50:31.470 than 1. 00:50:37.590 --> 00:50:41.050 Because if that's true, I'm taking a combination 00:50:41.050 --> 00:50:45.030 with positive coefficients adding to 1, I can't go wrong. 00:50:45.030 --> 00:50:45.530 OK. 00:50:45.530 --> 00:50:51.290 And in fact, that's pretty familiar. 00:50:55.450 --> 00:50:58.950 It combines the two conditions we know. 00:50:58.950 --> 00:51:01.200 The condition that r should be less or equal 00:51:01.200 --> 00:51:03.840 1, small r, the wave Courant-Friedrichs-Lewy 00:51:03.840 --> 00:51:06.300 condition. 00:51:06.300 --> 00:51:10.260 And the condition that capital R should 00:51:10.260 --> 00:51:14.010 be less or equal to half the diffusion condition. 00:51:14.010 --> 00:51:19.440 Now they're both -- now I'm forcing them even more 00:51:19.440 --> 00:51:23.251 by saying that the sum has to be less or equal 1. 00:51:23.251 --> 00:51:23.750 OK. 00:51:26.720 --> 00:51:35.600 So let me leave the explicit upwind at this point, 00:51:35.600 --> 00:51:40.270 and ask you what other method might you want to try? 00:51:40.270 --> 00:51:43.100 Well, it would be natural to think 00:51:43.100 --> 00:51:51.330 of explicit centered, where you would center this 00:51:51.330 --> 00:51:54.120 and we could follow that one through. 00:51:56.680 --> 00:51:58.750 I guess I'll do that next time very quickly, 00:51:58.750 --> 00:52:01.340 and the notes will do it. 00:52:01.340 --> 00:52:05.650 And I'm not sure which of those two ways is better. 00:52:05.650 --> 00:52:11.080 Or, you're going to say go implicit. 00:52:11.080 --> 00:52:14.520 And now I'm going to ask, all right, but what part of it 00:52:14.520 --> 00:52:15.890 do you want me to make implicit? 00:52:19.720 --> 00:52:22.000 Maybe all I need to make implicit 00:52:22.000 --> 00:52:25.380 is the heat equation part. 00:52:28.860 --> 00:52:35.310 And I can keep this explicit perhaps. 00:52:35.310 --> 00:52:39.810 And I could even think of splitting the step. 00:52:39.810 --> 00:52:41.370 This is something we haven't done. 00:52:41.370 --> 00:52:46.900 I could think of taking a wave step and then a diffusion step. 00:52:49.490 --> 00:52:50.320 Whatever. 00:52:50.320 --> 00:52:53.500 Putting them together into one step, where the diffusion 00:52:53.500 --> 00:53:02.210 step is changed to implicit, because of the great stability 00:53:02.210 --> 00:53:05.540 that that provides. 00:53:05.540 --> 00:53:06.620 OK. 00:53:06.620 --> 00:53:10.370 I think you're seeing the reality 00:53:10.370 --> 00:53:22.190 of scientific computation even at this fundamental model 00:53:22.190 --> 00:53:30.130 level of a model with two quite different physical effects 00:53:30.130 --> 00:53:37.530 of motion with finite speed, energy preserving, magnitude 00:53:37.530 --> 00:53:43.140 of G tending to be right at 1; and diffusion 00:53:43.140 --> 00:53:52.510 at infinite speed, giving us a growth factor that 00:53:52.510 --> 00:54:00.605 tends to be real and below 1, but the price 00:54:00.605 --> 00:54:05.020 is it pushes us toward implicit methods. 00:54:05.020 --> 00:54:10.490 I think that's really -- this one choice that I've worked 00:54:10.490 --> 00:54:14.370 through, and the others that I've mention, 00:54:14.370 --> 00:54:16.850 present a real problem. 00:54:16.850 --> 00:54:22.980 I'll say another word about this topic next time 00:54:22.980 --> 00:54:27.390 and then we'll pretty soon move to non-linear 00:54:27.390 --> 00:54:35.980 equations and more realistic, more difficult applications. 00:54:35.980 --> 00:54:36.670 OK. 00:54:36.670 --> 00:54:37.920 Thanks.