WEBVTT 00:00:01.540 --> 00:00:03.910 The following content is provided under a Creative 00:00:03.910 --> 00:00:05.300 Commons license. 00:00:05.300 --> 00:00:07.510 Your support will help MIT OpenCourseWare 00:00:07.510 --> 00:00:11.600 continue to offer high quality educational resources for free. 00:00:11.600 --> 00:00:14.140 To make a donation, or to view additional materials 00:00:14.140 --> 00:00:18.100 from hundreds of MIT courses, visit MIT OpenCourseWare 00:00:18.100 --> 00:00:18.980 at ocw.mit.edu. 00:00:23.125 --> 00:00:24.750 PROFESSOR: So today we're going to talk 00:00:24.750 --> 00:00:26.160 about partial differential equations. 00:00:26.160 --> 00:00:28.035 Let me just give a quick review where we are. 00:00:28.035 --> 00:00:31.640 So we're going to have a PDE lecture today. 00:00:31.640 --> 00:00:36.450 Then Professor Swan is going to do a review on Friday. 00:00:36.450 --> 00:00:40.980 There's no lecture on Monday, but you have the quiz, too, 00:00:40.980 --> 00:00:42.780 Monday night. 00:00:42.780 --> 00:00:46.452 Wednesday, next week is Veterans Day so it's a holiday. 00:00:46.452 --> 00:00:47.535 There's no classes at MIT. 00:00:50.042 --> 00:00:51.750 So the next time I'll be lecturing to you 00:00:51.750 --> 00:00:55.410 will be a week from Friday. 00:00:55.410 --> 00:00:58.340 The homework following the quiz will be posted 00:00:58.340 --> 00:01:01.620 and so you guys can get started on that. 00:01:01.620 --> 00:01:04.140 That homework involves a COMSOL problem 00:01:04.140 --> 00:01:06.480 which, for those of you who had trouble with COMSOL, 00:01:06.480 --> 00:01:07.680 might take some extra time, since it's 00:01:07.680 --> 00:01:08.971 the first time you're doing it. 00:01:08.971 --> 00:01:11.940 So I would suggest you at least give it a try early. 00:01:11.940 --> 00:01:14.010 I'm sure the TAs would be very happy to help you 00:01:14.010 --> 00:01:17.190 to force COMSOL to cooperate. 00:01:17.190 --> 00:01:18.540 Bend COMSOL to your will. 00:01:21.270 --> 00:01:24.990 All right, so we're going to talk about partial differential 00:01:24.990 --> 00:01:25.800 equations. 00:01:25.800 --> 00:01:28.002 And, as it sounds, what's different 00:01:28.002 --> 00:01:29.710 about a partial differential equation is, 00:01:29.710 --> 00:01:32.760 it has partial derivatives in it. 00:01:32.760 --> 00:01:35.880 And so we have a lot of equations like this. 00:01:35.880 --> 00:01:40.275 So for example, I work a lot on this equation. 00:02:01.302 --> 00:02:02.510 It's a poor choice, isn't it? 00:02:25.770 --> 00:02:30.120 OK, so this is the conservation equation for a chemical species 00:02:30.120 --> 00:02:32.710 in a reacting flow. 00:02:32.710 --> 00:02:40.870 So the change in the mass fraction of the species z, 00:02:40.870 --> 00:02:43.510 at a point, xyz at time, t. 00:02:43.510 --> 00:02:47.290 That differential is given by a diffusion term, a convection 00:02:47.290 --> 00:02:49.010 term, and a reaction term. 00:02:49.010 --> 00:02:50.990 So this is a pretty famous equation. 00:02:50.990 --> 00:02:53.950 This is part of the Navier-Stokes equations 00:02:53.950 --> 00:02:56.980 for the whole reactive flow system, which 00:02:56.980 --> 00:02:58.600 I hope you guys have seen already, 00:02:58.600 --> 00:03:01.784 and for sure you will see, if you haven't seen them yet. 00:03:01.784 --> 00:03:04.200 And there's a lot of other partial differential equations. 00:03:04.200 --> 00:03:05.740 I work a lot, also, on the Schrodinger equation. 00:03:05.740 --> 00:03:07.650 That's another partial differential equation. 00:03:07.650 --> 00:03:08.790 And what's special about the partial differential 00:03:08.790 --> 00:03:11.720 equations is that, in this case, this partial derivative 00:03:11.720 --> 00:03:15.820 is respect to time, holding all the spatial coordinates fixed. 00:03:15.820 --> 00:03:17.710 And these partial derivatives are respect 00:03:17.710 --> 00:03:22.280 to space, holding time fixed. 00:03:22.280 --> 00:03:24.800 A very important thing about partial derivatives, 00:03:24.800 --> 00:03:27.450 which you probably have encountered in 1040 already, 00:03:27.450 --> 00:03:29.650 is, it really depends what you hold fixed. 00:03:32.890 --> 00:03:34.420 And you get different results if you 00:03:34.420 --> 00:03:35.920 hold different things fixed. 00:03:35.920 --> 00:03:38.560 Now the implication, when you write an equation like this, 00:03:38.560 --> 00:03:42.090 is that whatever partials you see-- 00:03:42.090 --> 00:03:43.810 When I wrote this, if I have this 00:03:43.810 --> 00:03:45.670 and it has a term that's like d squared, 00:03:45.670 --> 00:03:49.400 dx squared, is one of the terms in this thing. 00:03:49.400 --> 00:03:50.340 This says t. 00:03:50.340 --> 00:03:51.200 This says x. 00:03:53.930 --> 00:03:57.020 The convention is, oh boy, I better hold x fixed. 00:03:57.020 --> 00:03:58.970 And here must be, I must hold t fixed. 00:03:58.970 --> 00:04:02.180 Because there's another one in the same equation. 00:04:02.180 --> 00:04:04.490 Now, you can change coordinates however you want. 00:04:04.490 --> 00:04:06.156 Just really watch out that you carefully 00:04:06.156 --> 00:04:08.660 follow the rules about what you hold 00:04:08.660 --> 00:04:10.040 fixed, when you change things. 00:04:10.040 --> 00:04:13.249 Otherwise you can cause all kinds of craziness. 00:04:13.249 --> 00:04:15.290 Now you guys probably took multivariable calculus 00:04:15.290 --> 00:04:17.440 at some point in your distant past. 00:04:17.440 --> 00:04:20.440 And they told you very carefully what the rules were. 00:04:20.440 --> 00:04:21.680 So follow them. 00:04:21.680 --> 00:04:24.235 And I suppose, when you're doing thermo right now? 00:04:24.235 --> 00:04:26.610 Have you started to do all this partial derivative stuff? 00:04:26.610 --> 00:04:28.485 Yeah, so you've seen how completely confusing 00:04:28.485 --> 00:04:30.984 it can be with negative signs showing up all over the place. 00:04:30.984 --> 00:04:32.330 And other things like this. 00:04:32.330 --> 00:04:34.230 I don't know if you've encountered this yet? 00:04:34.230 --> 00:04:37.890 At least, I thought it was confusing when I took it. 00:04:37.890 --> 00:04:39.620 So just be aware. 00:04:39.620 --> 00:04:42.640 It's the same thing. 00:04:42.640 --> 00:04:45.740 Now, often you do want to change the equations. 00:04:45.740 --> 00:04:47.740 So you can write down an equation like this. 00:04:47.740 --> 00:04:50.210 But, for example, if you have a special symmetry 00:04:50.210 --> 00:04:52.900 of the problem, like it's cylindrical, for example, 00:04:52.900 --> 00:04:55.391 then you might want to change the cylindrical coordinates. 00:04:55.391 --> 00:04:56.890 And then you have to just be careful 00:04:56.890 --> 00:04:59.150 that you write that correct Laplacian and cylindrical 00:04:59.150 --> 00:05:00.070 coordinates. 00:05:00.070 --> 00:05:03.310 Make sure it doesn't mess up anything about what 00:05:03.310 --> 00:05:05.350 you're holding fixed. 00:05:05.350 --> 00:05:08.550 Sometimes you might want to be crafty and use 00:05:08.550 --> 00:05:09.550 some weirdo coordinates. 00:05:09.550 --> 00:05:12.100 For example, if you wanted to solve the Schrodinger equation 00:05:12.100 --> 00:05:14.404 for h2 plus, it turns out there's 00:05:14.404 --> 00:05:15.820 a special coordinate system called 00:05:15.820 --> 00:05:17.480 the elliptic coordinate system. 00:05:17.480 --> 00:05:19.771 And in that coordinate system, the partial differential 00:05:19.771 --> 00:05:20.830 equation is separable. 00:05:20.830 --> 00:05:22.040 And so, if you do it in that coordinate system, 00:05:22.040 --> 00:05:23.831 you can actually solve analytical solution, 00:05:23.831 --> 00:05:26.860 which is pretty unusual for the Schroedinger equation. 00:05:26.860 --> 00:05:29.830 But it's a pretty goofball coordinate system 00:05:29.830 --> 00:05:32.770 and you really have to watch out, when you convert, 00:05:32.770 --> 00:05:34.990 to make sure you do it just right. 00:05:34.990 --> 00:05:37.390 But it's just the rules, chain rule and the rules 00:05:37.390 --> 00:05:40.032 about how to keep track of what's being held constant. 00:05:40.032 --> 00:05:41.740 And when you change, will tell them, how, 00:05:41.740 --> 00:05:43.323 being held constant, how to change it. 00:05:46.820 --> 00:05:50.710 Now in principle, this problem is 00:05:50.710 --> 00:05:52.870 like no different than the ODE BBP 00:05:52.870 --> 00:05:54.370 problems we were just doing. 00:05:54.370 --> 00:05:57.160 You just have a few more terms. 00:05:57.160 --> 00:05:59.740 And so, you can use finite element. 00:05:59.740 --> 00:06:02.200 You can use finite differences. 00:06:02.200 --> 00:06:05.230 You can use basis set methods. 00:06:05.230 --> 00:06:07.360 And it's really no different. 00:06:07.360 --> 00:06:11.410 So all those methods basically look the same. 00:06:11.410 --> 00:06:16.300 You just get equations that have more unknowns in them. 00:06:16.300 --> 00:06:19.720 Because, for example, if you're using finite element methods 00:06:19.720 --> 00:06:23.260 or basis set methods with a local basis set, 00:06:23.260 --> 00:06:27.310 you need to put points both in the x direction and the y 00:06:27.310 --> 00:06:29.320 direction, and in this case, the z direction, 00:06:29.320 --> 00:06:31.270 and then maybe also in the time direction. 00:06:31.270 --> 00:06:33.010 So now you have a lot of mesh points 00:06:33.010 --> 00:06:34.700 in all different directions. 00:06:34.700 --> 00:06:36.730 And so you get a lot of mesh points. 00:06:36.730 --> 00:06:39.437 And so fundamentally, there's no problem, 00:06:39.437 --> 00:06:41.020 but in practice, there's a big problem 00:06:41.020 --> 00:06:43.760 if the number of unknowns gets to be so large. 00:06:43.760 --> 00:06:45.170 That's the main problem. 00:06:47.901 --> 00:06:48.650 So think about it. 00:06:48.650 --> 00:06:50.620 You have some kind of equation like this. 00:06:50.620 --> 00:06:53.109 It's in three spatial dimensions and in time. 00:06:53.109 --> 00:06:54.650 How many points do you think you need 00:06:54.650 --> 00:07:00.150 to discretize your solution in the spatial coordinate? 00:07:00.150 --> 00:07:02.550 There might be somewhere, might be 100, I don't know, 00:07:02.550 --> 00:07:04.310 points might be enough. 00:07:04.310 --> 00:07:05.920 And so now I have x and y and z. 00:07:05.920 --> 00:07:11.242 So I'll need, well, just in COMSOL, you saw this on Monday. 00:07:11.242 --> 00:07:12.700 Suppose I just have two dimensions. 00:07:12.700 --> 00:07:16.790 Suppose I just have x and z. 00:07:16.790 --> 00:07:19.610 And I start putting points in. 00:07:19.610 --> 00:07:22.044 Well, I want to have some points along this way. 00:07:22.044 --> 00:07:23.960 And I want to have some points along this way. 00:07:26.820 --> 00:07:29.360 And so, I actually have an unknown here and an unknown 00:07:29.360 --> 00:07:31.270 here, and one there and one there, 00:07:31.270 --> 00:07:33.080 and one there and one there. 00:07:33.080 --> 00:07:34.480 There's quite a few of them. 00:07:34.480 --> 00:07:36.020 I think of how many points I have. 00:07:36.020 --> 00:07:41.470 Each of those is a point where I have a little basis function. 00:07:41.470 --> 00:07:45.179 And each one of those has a magnitude of how much weight 00:07:45.179 --> 00:07:46.720 I have on that basis function, things 00:07:46.720 --> 00:07:49.199 they called the d's before in the basis functions. 00:07:49.199 --> 00:07:50.240 And I have a lot of them. 00:07:50.240 --> 00:07:52.080 Now, it's going to be, if I have 100 points this way 00:07:52.080 --> 00:07:54.850 and 100 points this way, I have 10,000 of these little points. 00:07:54.850 --> 00:07:57.910 That means I have 10,000 unknowns I have to solve for. 00:07:57.910 --> 00:08:00.612 10,000 unknowns are a lot harder to solve for than 100 unknowns. 00:08:00.612 --> 00:08:03.070 and if you saw, on some of the problems we did with the ODE 00:08:03.070 --> 00:08:04.759 BBP problems, even with 100 unknowns, 00:08:04.759 --> 00:08:06.550 we were having a lot of trouble, sometimes, 00:08:06.550 --> 00:08:08.380 to get the real solution. 00:08:08.380 --> 00:08:09.520 All right? 00:08:09.520 --> 00:08:11.680 So you can just see how it just gets to be, like, 00:08:11.680 --> 00:08:12.850 a big numerical problem. 00:08:12.850 --> 00:08:13.891 So that's two dimensions. 00:08:13.891 --> 00:08:16.030 The third dimension, you'll get another 100 times 00:08:16.030 --> 00:08:16.750 as many points. 00:08:16.750 --> 00:08:18.860 You might have a million unknowns. 00:08:18.860 --> 00:08:20.610 And once you get up to a million unknowns, 00:08:20.610 --> 00:08:22.960 we're starting to talk serious math, now. 00:08:22.960 --> 00:08:24.670 OK? 00:08:24.670 --> 00:08:27.720 Even for your good computers that you have. 00:08:27.720 --> 00:08:31.470 In addition, it's not just you have the number of mesh points, 00:08:31.470 --> 00:08:34.600 but it's the number of mesh points 00:08:34.600 --> 00:08:39.390 times the number of variables at each mesh point. 00:08:39.390 --> 00:08:44.039 So in the Navier-Stokes equations, what variables 00:08:44.039 --> 00:08:44.858 do you have? 00:08:44.858 --> 00:08:49.250 AUDIENCE: [INAUDIBLE] 00:08:49.250 --> 00:08:50.130 PROFESSOR: Louder. 00:08:50.130 --> 00:08:51.050 AUDIENCE: [INAUDIBLE] 00:08:51.050 --> 00:08:52.070 PROFESSOR: Velocity and pressure, OK. 00:08:52.070 --> 00:08:53.330 So you have pressure. 00:08:53.330 --> 00:08:57.277 We have the different components of velocity. 00:08:57.277 --> 00:08:58.275 What else we got? 00:09:02.766 --> 00:09:04.263 AUDIENCE: [INAUDIBLE] 00:09:04.263 --> 00:09:06.259 PROFESSOR: Sorry, what? 00:09:06.259 --> 00:09:07.840 AUDIENCE: [INAUDIBLE] 00:09:07.840 --> 00:09:09.140 PROFESSOR: Yeah, though, yeah. 00:09:09.140 --> 00:09:11.340 Density or temperature. 00:09:11.340 --> 00:09:12.550 Yeah, temperature maybe, so. 00:09:12.550 --> 00:09:14.830 So maybe temperature would be what you might use. 00:09:14.830 --> 00:09:17.980 And then you'd have a state function maybe with these guys. 00:09:17.980 --> 00:09:20.826 What else would you have? 00:09:20.826 --> 00:09:22.700 So then you'd have all these species, right,? 00:09:22.700 --> 00:09:24.330 The mesh fractions are all the species. 00:09:24.330 --> 00:09:26.400 So you have, like, y, species one. 00:09:26.400 --> 00:09:27.925 Mesh fraction, species two. 00:09:27.925 --> 00:09:28.850 However many you got. 00:09:31.490 --> 00:09:33.680 So that's a pretty big set of variables 00:09:33.680 --> 00:09:37.280 and you need to know all these numbers at each point. 00:09:37.280 --> 00:09:38.560 So this might be-- 00:09:38.560 --> 00:09:39.560 how many have we've got? 00:09:39.560 --> 00:09:42.571 One, two, three, four, five, plus however many species 00:09:42.571 --> 00:09:43.070 we have? 00:09:43.070 --> 00:09:49.780 So this is n species plus 5 times the mesh. 00:09:49.780 --> 00:09:53.310 And I told you the mesh was what, a million? 00:09:53.310 --> 00:09:56.640 So maybe I have 10 million unknowns. 00:09:59.340 --> 00:10:05.930 So the issue, to a large extent, has to do with just, 00:10:05.930 --> 00:10:07.160 can I solve it? 00:10:07.160 --> 00:10:11.050 Can I solve a system of 10 million unknowns? 00:10:11.050 --> 00:10:13.360 10 million equations, 10 million unknowns. 00:10:13.360 --> 00:10:15.652 All right, so that's going to be one big issue. 00:10:15.652 --> 00:10:17.860 Now the other issues, we still have the same problems 00:10:17.860 --> 00:10:20.680 we had in ordinary differential equations, 00:10:20.680 --> 00:10:22.702 that we have to worry about numerical stability. 00:10:22.702 --> 00:10:24.910 We have to worry about actually, is the whole problem 00:10:24.910 --> 00:10:26.410 even stable to begin with? 00:10:26.410 --> 00:10:27.860 Physically? 00:10:27.860 --> 00:10:30.040 Because if it's not stable physically, then probably 00:10:30.040 --> 00:10:33.100 we're not going to get a stable solution. 00:10:33.100 --> 00:10:35.800 So you have all those kinds of problems 00:10:35.800 --> 00:10:37.780 that we had before in the ODE case, 00:10:37.780 --> 00:10:40.452 and now it's just sort of amplified by the fact, 00:10:40.452 --> 00:10:42.910 now we have boundary conditions in more than one dimension. 00:10:42.910 --> 00:10:45.311 So we have multiple chances to screw it up. 00:10:45.311 --> 00:10:47.560 And we're trying to figure out how to set the boundary 00:10:47.560 --> 00:10:49.240 conditions correctly. 00:10:49.240 --> 00:10:51.310 And we have to worry about stability in all 00:10:51.310 --> 00:10:54.030 the different directions. 00:10:54.030 --> 00:10:55.780 So those are the different kinds of things 00:10:55.780 --> 00:10:56.960 we have to worry about. 00:10:56.960 --> 00:11:00.008 But it's all the same problems you had before, just amplified. 00:11:03.840 --> 00:11:08.280 Now, because of this problem of the very large number 00:11:08.280 --> 00:11:12.900 of unknowns, currently, if you have 00:11:12.900 --> 00:11:16.680 a problem that has two dimensions, 00:11:16.680 --> 00:11:18.690 you can kind of very routinely solve it. 00:11:18.690 --> 00:11:20.880 And things like COMSOL on your laptop, 00:11:20.880 --> 00:11:22.410 you can usually do those. 00:11:22.410 --> 00:11:24.540 When you get up to three dimensions, 00:11:24.540 --> 00:11:28.290 you might be able to solve it or might not be able to solve it. 00:11:28.290 --> 00:11:30.930 And it could be, like, a really big challenge sometimes 00:11:30.930 --> 00:11:31.930 to solve it. 00:11:31.930 --> 00:11:34.430 And then in reality, we have problems like the Navier-Stokes 00:11:34.430 --> 00:11:35.513 is four dimensions, right? 00:11:35.513 --> 00:11:36.940 There's time, as well. 00:11:36.940 --> 00:11:38.646 So then that gets to be really hard. 00:11:38.646 --> 00:11:40.020 And then the Schrodinger equation 00:11:40.020 --> 00:11:43.590 has three times the number of electrons dimensions in it. 00:11:43.590 --> 00:11:46.470 So that gets impossibly hard. 00:11:46.470 --> 00:11:49.050 And so, specialized methods are developed 00:11:49.050 --> 00:11:52.105 to handle those systems with a lot of dimensionality. 00:11:52.105 --> 00:11:53.730 So, like, for the Schrodinger equation, 00:11:53.730 --> 00:11:56.100 people almost always do it with basis set methods. 00:11:56.100 --> 00:11:58.320 And there's been a huge effort, over decades, 00:11:58.320 --> 00:12:00.930 to figure out really clever basis sets that 00:12:00.930 --> 00:12:02.940 are really close to the real solutions, 00:12:02.940 --> 00:12:06.300 so that you don't need to use too many basis functions to get 00:12:06.300 --> 00:12:07.590 the good solution. 00:12:07.590 --> 00:12:09.881 So that's the way to keep the number of variables down. 00:12:12.054 --> 00:12:15.060 In the Navier-Stokes world, I'll show you 00:12:15.060 --> 00:12:17.010 what people do to try to deal with that. 00:12:19.650 --> 00:12:20.959 There's many other problems. 00:12:20.959 --> 00:12:22.750 I don't know all the problems in the world. 00:12:22.750 --> 00:12:26.490 But in general, when you have three or more dimensions 00:12:26.490 --> 00:12:30.004 in your PDE system, you're looking up special tricks. 00:12:30.004 --> 00:12:31.920 You're looking up, how do people in this field 00:12:31.920 --> 00:12:33.840 deal with this particular PDE? 00:12:33.840 --> 00:12:36.400 Special ways that are good for that kind of PDE. 00:12:39.265 --> 00:12:40.640 All right, so I said, one problem 00:12:40.640 --> 00:12:42.848 was that you could have, that the PDE system, itself, 00:12:42.848 --> 00:12:44.250 can be intrinsically unstable. 00:12:44.250 --> 00:12:48.440 So one example I know like that is, detonation problems. 00:12:48.440 --> 00:12:52.820 So if I have a system that has some explosive in it, 00:12:52.820 --> 00:12:55.130 or a mixture of hydrogen and oxygen even, 00:12:55.130 --> 00:12:59.720 and I have it in a tube or something, 00:12:59.720 --> 00:13:01.137 if I change my initial conditions, 00:13:01.137 --> 00:13:03.386 I can go from the case where it's stable, that it just 00:13:03.386 --> 00:13:05.100 sits there, a hydrogen and oxygen. 00:13:05.100 --> 00:13:08.239 Or it could turn into a regular burning. 00:13:08.239 --> 00:13:10.530 Or if I change the conditions a little bit differently, 00:13:10.530 --> 00:13:12.470 it could change into a detonation, where it actually 00:13:12.470 --> 00:13:14.100 sends a shockwave down the tube faster 00:13:14.100 --> 00:13:16.917 than the speed of sound, totally different situation. 00:13:16.917 --> 00:13:18.500 So those kind of systems can be really 00:13:18.500 --> 00:13:22.670 sensitive to the physical initial conditions. 00:13:22.670 --> 00:13:25.130 A very small spark can make a really big difference 00:13:25.130 --> 00:13:26.600 in a system like that. 00:13:26.600 --> 00:13:29.660 And also, in your numerical solution, 00:13:29.660 --> 00:13:31.820 if the thing is physically sensitive, 00:13:31.820 --> 00:13:34.134 it'll typically be sensitive also to numerical noise. 00:13:34.134 --> 00:13:36.050 If you get numerical noise at some point, that 00:13:36.050 --> 00:13:39.080 might be enough to kick it over from, say, not 00:13:39.080 --> 00:13:41.642 igniting at all to having a detonation, 00:13:41.642 --> 00:13:43.350 which is a completely different solution. 00:13:48.390 --> 00:13:51.270 But another case, like in biology, 00:13:51.270 --> 00:13:54.120 if you're trying to study the growth of a cell culture 00:13:54.120 --> 00:13:58.380 or the growth of a tumor, as you know, 00:13:58.380 --> 00:14:02.521 if somebody has a tumor and one of the cancer cells mutates 00:14:02.521 --> 00:14:04.020 and does something different, it can 00:14:04.020 --> 00:14:05.760 have a really different outcome of what 00:14:05.760 --> 00:14:07.920 happens to the patient, what happens to the tumor. 00:14:07.920 --> 00:14:10.210 The whole morphology can change completely. 00:14:10.210 --> 00:14:11.920 Everything about it can change. 00:14:11.920 --> 00:14:13.586 So that's the kind of system, also, that 00:14:13.586 --> 00:14:16.510 can be very sensitive to small details. 00:14:16.510 --> 00:14:18.450 A small change in one cell, for example, 00:14:18.450 --> 00:14:20.230 could be a really big difference. 00:14:20.230 --> 00:14:22.087 So this is not just in combustion problems 00:14:22.087 --> 00:14:23.670 you have these kinds of instabilities. 00:14:23.670 --> 00:14:26.087 All kinds of problems have this kind of thing. 00:14:26.087 --> 00:14:27.920 You can have just regular chemical reactors, 00:14:27.920 --> 00:14:30.650 like in multiple steady states, little tiny fluctuations 00:14:30.650 --> 00:14:32.960 can make it jump from one steady state to another one. 00:14:32.960 --> 00:14:34.751 That's another very famous kind of problem. 00:14:38.160 --> 00:14:40.640 There's another class of PDEs where 00:14:40.640 --> 00:14:45.830 it's stable in some sense, but it causes a problem. 00:14:45.830 --> 00:14:48.740 We call those hyperbolic. 00:14:48.740 --> 00:14:51.620 And those are like wave equations. 00:14:51.620 --> 00:14:53.120 The solutions are propagating waves. 00:14:53.120 --> 00:14:54.953 So the equations of acoustics, the equations 00:14:54.953 --> 00:14:57.560 of electromagnetism, are wave equations 00:14:57.560 --> 00:15:00.590 where, basically, a wave propagates more or less 00:15:00.590 --> 00:15:03.000 without any damping. 00:15:03.000 --> 00:15:05.540 So what happens, then, is you have some numerical noise that 00:15:05.540 --> 00:15:07.970 will make a little wavelet that will propagate, more or less, 00:15:07.970 --> 00:15:09.080 without any damping, too. 00:15:09.080 --> 00:15:11.660 And it'll keep bouncing back and forth inside your domain 00:15:11.660 --> 00:15:13.010 numerically. 00:15:13.010 --> 00:15:14.870 And if you keep introducing numerical noise, 00:15:14.870 --> 00:15:16.916 eventually the whole thing is just full of noise. 00:15:16.916 --> 00:15:18.290 And it doesn't have much relation 00:15:18.290 --> 00:15:20.240 to the real physical solution. 00:15:20.240 --> 00:15:23.810 So hyperbolic systems of differential equations 00:15:23.810 --> 00:15:25.490 need special solution methods. 00:15:25.490 --> 00:15:27.994 That in the numerical method, it's 00:15:27.994 --> 00:15:30.410 carefully trying to damp out that kind of propagating wave 00:15:30.410 --> 00:15:32.060 that's coming from the noise. 00:15:32.060 --> 00:15:33.680 And they set them up a special way. 00:15:33.680 --> 00:15:36.221 I'm not going to talk about how you solve them in this class, 00:15:36.221 --> 00:15:38.570 but there's a special, whole group of solution methods 00:15:38.570 --> 00:15:40.695 that people use who are solving acoustics equations 00:15:40.695 --> 00:15:44.540 and solving electromagnetic equations. 00:15:44.540 --> 00:15:47.219 If you get into solving shockwave equations-- 00:15:47.219 --> 00:15:49.010 I actually work for Professor [? Besant. ?] 00:15:49.010 --> 00:15:51.517 He has, like, electrochemical shocks in some of his systems. 00:15:51.517 --> 00:15:52.850 You have the same kind of thing. 00:15:52.850 --> 00:15:55.141 You have to use special mathematics to correctly handle 00:15:55.141 --> 00:15:58.166 that, special numerical tricks. 00:15:58.166 --> 00:16:00.290 But we're not going to get into that in this class. 00:16:00.290 --> 00:16:01.664 If you get into those things, you 00:16:01.664 --> 00:16:03.760 should take the PDE class for hyperbolic equations 00:16:03.760 --> 00:16:05.130 and they'll tell you just all about the solvers 00:16:05.130 --> 00:16:07.310 for those kinds of things and special tricks. 00:16:07.310 --> 00:16:08.976 And there's journal papers all about it. 00:16:11.352 --> 00:16:13.060 Here, what we're going to focus primarily 00:16:13.060 --> 00:16:14.800 on two kinds of problems: elliptic problems 00:16:14.800 --> 00:16:15.758 and parabolic problems. 00:16:18.640 --> 00:16:20.260 And those problems are ones where 00:16:20.260 --> 00:16:23.230 there's some dissipation that's very significant. 00:16:23.230 --> 00:16:26.290 And so, in an elliptic problem, you 00:16:26.290 --> 00:16:31.750 introduce some numerical noise and as you solve it, 00:16:31.750 --> 00:16:34.000 the numerical noise kind of goes away. 00:16:34.000 --> 00:16:36.524 Sort of like intrinsically stable all throughout. 00:16:36.524 --> 00:16:38.190 Now the real definition of these things, 00:16:38.190 --> 00:16:40.350 parabolic, hyperbolic, elliptic, has 00:16:40.350 --> 00:16:42.760 to do with flow of information. 00:16:42.760 --> 00:16:48.540 So in a hyperbolic system, information-- 00:16:48.540 --> 00:16:51.090 there's regions of the domain that you can make a change here 00:16:51.090 --> 00:16:54.460 and it doesn't make any effect on some domain part over there. 00:16:54.460 --> 00:16:56.250 And so that causes-- because that 00:16:56.250 --> 00:16:59.860 is the real situation, that also propagates 00:16:59.860 --> 00:17:02.330 into the numerical solution method. 00:17:02.330 --> 00:17:04.450 So for example, if I'm modeling a shockwave that's 00:17:04.450 --> 00:17:06.460 moving faster than the speed of sound, 00:17:06.460 --> 00:17:08.319 I can introduce any pressure fluctuation 00:17:08.319 --> 00:17:10.089 I want behind the shockwave, and it 00:17:10.089 --> 00:17:12.230 has no effect on the shockwave. 00:17:12.230 --> 00:17:14.920 Because the pressure waves from that pressure fluctuation 00:17:14.920 --> 00:17:15.760 are traveling at the speed of sound. 00:17:15.760 --> 00:17:17.410 They won't catch up to the shock. 00:17:17.410 --> 00:17:19.786 So they will have no effect whatsoever. 00:17:19.786 --> 00:17:21.369 So if I was judging what was happening 00:17:21.369 --> 00:17:23.202 by looking at what's happening in the shock, 00:17:23.202 --> 00:17:28.830 I can introduce any kind of random noise I want back over, 00:17:28.830 --> 00:17:30.260 downstream of the shock, I guess. 00:17:30.260 --> 00:17:31.050 Yeah? 00:17:31.050 --> 00:17:33.212 And it won't make any difference at all. 00:17:33.212 --> 00:17:34.836 And so, that's part of the reason why I 00:17:34.836 --> 00:17:37.110 need a special kind of solver. 00:17:37.110 --> 00:17:39.900 If I have elliptic equations, every point 00:17:39.900 --> 00:17:42.090 affects every other point. 00:17:42.090 --> 00:17:45.000 A famous case for that is, like, this steady state heat 00:17:45.000 --> 00:17:46.500 equation. 00:17:46.500 --> 00:17:50.310 You have some temperature source here, a heat source here, 00:17:50.310 --> 00:17:53.640 and you have some coolness source over here. 00:17:53.640 --> 00:17:57.585 And there's some heat flowing from the hot to the cold. 00:17:57.585 --> 00:18:00.104 And the heat flow is kind of slow enough 00:18:00.104 --> 00:18:02.520 that everything sort of affects everything else, since you 00:18:02.520 --> 00:18:04.890 actually have, the heat is trying to flow by, 00:18:04.890 --> 00:18:06.660 sort of, every path. 00:18:06.660 --> 00:18:10.282 And there's some insulation that's resisting its flow. 00:18:10.282 --> 00:18:11.740 And it has been in the connectivity 00:18:11.740 --> 00:18:12.869 of flow a certain way. 00:18:12.869 --> 00:18:14.910 And those kinds of equations are called elliptic. 00:18:14.910 --> 00:18:16.830 And the methods that we use for the relaxation mesh 00:18:16.830 --> 00:18:19.246 that we showed are really good for those kinds of methods, 00:18:19.246 --> 00:18:20.550 for those kind of problems. 00:18:20.550 --> 00:18:22.841 And then another kind of problem that we get into a lot 00:18:22.841 --> 00:18:23.890 is like this one. 00:18:23.890 --> 00:18:25.800 This is called parabolic. 00:18:25.800 --> 00:18:28.230 And typically, in the cases we see that are parabolic, 00:18:28.230 --> 00:18:32.390 we have time as part of the problem. 00:18:32.390 --> 00:18:34.374 Always what we think should happen 00:18:34.374 --> 00:18:35.790 is that what happens in the future 00:18:35.790 --> 00:18:40.170 depends on what happened in the past, but not vice versa. 00:18:40.170 --> 00:18:45.060 So you would expect that when you're 00:18:45.060 --> 00:18:48.662 computing what's happening at the next time point, 00:18:48.662 --> 00:18:50.870 it's going to depend on what happened at earlier time 00:18:50.870 --> 00:18:51.590 points. 00:18:51.590 --> 00:18:53.846 But if you try to do it the other way around, 00:18:53.846 --> 00:18:55.220 it would be kind of weird, right? 00:18:55.220 --> 00:18:57.011 You wouldn't expect that what in the future 00:18:57.011 --> 00:18:59.210 really affected stuff in the past. 00:18:59.210 --> 00:19:04.370 So it naturally leads us to sort of pose it as an ODE IVP, 00:19:04.370 --> 00:19:05.510 or IVP kind of problem. 00:19:05.510 --> 00:19:07.844 And we know initially, some time zero something, 00:19:07.844 --> 00:19:09.260 and we're trying to calculate what 00:19:09.260 --> 00:19:11.660 happens at some later time. 00:19:11.660 --> 00:19:14.210 And so those kinds of problems we call parabolic. 00:19:14.210 --> 00:19:15.700 And the methods that we're going to talk about mostly 00:19:15.700 --> 00:19:17.908 are ones for solving parabolic and elliptic problems. 00:19:21.120 --> 00:19:24.894 And you can read, in the middle of chapter 6, 00:19:24.894 --> 00:19:25.810 has a nice discussion. 00:19:25.810 --> 00:19:27.580 Just has the mathematical definitions 00:19:27.580 --> 00:19:28.621 of what these things are. 00:19:32.660 --> 00:19:36.710 Now, even in a steady state problem, 00:19:36.710 --> 00:19:39.860 the problem can have a sort of directionality to it. 00:19:39.860 --> 00:19:42.170 So if you have a problem that's very convective. 00:19:42.170 --> 00:19:43.820 So you have a flow in a pipe. 00:19:43.820 --> 00:19:46.362 And you have a nice, fast velocity flow down the pipe. 00:19:46.362 --> 00:19:48.320 What's happening downwind is very much affected 00:19:48.320 --> 00:19:51.750 by what happened upwind, but not vice versa. 00:19:51.750 --> 00:19:54.710 So maybe you're doing some reaction upwind. 00:19:54.710 --> 00:19:56.639 Say you're burning some hydrogen. 00:19:56.639 --> 00:19:58.430 You're going to get water vapor made upwind 00:19:58.430 --> 00:19:59.763 and it's going to flow downwind. 00:19:59.763 --> 00:20:03.440 And you're going to have steam downwind. 00:20:03.440 --> 00:20:06.117 If I squirt a little methane in downwind, 00:20:06.117 --> 00:20:08.450 it doesn't actually do anything to what happened upwind. 00:20:08.450 --> 00:20:10.491 I still have a hydrogen and oxygen flame up there 00:20:10.491 --> 00:20:12.260 and it's been burning and making steam. 00:20:12.260 --> 00:20:13.700 Do you see what I mean? 00:20:13.700 --> 00:20:15.470 So it has a directionality, even though I 00:20:15.470 --> 00:20:18.011 might do it as a steady state problem and not have time in it 00:20:18.011 --> 00:20:18.704 at all. 00:20:18.704 --> 00:20:20.870 And you've probably already done this, where you can 00:20:20.870 --> 00:20:23.300 convert, like, a flow reactor. 00:20:23.300 --> 00:20:26.464 You can write it, sort of, as a time thing or as a space thing. 00:20:26.464 --> 00:20:28.880 Time and space are so related by the velocity of the flow. 00:20:31.710 --> 00:20:34.050 We'll have similar kinds of phenomenon then. 00:20:34.050 --> 00:20:36.510 It's almost like a parabolic time system 00:20:36.510 --> 00:20:39.110 that, if the flow is fast enough, 00:20:39.110 --> 00:20:42.480 the fusion backup, anything moving upstream is negligible. 00:20:42.480 --> 00:20:45.640 It's all moving from upstream to downstream. 00:20:45.640 --> 00:20:49.220 And so then, you'll have the same kind of issues. 00:20:49.220 --> 00:20:51.060 In fact, you may even try to solve it 00:20:51.060 --> 00:20:54.750 by just starting it at the upstream end and competing 00:20:54.750 --> 00:20:57.872 stuff and then propagating what that output from that first one 00:20:57.872 --> 00:20:59.580 does to the next one downstream, and then 00:20:59.580 --> 00:21:00.538 do them, one at a time. 00:21:00.538 --> 00:21:02.130 That would be one possible way to try 00:21:02.130 --> 00:21:05.430 to solve those kinds of problems that are highly convective. 00:21:05.430 --> 00:21:08.572 As you slow the velocity down, then the diffusion 00:21:08.572 --> 00:21:10.530 will start to fight the velocity more and more. 00:21:10.530 --> 00:21:13.000 If you make the velocity very slow, 00:21:13.000 --> 00:21:15.030 the Peclet number very low, then you really 00:21:15.030 --> 00:21:17.520 have to worry about things diffusing back upstream. 00:21:17.520 --> 00:21:19.410 And then it's a different situation. 00:21:19.410 --> 00:21:23.260 And that would be more like the elliptic problems we have. 00:21:23.260 --> 00:21:25.740 Because things downstream would actually affect upstream, 00:21:25.740 --> 00:21:27.090 if the flow is really small. 00:21:27.090 --> 00:21:28.400 But typically for our problems, the velocity 00:21:28.400 --> 00:21:29.580 is always really high. 00:21:29.580 --> 00:21:32.434 Peclet numbers are really large and so not much goes 00:21:32.434 --> 00:21:33.600 from downstream to upstream. 00:21:33.600 --> 00:21:38.030 It's almost all upstream to downstream. 00:21:38.030 --> 00:21:41.235 This leads to a famous situation. 00:21:41.235 --> 00:21:50.627 If you look in the textbook at figures 6.7, 6.8, 00:21:50.627 --> 00:21:52.085 it's kind of a very famous problem. 00:21:54.730 --> 00:22:03.640 That if you just do this problem. 00:22:23.619 --> 00:22:24.910 So this is just an ODE problem. 00:22:29.760 --> 00:22:32.730 So I just have convection and diffusion, no reaction, 00:22:32.730 --> 00:22:33.900 very simple. 00:22:33.900 --> 00:22:35.070 This problem. 00:22:35.070 --> 00:22:38.100 You'd think that I could just use my finite differences, 00:22:38.100 --> 00:22:40.540 where I would put, I could use the center difference 00:22:40.540 --> 00:22:47.860 for here and here and change this into this obvious looking 00:22:47.860 --> 00:22:48.360 thing. 00:23:05.669 --> 00:23:08.210 You think that that would be a reasonable thing to do, right? 00:23:08.210 --> 00:23:09.380 And you can do the same one here. 00:23:09.380 --> 00:23:10.004 You can write-- 00:23:17.556 --> 00:23:18.930 I'm terrible with minuses, sorry. 00:23:27.620 --> 00:23:30.900 So these, the finite difference formulas. 00:23:30.900 --> 00:23:34.640 And so, you can write these like this. 00:23:34.640 --> 00:23:36.430 This stays equal to zero. 00:23:36.430 --> 00:23:39.640 And now it's just an algebraic problem to solve. 00:23:39.640 --> 00:23:41.530 A system of equations, you have equations 00:23:41.530 --> 00:23:43.430 like this for every value of n. 00:23:43.430 --> 00:23:46.290 And there's the mesh points. 00:23:46.290 --> 00:23:48.420 You guys have done this before, yes? 00:23:48.420 --> 00:23:49.170 Yes. 00:23:49.170 --> 00:23:51.870 OK, so famous problem is, if you actually 00:23:51.870 --> 00:23:57.740 try to solve it this way, unless you make delta z really tiny, 00:23:57.740 --> 00:23:58.880 you get crazy stuff. 00:23:58.880 --> 00:24:02.210 So on figure 6.7, they actually show 00:24:02.210 --> 00:24:08.780 what happens for different values of delta z, which 00:24:08.780 --> 00:24:11.930 is different values of what they call the local Peclet number. 00:24:11.930 --> 00:24:18.430 So there's a thing called local Peclet number, which is-- 00:24:23.680 --> 00:24:26.090 OK? 00:24:26.090 --> 00:24:29.120 If you make the local Peclet number too large, actually 00:24:29.120 --> 00:24:30.950 anywhere bigger than 2, and you try 00:24:30.950 --> 00:24:32.450 to solve this equation, what you get 00:24:32.450 --> 00:24:35.210 is oscillations, unphysical oscillations. 00:24:35.210 --> 00:24:37.800 They'll make the phi go negative. 00:24:37.800 --> 00:24:38.310 It's crazy. 00:24:41.107 --> 00:24:43.440 It looks like the simplest equation in the world, right? 00:24:43.440 --> 00:24:47.320 It's a linear equation, so what's the problem with this? 00:24:47.320 --> 00:24:49.250 But it doesn't work. 00:24:49.250 --> 00:24:53.590 And so, then you could think, well, why doesn't it work? 00:24:53.590 --> 00:24:55.660 And there's, like, multiple ways to explain this 00:24:55.660 --> 00:24:58.249 about why this doesn't work. 00:24:58.249 --> 00:25:00.040 I think that the best way to think about it 00:25:00.040 --> 00:25:01.400 is about the information flow. 00:25:01.400 --> 00:25:04.400 So if the local Peclet number is large, 00:25:04.400 --> 00:25:08.650 that means that phi is just flowing downwind 00:25:08.650 --> 00:25:10.660 and diffusion is not doing much. 00:25:10.660 --> 00:25:12.920 Because the Peclet number is big. 00:25:12.920 --> 00:25:18.140 And so, you really shouldn't use this formula 00:25:18.140 --> 00:25:19.970 to calculate the flow. 00:25:19.970 --> 00:25:24.769 Because what happens downwind, at point n plus one, I guess. 00:25:24.769 --> 00:25:26.060 Here we have the flow this way. 00:25:29.240 --> 00:25:32.520 So here's n minus one. 00:25:32.520 --> 00:25:34.770 Here's my point n, n plus 1. 00:25:38.470 --> 00:25:42.280 I really should not include, is this right, 00:25:42.280 --> 00:25:44.991 n plus 1 in this formula. 00:25:44.991 --> 00:25:46.240 Hope I have those signs right. 00:25:46.240 --> 00:25:49.266 I apologize if I had them backwards. 00:25:49.266 --> 00:25:50.890 Because I don't think that what happens 00:25:50.890 --> 00:25:54.320 downstream should really affect this at all. 00:25:54.320 --> 00:25:57.070 So really, I would do better to change from this center 00:25:57.070 --> 00:26:01.910 difference formula to what's called the upwind difference 00:26:01.910 --> 00:26:02.410 formula. 00:26:09.280 --> 00:26:10.930 Where they would say, OK, instead, 00:26:10.930 --> 00:26:23.810 let's use phi of n minus, the phi of n minus 1 over delta z. 00:26:23.810 --> 00:26:27.346 Just use that instead of using this formula. 00:26:27.346 --> 00:26:28.720 Now, you wouldn't think that that 00:26:28.720 --> 00:26:30.742 would make much difference. 00:26:30.742 --> 00:26:32.200 But it makes a gigantic difference. 00:26:32.200 --> 00:26:35.290 And so if you look at figure 6.8 in the textbook, 00:26:35.290 --> 00:26:40.180 you see you get perfectly stable solutions when you do that. 00:26:40.180 --> 00:26:41.740 But where you do this one, you get 00:26:41.740 --> 00:26:43.480 crazy, oscillatory, unphysical solutions, 00:26:43.480 --> 00:26:47.030 unless you choose delta z really tiny. 00:26:47.030 --> 00:26:50.487 Now, there's other ways to look at this. 00:26:50.487 --> 00:26:52.320 So I think this is the correct way, the best 00:26:52.320 --> 00:26:55.210 way to think about it is, it has to do with information flow. 00:26:55.210 --> 00:26:57.210 If you try to make your solution depend on stuff 00:26:57.210 --> 00:26:58.666 that it doesn't really depend on. 00:26:58.666 --> 00:27:00.540 So upwind does not really depend on downwind. 00:27:00.540 --> 00:27:03.360 If you force it to, by choosing to include 00:27:03.360 --> 00:27:07.050 that information here, you end up with crazy stuff. 00:27:07.050 --> 00:27:09.610 So that's one way to look at it. 00:27:09.610 --> 00:27:12.000 Another way to look at it is that, if you really 00:27:12.000 --> 00:27:14.260 want to compute these derivatives, 00:27:14.260 --> 00:27:15.960 you've got to keep delta z small. 00:27:15.960 --> 00:27:17.640 Because we're taking a differential. 00:27:17.640 --> 00:27:20.010 We're turning it into a finite difference. 00:27:20.010 --> 00:27:22.640 And if you take your delta z too big, 00:27:22.640 --> 00:27:27.680 then it's a terrible approximation to do this. 00:27:27.680 --> 00:27:28.639 And you can look at it. 00:27:28.639 --> 00:27:30.013 You can see that what we're doing 00:27:30.013 --> 00:27:32.219 is, when we have this problem where the local Peclet 00:27:32.219 --> 00:27:34.510 number is getting too large, is we're actually choosing 00:27:34.510 --> 00:27:35.885 delta z bigger than the, sort of, 00:27:35.885 --> 00:27:38.340 the thickness of the solution. 00:27:38.340 --> 00:27:41.270 So if you look at the solution, the analytical solution 00:27:41.270 --> 00:27:44.590 of this problem, it's like this. 00:27:52.190 --> 00:27:54.930 Where this has a very, very thin layer at the end. 00:27:54.930 --> 00:27:58.720 It's like the flow has pushed all your stuff to one side. 00:27:58.720 --> 00:28:03.520 And if you choose your delta z to be from here to here, 00:28:03.520 --> 00:28:05.530 that's how big your delta z is, and then you're 00:28:05.530 --> 00:28:11.680 computing the derivative of this point halfway along from here 00:28:11.680 --> 00:28:13.395 to where this is. 00:28:13.395 --> 00:28:15.180 You compute your derivative like that. 00:28:15.180 --> 00:28:18.520 It's not a very good, not a very accurate representation of what 00:28:18.520 --> 00:28:19.912 the real derivative is here. 00:28:19.912 --> 00:28:21.370 And in fact, you really should have 00:28:21.370 --> 00:28:23.411 a bunch of points here and here and here and here 00:28:23.411 --> 00:28:26.789 and here to really resolve the sharp edge. 00:28:26.789 --> 00:28:28.330 And so you did something really crazy 00:28:28.330 --> 00:28:31.260 to use a big value of delta z to start with. 00:28:31.260 --> 00:28:34.070 So that's another way to look at this problem. 00:28:34.070 --> 00:28:38.062 Now, the fix, by doing this upwind differencing, 00:28:38.062 --> 00:28:40.020 I think the best way to look at this is saying, 00:28:40.020 --> 00:28:42.020 well, I'm just going to make sure I only depend, 00:28:42.020 --> 00:28:45.069 I only make the equations depend on what's upwind, 00:28:45.069 --> 00:28:46.110 with the convective term. 00:28:46.110 --> 00:28:47.850 Because there's nothing convecting from downwind. 00:28:47.850 --> 00:28:48.930 So just leave that out. 00:28:48.930 --> 00:28:51.270 So that's a conceptual way to think about it. 00:28:51.270 --> 00:28:53.760 Another way to think about it, which is in the text, 00:28:53.760 --> 00:28:57.300 is that by doing this, you're introducing 00:28:57.300 --> 00:28:59.550 another thing called numerical diffusion, where 00:28:59.550 --> 00:29:03.120 you're actually, effectively increasing the diffusivity. 00:29:03.120 --> 00:29:05.317 Because this is a very poor approximation 00:29:05.317 --> 00:29:06.150 to the differential. 00:29:06.150 --> 00:29:10.300 This is a asymmetrical finite difference formula. 00:29:10.300 --> 00:29:13.260 It's not a very good value estimate of the derivative, 00:29:13.260 --> 00:29:15.930 if delta z is big. 00:29:15.930 --> 00:29:18.829 But you can write it out carefully and write this out 00:29:18.829 --> 00:29:20.370 and say, oh, this actually is sort of 00:29:20.370 --> 00:29:24.690 like this plus an effective diffusivity 00:29:24.690 --> 00:29:27.960 chosen very carefully, so the terms cancel out just right. 00:29:27.960 --> 00:29:29.490 And it turns out, if you look at it 00:29:29.490 --> 00:29:31.650 that way, the effective diffusion you've 00:29:31.650 --> 00:29:34.320 added, by using this instead of that, 00:29:34.320 --> 00:29:38.920 is just enough to make the local Peclet number, 00:29:38.920 --> 00:29:41.440 if you compute it with the effects of this d effective, 00:29:41.440 --> 00:29:42.940 to stay less than 2. 00:29:42.940 --> 00:29:45.130 Which is what you need as a condition to make 00:29:45.130 --> 00:29:46.130 this numerically stable. 00:29:49.040 --> 00:29:51.022 So I strongly suggest you read that. 00:29:51.022 --> 00:29:52.730 It's only two pages long in the textbook. 00:29:52.730 --> 00:29:53.180 It's very clear. 00:29:53.180 --> 00:29:55.638 Just read, just look at it if you haven't seen this before. 00:30:01.347 --> 00:30:03.930 There's a similar thing and it might be relevant to the COMSOL 00:30:03.930 --> 00:30:09.780 problem that Kristyn showed you on Monday, is you can-- 00:30:12.450 --> 00:30:17.430 In that case, there was a region of the boundary here 00:30:17.430 --> 00:30:20.650 that had some concentration c naught, 00:30:20.650 --> 00:30:24.790 and then there was, over here, the concentration was zero. 00:30:24.790 --> 00:30:25.890 The boundary is zero. 00:30:25.890 --> 00:30:31.420 Here, it's some number, 17, whatever number it was. 00:30:31.420 --> 00:30:39.470 And if you think of how to model what's going on, 00:30:39.470 --> 00:30:44.120 one way to look at it is, I have a diffusive flux coming in 00:30:44.120 --> 00:30:48.230 from the drug patch, diffusing the drug into the flow. 00:30:48.230 --> 00:30:52.870 And I have this big, giant flow, flowing along here. 00:30:52.870 --> 00:31:00.900 And the flow gets slower as I get closer to the wall because 00:31:00.900 --> 00:31:03.250 of the friction with the wall. 00:31:03.250 --> 00:31:06.290 And so if I looked somewhere right in here, 00:31:06.290 --> 00:31:08.170 I could just draw a little control volume, 00:31:08.170 --> 00:31:09.280 and look what's happening. 00:31:09.280 --> 00:31:11.262 I have a diffusive flux coming in here. 00:31:11.262 --> 00:31:12.970 I don't have any drug up here, so there's 00:31:12.970 --> 00:31:13.970 no drug coming in here. 00:31:13.970 --> 00:31:15.770 But I have drug leaving. 00:31:15.770 --> 00:31:19.950 So basically, it's coming up and it's making a right-hand turn. 00:31:19.950 --> 00:31:21.927 And so then I could think, well, if I 00:31:21.927 --> 00:31:23.510 was going to try to figure out, what's 00:31:23.510 --> 00:31:26.250 the steady state concentration of drug in there, I could say, 00:31:26.250 --> 00:31:27.120 well, there's a certain amount of drug 00:31:27.120 --> 00:31:28.744 entering from the diffusive flux, which 00:31:28.744 --> 00:31:30.885 would be just this times dcdy. 00:31:34.320 --> 00:31:35.741 This is the y direction. 00:31:38.970 --> 00:31:41.280 And then this part over here would 00:31:41.280 --> 00:31:46.720 be the velocity, actually, just times the concentration 00:31:46.720 --> 00:31:47.900 I have in here, right? 00:31:47.900 --> 00:31:52.160 So how much is flowing out. 00:31:52.160 --> 00:31:54.570 And it's probably, my units are screwed up so there's-- 00:31:54.570 --> 00:31:55.700 Well, maybe not. 00:31:55.700 --> 00:31:56.200 That's OK. 00:31:56.200 --> 00:31:57.610 Is that all right? 00:32:00.430 --> 00:32:01.725 Units are screwed up. 00:32:01.725 --> 00:32:02.850 You guys will get it right. 00:32:05.360 --> 00:32:06.946 It's probably an area. 00:32:06.946 --> 00:32:08.770 Something to do with this size right here. 00:32:08.770 --> 00:32:09.561 This one here, too. 00:32:14.640 --> 00:32:17.570 Is that right? 00:32:17.570 --> 00:32:20.747 So you could compute what the steady state 00:32:20.747 --> 00:32:23.330 concentration would be if it's just coming in and flowing out. 00:32:23.330 --> 00:32:25.320 And so, this is the finite volume view 00:32:25.320 --> 00:32:27.620 of this kind of problem. 00:32:27.620 --> 00:32:30.464 And so you can write the equations that way, instead. 00:32:30.464 --> 00:32:32.880 And what you really care about are what the velocity flows 00:32:32.880 --> 00:32:34.879 and the diffusive fluxes are on these boundaries 00:32:34.879 --> 00:32:36.040 around the point. 00:32:36.040 --> 00:32:38.550 So you don't try to compute things right at the point. 00:32:38.550 --> 00:32:39.940 You compute around it. 00:32:39.940 --> 00:32:42.200 And one way to look at this is, we just 00:32:42.200 --> 00:32:43.310 did finite differencing. 00:32:43.310 --> 00:32:45.750 We were looking at what was happening coming in. 00:32:49.390 --> 00:32:52.160 And then, this is actually a center difference 00:32:52.160 --> 00:32:56.360 around this interface, which is halfway between the two points. 00:32:56.360 --> 00:32:59.450 And so actually, this is a good approximation 00:32:59.450 --> 00:33:01.925 for the finite volume point of view. 00:33:01.925 --> 00:33:04.090 This is all right? 00:33:04.090 --> 00:33:05.810 Now, these are all different ways 00:33:05.810 --> 00:33:07.220 to write down the equations. 00:33:07.220 --> 00:33:10.100 All of them are correct in the limit 00:33:10.100 --> 00:33:11.694 that delta z goes to zero. 00:33:11.694 --> 00:33:13.610 If you have an infinite number of mesh points, 00:33:13.610 --> 00:33:15.350 all these are exactly the same. 00:33:15.350 --> 00:33:18.710 But they lead to really different numerical properties 00:33:18.710 --> 00:33:21.360 when delta z is large. 00:33:21.360 --> 00:33:23.910 And they're all inaccurate when delta z is large, 00:33:23.910 --> 00:33:26.390 so it's all about what happens with the inaccuracies. 00:33:30.910 --> 00:33:35.650 Now, we can go ahead and take any problem, even 00:33:35.650 --> 00:33:38.980 really complicated 3-D problems or 4-D problems like this, 00:33:38.980 --> 00:33:42.164 and we can write them with relaxation methods. 00:33:42.164 --> 00:33:43.955 And what we're basically doing in a problem 00:33:43.955 --> 00:33:48.430 like this is turning it into some function of y 00:33:48.430 --> 00:33:53.460 is equal to zero, where y is the value of all the unknowns. 00:33:57.320 --> 00:34:01.230 And there's a lot of equations. 00:34:01.230 --> 00:34:04.110 So I might have 100 million equations 00:34:04.110 --> 00:34:08.130 and 100 million unknowns that are the values, 00:34:08.130 --> 00:34:12.179 or every state variable at every state point at all times. 00:34:12.179 --> 00:34:14.389 And I can write it down, by finding differences, 00:34:14.389 --> 00:34:16.770 by finding elements, by finding volumes. 00:34:16.770 --> 00:34:20.239 Any way I want, I can write down an expression like this. 00:34:20.239 --> 00:34:21.929 And I know in the limit, as I make 00:34:21.929 --> 00:34:24.370 the number of elements of y go to infinity, 00:34:24.370 --> 00:34:27.567 it will actually be the real solution. 00:34:27.567 --> 00:34:29.400 And the problem is, my computer can't handle 00:34:29.400 --> 00:34:30.870 infinite number of unknowns. 00:34:30.870 --> 00:34:32.536 In fact, it's even going to have trouble 00:34:32.536 --> 00:34:33.719 with a 100 million unknowns. 00:34:33.719 --> 00:34:35.622 And so, I have to figure out what to do. 00:34:35.622 --> 00:34:37.080 So how would I solve this normally, 00:34:37.080 --> 00:34:39.570 is I would take the Jacobian of this and I would do, 00:34:39.570 --> 00:34:42.310 I'd take an initial guess and I'd update it and have 00:34:42.310 --> 00:34:48.889 the Jacobian times my change from the initial guess 00:34:48.889 --> 00:34:52.424 to the equal negative of f, right? 00:34:52.424 --> 00:34:54.090 You guys have done this a million times. 00:34:54.090 --> 00:34:57.030 As Newton-Raphson is how you'd find 00:34:57.030 --> 00:34:58.675 improved from an initial guess. 00:34:58.675 --> 00:35:00.300 You have to provide some initial guess. 00:35:00.300 --> 00:35:02.110 This is actually pretty hard if you have 100 million unknowns. 00:35:02.110 --> 00:35:05.000 You have to provide the value of the initial guess for every one 00:35:05.000 --> 00:35:06.300 of the 100 million. 00:35:06.300 --> 00:35:08.309 But somehow you did it. 00:35:08.309 --> 00:35:10.350 And now you want to refine that guess and this is 00:35:10.350 --> 00:35:12.960 the formula for it. 00:35:12.960 --> 00:35:15.000 And so, you just use backslash, right? 00:35:15.000 --> 00:35:19.320 But the problem here is that now f is a 100 million long vector 00:35:19.320 --> 00:35:21.970 and this is a 100 million squared matrix. 00:35:21.970 --> 00:35:23.820 And a 100 million squared is pretty big, 00:35:23.820 --> 00:35:26.050 10 to 16th elements in it. 00:35:26.050 --> 00:35:28.050 So I'm not going to just write that matrix 00:35:28.050 --> 00:35:30.060 down directly like this. 00:35:30.060 --> 00:35:33.630 Now, we cleverly chose local basis functions, 00:35:33.630 --> 00:35:37.480 so this is going to be sparse, this Jacobian. 00:35:37.480 --> 00:35:39.950 So most of the numbers, are those 10 to the 16th numbers 00:35:39.950 --> 00:35:41.942 in this matrix are zero. 00:35:41.942 --> 00:35:43.400 So we don't have to represent them. 00:35:43.400 --> 00:35:45.370 But there still might be quite a few. 00:35:45.370 --> 00:35:47.120 Probably the diagonal events are non-zero, 00:35:47.120 --> 00:35:49.530 so that's like 10 to the 8th of them right there. 00:35:49.530 --> 00:35:50.930 And we've got, probably, a couple of other bands. 00:35:50.930 --> 00:35:53.390 So I don't know, might be 10 to the 9th non-zero elements 00:35:53.390 --> 00:35:56.240 inside this thing, which is pretty nutty. 00:35:56.240 --> 00:35:58.790 And depending on how much memory they provided in the laptop 00:35:58.790 --> 00:36:01.370 they gave you, you might be filling up the memory 00:36:01.370 --> 00:36:02.961 already just to write down j. 00:36:02.961 --> 00:36:05.210 And certainly, if you tried to do Gaussian elimination 00:36:05.210 --> 00:36:08.477 to solve this, you're going to have a problem because you 00:36:08.477 --> 00:36:09.560 get what's called fill in. 00:36:09.560 --> 00:36:12.510 Do you guys remember this? 00:36:12.510 --> 00:36:14.900 Even if you start from a very sparse Jacobian, 00:36:14.900 --> 00:36:17.360 as you run the Gaussian elimination steps 00:36:17.360 --> 00:36:19.850 the number of non-zero entries in the intermediate matrices 00:36:19.850 --> 00:36:21.180 gets larger and larger. 00:36:21.180 --> 00:36:22.359 Remember this? 00:36:22.359 --> 00:36:24.900 And so, even if, initially, you can write down j and store it 00:36:24.900 --> 00:36:26.710 in your computer, you could easily 00:36:26.710 --> 00:36:31.330 overwhelm it with the third iteration or something. 00:36:31.330 --> 00:36:33.560 So you're not going to just use Gaussian elimination. 00:36:33.560 --> 00:36:34.580 So what do you have to do? 00:36:34.580 --> 00:36:35.663 You want to find a method. 00:36:35.663 --> 00:36:38.518 They're called direct. 00:36:38.518 --> 00:36:43.090 Direct methods to solve this kind of problem, 00:36:43.090 --> 00:36:44.730 in order to get your increment, which 00:36:44.730 --> 00:36:47.260 is going to be your update, to get a better estimate of what 00:36:47.260 --> 00:36:49.050 the solution is. 00:36:49.050 --> 00:36:51.870 And the direct methods are ones that you don't actually 00:36:51.870 --> 00:36:53.940 store the gigantic matrix. 00:36:53.940 --> 00:36:57.260 So you're trading off CPU time for memory. 00:36:57.260 --> 00:37:00.450 So you don't have enough memory to store the whole thing. 00:37:00.450 --> 00:37:03.810 Which we, normally you would do this with Gaussian elimination 00:37:03.810 --> 00:37:04.760 and huge steps. 00:37:04.760 --> 00:37:06.712 It would solve it perfectly, right? 00:37:06.712 --> 00:37:08.670 Right, Gaussian elimination's great, backslash. 00:37:08.670 --> 00:37:10.335 You guys have used it a few times. 00:37:10.335 --> 00:37:11.460 It's a really nice program. 00:37:11.460 --> 00:37:12.660 It always works. 00:37:12.660 --> 00:37:14.660 It's really nice. 00:37:14.660 --> 00:37:17.160 But you're going to give up that because you can't afford it 00:37:17.160 --> 00:37:19.260 because the memory, it's consuming too much memory 00:37:19.260 --> 00:37:21.093 and your cheapo department didn't buy enough 00:37:21.093 --> 00:37:23.940 RAM in your computer for you. 00:37:23.940 --> 00:37:27.540 So you're going to instead try a direct method which 00:37:27.540 --> 00:37:37.830 is, so it's trading off CPU time versus RAM. 00:37:37.830 --> 00:37:40.080 So you're going to reduce how much memory you actually 00:37:40.080 --> 00:37:42.754 actively have, but you're going to expect that it's 00:37:42.754 --> 00:37:43.920 going to cost more CPU time. 00:37:43.920 --> 00:37:45.590 Because otherwise, everybody would do something else besides 00:37:45.590 --> 00:37:48.360 Gaussian elimination all the time, and they don't. 00:37:48.360 --> 00:37:52.080 All right, so what's the method we're going to use? 00:37:52.080 --> 00:37:55.260 It turns out that variance on the conjugate gradient method 00:37:55.260 --> 00:37:56.310 turned out really well. 00:37:56.310 --> 00:38:00.510 So we mentioned this earlier, briefly. 00:38:00.510 --> 00:38:05.920 What you do is, instead of trying to solve this, 00:38:05.920 --> 00:38:08.140 you try to solve this. 00:38:26.060 --> 00:38:28.660 So you just try to minimize that, right? 00:38:28.660 --> 00:38:32.310 You know at the solution this is going to be zero. 00:38:32.310 --> 00:38:35.800 Right, jy plus [INAUDIBLE] zero when this is solved. 00:38:35.800 --> 00:38:38.740 So you try and vary and find the delta y 00:38:38.740 --> 00:38:42.390 that makes this go to zero by minimizing it. 00:38:42.390 --> 00:38:44.070 And then this method, it turns out 00:38:44.070 --> 00:38:47.640 that the conjugate gradient, if everything works great, 00:38:47.640 --> 00:38:50.220 this is guaranteed in n iterations 00:38:50.220 --> 00:38:53.820 to go directly to the solution of this kind of problem. 00:38:53.820 --> 00:38:56.320 Now n is large, now, because we have 10 to the 8th unknowns. 00:38:56.320 --> 00:38:58.190 So that's 10 to the 8th iterations. 00:38:58.190 --> 00:38:59.330 You do anything in 10 to the 8th iterations, 00:38:59.330 --> 00:39:01.300 you're not really sure it's really going to work. 00:39:01.300 --> 00:39:02.133 So that's a warning. 00:39:02.133 --> 00:39:05.950 But at least it should get, even after a few iterations, 00:39:05.950 --> 00:39:10.270 you should get a smaller value of this than you had before. 00:39:10.270 --> 00:39:11.380 And so this is the method. 00:39:11.380 --> 00:39:12.520 And what's really nice about this method, 00:39:12.520 --> 00:39:14.470 if you look into the details of how it works, 00:39:14.470 --> 00:39:17.740 it never has to store any intermediate matrices. 00:39:17.740 --> 00:39:20.650 So all it needs is the capability 00:39:20.650 --> 00:39:24.730 to multiply your matrix times some vector. 00:39:24.730 --> 00:39:29.970 And from that multiplication, it figures out a step direction. 00:39:29.970 --> 00:39:33.370 That's the trick of this method. 00:39:33.370 --> 00:39:36.422 And because it only has to do a forward multiplication, 00:39:36.422 --> 00:39:37.880 you don't have to actually store j. 00:39:37.880 --> 00:39:40.171 You can just compute elements of j, as needed, in order 00:39:40.171 --> 00:39:43.240 to evaluate this top product. 00:39:43.240 --> 00:39:46.730 This j times v. And you don't have to ever store 00:39:46.730 --> 00:39:48.390 the whole j at once. 00:39:48.390 --> 00:39:50.620 So it's really great for RAM to do this method. 00:39:53.330 --> 00:39:56.350 Now, as I said, when you do 10 to the 8th iterations, 00:39:56.350 --> 00:39:58.400 things don't always work out so well. 00:39:58.400 --> 00:40:00.220 And so people have found that when you go more than maybe 00:40:00.220 --> 00:40:01.844 10 or 20 iterations like this, you tend 00:40:01.844 --> 00:40:04.300 to pick up numerical problems. 00:40:04.300 --> 00:40:06.527 And so, they've worked out better methods. 00:40:06.527 --> 00:40:08.860 And there's one that people use a lot for these problems 00:40:08.860 --> 00:40:12.820 called bi cg stab. 00:40:12.820 --> 00:40:16.800 Which is biconjugate gradient stabilized. 00:40:16.800 --> 00:40:19.540 The applied mathematicians went crazy trying to figure out 00:40:19.540 --> 00:40:20.617 better methods like this. 00:40:20.617 --> 00:40:22.450 And this is the one that, I think currently, 00:40:22.450 --> 00:40:23.570 people think is the best. 00:40:23.570 --> 00:40:25.040 Though probably, there's a lot of research on this, 00:40:25.040 --> 00:40:26.623 so maybe there's even better ones now. 00:40:26.623 --> 00:40:29.870 But inside Matlab, I think this is the best on they have. 00:40:29.870 --> 00:40:34.770 And this is a method for doing this. 00:40:34.770 --> 00:40:35.830 Now, it's iterative. 00:40:35.830 --> 00:40:37.371 Now, let's remember what we're doing. 00:40:37.371 --> 00:40:38.910 We're trying to solve this problem. 00:40:38.910 --> 00:40:39.910 We started with a guess. 00:40:39.910 --> 00:40:42.535 We're going to break it up into an iterative problem like this, 00:40:42.535 --> 00:40:44.520 where we're going to do some steps, delta y. 00:40:44.520 --> 00:40:47.910 This is now doing iterative procedure 00:40:47.910 --> 00:40:49.710 in order to compute delta y. 00:40:49.710 --> 00:40:52.470 So we have an interim procedure and another iterative procedure 00:40:52.470 --> 00:40:54.499 inside the first iterative procedure. 00:40:54.499 --> 00:40:56.290 So this is going to be a lot more CPU time. 00:40:56.290 --> 00:40:57.940 So just to warn you. 00:41:01.024 --> 00:41:02.440 And it's also, it's not guaranteed 00:41:02.440 --> 00:41:04.356 like Gaussian elimination, to beautifully come 00:41:04.356 --> 00:41:06.097 to machine precision solution. 00:41:06.097 --> 00:41:07.430 It's going to get you something. 00:41:10.282 --> 00:41:11.740 But anyway, this is what people do. 00:41:11.740 --> 00:41:15.760 So that's one good thing to try. 00:41:15.760 --> 00:41:20.860 Now, this whole approach is sort of based on Newton-Raphson. 00:41:20.860 --> 00:41:24.310 We know that doesn't have the greatest radius of convergence. 00:41:24.310 --> 00:41:28.520 So you need to provide a pretty darn good initial guess. 00:41:28.520 --> 00:41:30.770 So how are we going to get a good initial guess if you 00:41:30.770 --> 00:41:32.210 have a 100 million unknowns? 00:41:35.250 --> 00:41:36.410 So what methods do we know? 00:41:36.410 --> 00:41:37.826 So one idea is you could do things 00:41:37.826 --> 00:41:40.210 like homotopy, stuff like that. 00:41:40.210 --> 00:41:42.830 If you could somehow take your PDE system, 00:41:42.830 --> 00:41:44.300 get rid of the nonlinear terms. 00:41:44.300 --> 00:41:47.741 You're really good at solving linear systems of PDEs. 00:41:47.741 --> 00:41:49.490 You start with that and then you gradually 00:41:49.490 --> 00:41:51.162 turn the non-linear term on. 00:41:51.162 --> 00:41:52.370 Maybe you could coax it over. 00:41:52.370 --> 00:41:56.180 So that's one possibility, if you're really good at PDEs. 00:41:56.180 --> 00:41:58.400 Another possibility that I've run 00:41:58.400 --> 00:42:03.140 into a lot is the problem you're trying to solve, for example, 00:42:03.140 --> 00:42:04.890 might be a steady state problem like this. 00:42:04.890 --> 00:42:07.040 Where this is zero and you try to figure out what 00:42:07.040 --> 00:42:11.900 the steady state situation is in a flow reactor, for example. 00:42:11.900 --> 00:42:14.240 And so, it's actually, your steady state problem 00:42:14.240 --> 00:42:17.010 came from a time dependent problem. 00:42:17.010 --> 00:42:19.640 And if you think that your time dependent problem really 00:42:19.640 --> 00:42:22.970 converges to the steady state problem, 00:42:22.970 --> 00:42:24.851 sort of without much concern about what 00:42:24.851 --> 00:42:26.350 the initial guess is, then you could 00:42:26.350 --> 00:42:28.550 put any random initial guess in and march it along 00:42:28.550 --> 00:42:32.300 in time enough, and eventually, it should end up 00:42:32.300 --> 00:42:35.420 at the solution you want. 00:42:35.420 --> 00:42:36.320 So that's the idea. 00:42:36.320 --> 00:42:37.903 So that's called a time marching idea. 00:42:37.903 --> 00:42:40.790 So there's one idea, sort of the homotopy continuation 00:42:40.790 --> 00:42:41.420 kind of idea. 00:42:41.420 --> 00:42:44.150 That's one idea about what to do for initial guesses. 00:42:44.150 --> 00:42:46.405 Another idea is the time marching idea. 00:42:46.405 --> 00:42:48.530 And the key about it is, you have to really believe 00:42:48.530 --> 00:42:50.900 that, physically, the system does march 00:42:50.900 --> 00:42:52.510 to the steady state you want. 00:42:52.510 --> 00:42:54.440 If it does, then this is a good idea. 00:42:54.440 --> 00:42:57.470 If it doesn't, then you're not going to end up where you want, 00:42:57.470 --> 00:42:59.990 because it's not going where you want. 00:42:59.990 --> 00:43:02.120 And so, you have to watch out. 00:43:02.120 --> 00:43:04.980 But in some situations, this is true, that things work. 00:43:04.980 --> 00:43:07.220 So I do a lot of combustion problems. 00:43:07.220 --> 00:43:11.390 A lot of these things, you light a spark on a stove burner. 00:43:11.390 --> 00:43:13.610 No matter how you light it, it ends up 00:43:13.610 --> 00:43:16.156 basically the same flame. 00:43:16.156 --> 00:43:17.030 So that's a good one. 00:43:17.030 --> 00:43:20.840 But I have some other ones where I light a candle 00:43:20.840 --> 00:43:23.840 and I'm in a convective field, a turbulent field, 00:43:23.840 --> 00:43:25.130 and it flickers. 00:43:25.130 --> 00:43:26.755 And no matter what, it always flickers. 00:43:26.755 --> 00:43:28.588 And I never get a steady state, so I'm never 00:43:28.588 --> 00:43:30.890 going to find a steady state solution with that one. 00:43:30.890 --> 00:43:34.717 And so I could time march until I retire. 00:43:34.717 --> 00:43:36.800 My computer's still burning CPU time all that time 00:43:36.800 --> 00:43:38.980 and it'll never get to a steady state solution. 00:43:38.980 --> 00:43:41.539 So you really have to know what the real situation is. 00:43:41.539 --> 00:43:43.580 But if you have one where it's going to converge, 00:43:43.580 --> 00:43:45.170 then this is a reasonable idea. 00:43:45.170 --> 00:43:48.900 Now let's talk about the time marching. 00:43:48.900 --> 00:43:53.436 The time marching, you have, say, dy, 00:43:53.436 --> 00:43:59.820 dt is equal to some f of y where this 00:43:59.820 --> 00:44:03.750 is the discretized spatial vergence. 00:44:03.750 --> 00:44:05.977 I took all the space dimensions and replaced them 00:44:05.977 --> 00:44:07.560 with the finite difference expressions 00:44:07.560 --> 00:44:08.820 for the derivatives. 00:44:08.820 --> 00:44:12.570 Or I did colocation or something to convert this 00:44:12.570 --> 00:44:15.912 into algebraic equations. 00:44:15.912 --> 00:44:18.120 So the right-hand side is now algebraic equations, no 00:44:18.120 --> 00:44:21.090 longer differentials. 00:44:21.090 --> 00:44:25.200 And now, I'm going to time march this and I just-- 00:44:25.200 --> 00:44:25.950 This is very long. 00:44:25.950 --> 00:44:28.470 It's a 100 million. 00:44:28.470 --> 00:44:30.720 So how are we going to time march that? 00:44:30.720 --> 00:44:33.690 Well there's kind of two schools of thought for this. 00:44:33.690 --> 00:44:38.090 One school of thought is, if I can use an explicit method, 00:44:38.090 --> 00:44:41.220 an explicit ODE solver, those are pretty nice. 00:44:41.220 --> 00:44:44.070 Because all I have to do is evaluate f and then 00:44:44.070 --> 00:44:47.610 I can march along and I never have to save any matrices. 00:44:47.610 --> 00:44:51.630 Because remember, look at those formulas, y, 00:44:51.630 --> 00:44:56.340 y new, by the explicit method, is 00:44:56.340 --> 00:45:00.320 equal to some function g of y old. 00:45:04.469 --> 00:45:05.760 And this is the update formula. 00:45:05.760 --> 00:45:06.760 It can be forward Euler. 00:45:06.760 --> 00:45:08.240 It could be RK 2. 00:45:08.240 --> 00:45:09.346 It could be whatever. 00:45:09.346 --> 00:45:10.970 But there's just some explicit formula. 00:45:10.970 --> 00:45:12.500 I plug in the old vector. 00:45:12.500 --> 00:45:13.730 I compute a new vector. 00:45:13.730 --> 00:45:18.050 I never had to do any inversions of matrices or anything. 00:45:18.050 --> 00:45:18.980 So that's good. 00:45:18.980 --> 00:45:21.530 And this is, indeed, the way that the Navier-Stokes 00:45:21.530 --> 00:45:25.310 equations are solved currently by the best possible solver. 00:45:25.310 --> 00:45:30.050 So the best solver in the world discretizes in space 00:45:30.050 --> 00:45:34.970 and then does ODE IVP using explicit method 00:45:34.970 --> 00:45:35.899 for marching it. 00:45:35.899 --> 00:45:37.440 They do problems that are pretty big. 00:45:37.440 --> 00:45:40.190 So they'll have like 100 million mesh points and then 100 00:45:40.190 --> 00:45:41.450 species. 00:45:41.450 --> 00:45:44.870 So they have 10 to the 10th unknowns. 00:45:44.870 --> 00:45:47.120 And at each time step, you're computing another vector 00:45:47.120 --> 00:45:49.310 that's 10 to the 10th long. 00:45:49.310 --> 00:45:51.270 So it's a pretty big problem. 00:45:51.270 --> 00:45:55.460 And they use, like, 40% of one of the national supercomputer 00:45:55.460 --> 00:45:58.310 centers will be running at one time for one job like this. 00:45:58.310 --> 00:45:59.701 But it works. 00:45:59.701 --> 00:46:01.700 You never have to store anything that's too big. 00:46:01.700 --> 00:46:04.070 The biggest thing you have to store are the vectors. 00:46:04.070 --> 00:46:05.030 Which, you're always going to have 00:46:05.030 --> 00:46:06.988 to store some object the size of your solution, 00:46:06.988 --> 00:46:09.230 right, if you're going to get your solution. 00:46:09.230 --> 00:46:10.880 So that's one possibility. 00:46:10.880 --> 00:46:16.240 Now, this is really good if an explicit solver can solve it. 00:46:16.240 --> 00:46:20.510 So that means it's just not stiff, your problem. 00:46:20.510 --> 00:46:24.140 And also, if you want a time accurate solution. 00:46:24.140 --> 00:46:26.810 So in that case, they actually want the time accurate 00:46:26.810 --> 00:46:27.405 solution. 00:46:27.405 --> 00:46:29.780 But suppose we're trying to really solve the steady state 00:46:29.780 --> 00:46:31.547 problem over here. 00:46:31.547 --> 00:46:33.630 Then we really don't care about the time accuracy. 00:46:33.630 --> 00:46:36.720 We're not going to report anything about our time steps. 00:46:36.720 --> 00:46:38.180 In the end, we're only going to report our steady state 00:46:38.180 --> 00:46:38.679 solution. 00:46:38.679 --> 00:46:40.044 That's all we cared about. 00:46:40.044 --> 00:46:41.460 So in that case, there's no reason 00:46:41.460 --> 00:46:44.580 to try to have a really high accuracy, explicit formula 00:46:44.580 --> 00:46:47.250 and try to make sure we all our time steps are exactly right. 00:46:47.250 --> 00:46:48.780 We're just going to throw away all those time points anyway. 00:46:48.780 --> 00:46:51.360 We'll just keep the final one as our initial guess 00:46:51.360 --> 00:46:54.630 for a Newton-Raphson solve to find the steady state. 00:46:54.630 --> 00:46:57.840 So in those cases, you might instead want to use-- 00:46:57.840 --> 00:46:59.850 So this is time accurate. 00:47:05.312 --> 00:47:06.270 There's another method. 00:47:06.270 --> 00:47:10.260 I don't know if I should call it time inaccurate. 00:47:10.260 --> 00:47:14.250 If you don't care, if the solution is really, 00:47:14.250 --> 00:47:16.320 what the solution yft is, all you 00:47:16.320 --> 00:47:19.770 care about is where you get to, then you can do other methods. 00:47:19.770 --> 00:47:22.230 And one of them is the backward Euler. 00:47:22.230 --> 00:47:30.380 So you can say, well, y new is equal to y old 00:47:30.380 --> 00:47:35.660 plus delta t times f of y new. 00:47:35.660 --> 00:47:39.150 Now, the advantage of this is, I can make delta t pretty large. 00:47:39.150 --> 00:47:41.250 This is guaranteed stable as long 00:47:41.250 --> 00:47:42.640 as the ODE system is stable. 00:47:42.640 --> 00:47:43.890 Remember? 00:47:43.890 --> 00:47:49.110 And so, this will go to this, go to the solution pretty well, 00:47:49.110 --> 00:47:50.745 to the steady state solution. 00:47:50.745 --> 00:47:52.620 It won't go there in the right amount of time 00:47:52.620 --> 00:47:54.760 because this is a very low order of formula. 00:47:54.760 --> 00:47:57.190 It's not really accurate. 00:47:57.190 --> 00:47:59.620 But it will end up to the real steady state solution. 00:47:59.620 --> 00:48:01.220 This is what's actually used a lot. 00:48:01.220 --> 00:48:04.210 However, the problem with this is, this is now 00:48:04.210 --> 00:48:07.530 an implicit equation. 00:48:07.530 --> 00:48:09.810 So we may have to solve it with some method 00:48:09.810 --> 00:48:11.904 like Newton-Raphson. 00:48:11.904 --> 00:48:13.570 But we got into this because we couldn't 00:48:13.570 --> 00:48:15.569 solve our Newton-Raphson steps because we didn't 00:48:15.569 --> 00:48:17.230 have a good initial guess. 00:48:17.230 --> 00:48:19.890 So then the question is, can we find a good initial guess 00:48:19.890 --> 00:48:21.839 for this problem? 00:48:21.839 --> 00:48:23.380 And in this case, you can because you 00:48:23.380 --> 00:48:27.970 can use the explicit formulas to again provide the initial guess 00:48:27.970 --> 00:48:30.110 for this implicit formula. 00:48:30.110 --> 00:48:33.220 So that's what they do, also. 00:48:33.220 --> 00:48:36.204 And so now, if we're doing an iterative procedure in order 00:48:36.204 --> 00:48:38.620 to get the initial guess for our other iterative procedure 00:48:38.620 --> 00:48:40.090 ever there, which we're going to break up 00:48:40.090 --> 00:48:42.715 into another iterative procedure inside it to solve every step. 00:48:45.910 --> 00:48:47.860 But this is currently what people do. 00:48:47.860 --> 00:48:49.360 So you guys are smart. 00:48:49.360 --> 00:48:51.610 Maybe you can figure a better way to do it. 00:48:51.610 --> 00:48:54.360 The world is waiting. 00:48:54.360 --> 00:48:55.360 I think I'll stop there. 00:48:57.880 --> 00:48:59.890 By the way, actually just names-- 00:48:59.890 --> 00:49:02.980 This kind of thing is called method of lines. 00:49:07.350 --> 00:49:09.630 And if so, if anybody ever says I solved something 00:49:09.630 --> 00:49:12.060 by method of lines, what they meant was, 00:49:12.060 --> 00:49:14.389 they discretized in some of the coordinates 00:49:14.389 --> 00:49:16.680 and they kept one of the coordinates as a differential. 00:49:16.680 --> 00:49:20.660 And then they solved it with an ODE solver at the end. 00:49:20.660 --> 00:49:23.630 And in certain systems like this, it's a smart thing to do. 00:49:23.630 --> 00:49:26.240 You need to be kind of lucky that the set up is right, 00:49:26.240 --> 00:49:27.900 so you can use it. 00:49:27.900 --> 00:49:29.960 But if you can, it can be pretty good. 00:49:29.960 --> 00:49:31.290 All right. 00:49:31.290 --> 00:49:33.190 See you on Friday.