1 00:00:01,540 --> 00:00:03,910 The following content is provided under a Creative 2 00:00:03,910 --> 00:00:05,300 Commons license. 3 00:00:05,300 --> 00:00:07,510 Your support will help MIT OpenCourseWare 4 00:00:07,510 --> 00:00:11,600 continue to offer high quality educational resources for free. 5 00:00:11,600 --> 00:00:14,140 To make a donation, or to view additional materials 6 00:00:14,140 --> 00:00:18,100 from hundreds of MIT courses, visit MIT OpenCourseWare 7 00:00:18,100 --> 00:00:18,980 at ocw.mit.edu. 8 00:00:23,125 --> 00:00:24,750 PROFESSOR: So today we're going to talk 9 00:00:24,750 --> 00:00:26,160 about partial differential equations. 10 00:00:26,160 --> 00:00:28,035 Let me just give a quick review where we are. 11 00:00:28,035 --> 00:00:31,640 So we're going to have a PDE lecture today. 12 00:00:31,640 --> 00:00:36,450 Then Professor Swan is going to do a review on Friday. 13 00:00:36,450 --> 00:00:40,980 There's no lecture on Monday, but you have the quiz, too, 14 00:00:40,980 --> 00:00:42,780 Monday night. 15 00:00:42,780 --> 00:00:46,452 Wednesday, next week is Veterans Day so it's a holiday. 16 00:00:46,452 --> 00:00:47,535 There's no classes at MIT. 17 00:00:50,042 --> 00:00:51,750 So the next time I'll be lecturing to you 18 00:00:51,750 --> 00:00:55,410 will be a week from Friday. 19 00:00:55,410 --> 00:00:58,340 The homework following the quiz will be posted 20 00:00:58,340 --> 00:01:01,620 and so you guys can get started on that. 21 00:01:01,620 --> 00:01:04,140 That homework involves a COMSOL problem 22 00:01:04,140 --> 00:01:06,480 which, for those of you who had trouble with COMSOL, 23 00:01:06,480 --> 00:01:07,680 might take some extra time, since it's 24 00:01:07,680 --> 00:01:08,971 the first time you're doing it. 25 00:01:08,971 --> 00:01:11,940 So I would suggest you at least give it a try early. 26 00:01:11,940 --> 00:01:14,010 I'm sure the TAs would be very happy to help you 27 00:01:14,010 --> 00:01:17,190 to force COMSOL to cooperate. 28 00:01:17,190 --> 00:01:18,540 Bend COMSOL to your will. 29 00:01:21,270 --> 00:01:24,990 All right, so we're going to talk about partial differential 30 00:01:24,990 --> 00:01:25,800 equations. 31 00:01:25,800 --> 00:01:28,002 And, as it sounds, what's different 32 00:01:28,002 --> 00:01:29,710 about a partial differential equation is, 33 00:01:29,710 --> 00:01:32,760 it has partial derivatives in it. 34 00:01:32,760 --> 00:01:35,880 And so we have a lot of equations like this. 35 00:01:35,880 --> 00:01:40,275 So for example, I work a lot on this equation. 36 00:02:01,302 --> 00:02:02,510 It's a poor choice, isn't it? 37 00:02:25,770 --> 00:02:30,120 OK, so this is the conservation equation for a chemical species 38 00:02:30,120 --> 00:02:32,710 in a reacting flow. 39 00:02:32,710 --> 00:02:40,870 So the change in the mass fraction of the species z, 40 00:02:40,870 --> 00:02:43,510 at a point, xyz at time, t. 41 00:02:43,510 --> 00:02:47,290 That differential is given by a diffusion term, a convection 42 00:02:47,290 --> 00:02:49,010 term, and a reaction term. 43 00:02:49,010 --> 00:02:50,990 So this is a pretty famous equation. 44 00:02:50,990 --> 00:02:53,950 This is part of the Navier-Stokes equations 45 00:02:53,950 --> 00:02:56,980 for the whole reactive flow system, which 46 00:02:56,980 --> 00:02:58,600 I hope you guys have seen already, 47 00:02:58,600 --> 00:03:01,784 and for sure you will see, if you haven't seen them yet. 48 00:03:01,784 --> 00:03:04,200 And there's a lot of other partial differential equations. 49 00:03:04,200 --> 00:03:05,740 I work a lot, also, on the Schrodinger equation. 50 00:03:05,740 --> 00:03:07,650 That's another partial differential equation. 51 00:03:07,650 --> 00:03:08,790 And what's special about the partial differential 52 00:03:08,790 --> 00:03:11,720 equations is that, in this case, this partial derivative 53 00:03:11,720 --> 00:03:15,820 is respect to time, holding all the spatial coordinates fixed. 54 00:03:15,820 --> 00:03:17,710 And these partial derivatives are respect 55 00:03:17,710 --> 00:03:22,280 to space, holding time fixed. 56 00:03:22,280 --> 00:03:24,800 A very important thing about partial derivatives, 57 00:03:24,800 --> 00:03:27,450 which you probably have encountered in 1040 already, 58 00:03:27,450 --> 00:03:29,650 is, it really depends what you hold fixed. 59 00:03:32,890 --> 00:03:34,420 And you get different results if you 60 00:03:34,420 --> 00:03:35,920 hold different things fixed. 61 00:03:35,920 --> 00:03:38,560 Now the implication, when you write an equation like this, 62 00:03:38,560 --> 00:03:42,090 is that whatever partials you see-- 63 00:03:42,090 --> 00:03:43,810 When I wrote this, if I have this 64 00:03:43,810 --> 00:03:45,670 and it has a term that's like d squared, 65 00:03:45,670 --> 00:03:49,400 dx squared, is one of the terms in this thing. 66 00:03:49,400 --> 00:03:50,340 This says t. 67 00:03:50,340 --> 00:03:51,200 This says x. 68 00:03:53,930 --> 00:03:57,020 The convention is, oh boy, I better hold x fixed. 69 00:03:57,020 --> 00:03:58,970 And here must be, I must hold t fixed. 70 00:03:58,970 --> 00:04:02,180 Because there's another one in the same equation. 71 00:04:02,180 --> 00:04:04,490 Now, you can change coordinates however you want. 72 00:04:04,490 --> 00:04:06,156 Just really watch out that you carefully 73 00:04:06,156 --> 00:04:08,660 follow the rules about what you hold 74 00:04:08,660 --> 00:04:10,040 fixed, when you change things. 75 00:04:10,040 --> 00:04:13,249 Otherwise you can cause all kinds of craziness. 76 00:04:13,249 --> 00:04:15,290 Now you guys probably took multivariable calculus 77 00:04:15,290 --> 00:04:17,440 at some point in your distant past. 78 00:04:17,440 --> 00:04:20,440 And they told you very carefully what the rules were. 79 00:04:20,440 --> 00:04:21,680 So follow them. 80 00:04:21,680 --> 00:04:24,235 And I suppose, when you're doing thermo right now? 81 00:04:24,235 --> 00:04:26,610 Have you started to do all this partial derivative stuff? 82 00:04:26,610 --> 00:04:28,485 Yeah, so you've seen how completely confusing 83 00:04:28,485 --> 00:04:30,984 it can be with negative signs showing up all over the place. 84 00:04:30,984 --> 00:04:32,330 And other things like this. 85 00:04:32,330 --> 00:04:34,230 I don't know if you've encountered this yet? 86 00:04:34,230 --> 00:04:37,890 At least, I thought it was confusing when I took it. 87 00:04:37,890 --> 00:04:39,620 So just be aware. 88 00:04:39,620 --> 00:04:42,640 It's the same thing. 89 00:04:42,640 --> 00:04:45,740 Now, often you do want to change the equations. 90 00:04:45,740 --> 00:04:47,740 So you can write down an equation like this. 91 00:04:47,740 --> 00:04:50,210 But, for example, if you have a special symmetry 92 00:04:50,210 --> 00:04:52,900 of the problem, like it's cylindrical, for example, 93 00:04:52,900 --> 00:04:55,391 then you might want to change the cylindrical coordinates. 94 00:04:55,391 --> 00:04:56,890 And then you have to just be careful 95 00:04:56,890 --> 00:04:59,150 that you write that correct Laplacian and cylindrical 96 00:04:59,150 --> 00:05:00,070 coordinates. 97 00:05:00,070 --> 00:05:03,310 Make sure it doesn't mess up anything about what 98 00:05:03,310 --> 00:05:05,350 you're holding fixed. 99 00:05:05,350 --> 00:05:08,550 Sometimes you might want to be crafty and use 100 00:05:08,550 --> 00:05:09,550 some weirdo coordinates. 101 00:05:09,550 --> 00:05:12,100 For example, if you wanted to solve the Schrodinger equation 102 00:05:12,100 --> 00:05:14,404 for h2 plus, it turns out there's 103 00:05:14,404 --> 00:05:15,820 a special coordinate system called 104 00:05:15,820 --> 00:05:17,480 the elliptic coordinate system. 105 00:05:17,480 --> 00:05:19,771 And in that coordinate system, the partial differential 106 00:05:19,771 --> 00:05:20,830 equation is separable. 107 00:05:20,830 --> 00:05:22,040 And so, if you do it in that coordinate system, 108 00:05:22,040 --> 00:05:23,831 you can actually solve analytical solution, 109 00:05:23,831 --> 00:05:26,860 which is pretty unusual for the Schroedinger equation. 110 00:05:26,860 --> 00:05:29,830 But it's a pretty goofball coordinate system 111 00:05:29,830 --> 00:05:32,770 and you really have to watch out, when you convert, 112 00:05:32,770 --> 00:05:34,990 to make sure you do it just right. 113 00:05:34,990 --> 00:05:37,390 But it's just the rules, chain rule and the rules 114 00:05:37,390 --> 00:05:40,032 about how to keep track of what's being held constant. 115 00:05:40,032 --> 00:05:41,740 And when you change, will tell them, how, 116 00:05:41,740 --> 00:05:43,323 being held constant, how to change it. 117 00:05:46,820 --> 00:05:50,710 Now in principle, this problem is 118 00:05:50,710 --> 00:05:52,870 like no different than the ODE BBP 119 00:05:52,870 --> 00:05:54,370 problems we were just doing. 120 00:05:54,370 --> 00:05:57,160 You just have a few more terms. 121 00:05:57,160 --> 00:05:59,740 And so, you can use finite element. 122 00:05:59,740 --> 00:06:02,200 You can use finite differences. 123 00:06:02,200 --> 00:06:05,230 You can use basis set methods. 124 00:06:05,230 --> 00:06:07,360 And it's really no different. 125 00:06:07,360 --> 00:06:11,410 So all those methods basically look the same. 126 00:06:11,410 --> 00:06:16,300 You just get equations that have more unknowns in them. 127 00:06:16,300 --> 00:06:19,720 Because, for example, if you're using finite element methods 128 00:06:19,720 --> 00:06:23,260 or basis set methods with a local basis set, 129 00:06:23,260 --> 00:06:27,310 you need to put points both in the x direction and the y 130 00:06:27,310 --> 00:06:29,320 direction, and in this case, the z direction, 131 00:06:29,320 --> 00:06:31,270 and then maybe also in the time direction. 132 00:06:31,270 --> 00:06:33,010 So now you have a lot of mesh points 133 00:06:33,010 --> 00:06:34,700 in all different directions. 134 00:06:34,700 --> 00:06:36,730 And so you get a lot of mesh points. 135 00:06:36,730 --> 00:06:39,437 And so fundamentally, there's no problem, 136 00:06:39,437 --> 00:06:41,020 but in practice, there's a big problem 137 00:06:41,020 --> 00:06:43,760 if the number of unknowns gets to be so large. 138 00:06:43,760 --> 00:06:45,170 That's the main problem. 139 00:06:47,901 --> 00:06:48,650 So think about it. 140 00:06:48,650 --> 00:06:50,620 You have some kind of equation like this. 141 00:06:50,620 --> 00:06:53,109 It's in three spatial dimensions and in time. 142 00:06:53,109 --> 00:06:54,650 How many points do you think you need 143 00:06:54,650 --> 00:07:00,150 to discretize your solution in the spatial coordinate? 144 00:07:00,150 --> 00:07:02,550 There might be somewhere, might be 100, I don't know, 145 00:07:02,550 --> 00:07:04,310 points might be enough. 146 00:07:04,310 --> 00:07:05,920 And so now I have x and y and z. 147 00:07:05,920 --> 00:07:11,242 So I'll need, well, just in COMSOL, you saw this on Monday. 148 00:07:11,242 --> 00:07:12,700 Suppose I just have two dimensions. 149 00:07:12,700 --> 00:07:16,790 Suppose I just have x and z. 150 00:07:16,790 --> 00:07:19,610 And I start putting points in. 151 00:07:19,610 --> 00:07:22,044 Well, I want to have some points along this way. 152 00:07:22,044 --> 00:07:23,960 And I want to have some points along this way. 153 00:07:26,820 --> 00:07:29,360 And so, I actually have an unknown here and an unknown 154 00:07:29,360 --> 00:07:31,270 here, and one there and one there, 155 00:07:31,270 --> 00:07:33,080 and one there and one there. 156 00:07:33,080 --> 00:07:34,480 There's quite a few of them. 157 00:07:34,480 --> 00:07:36,020 I think of how many points I have. 158 00:07:36,020 --> 00:07:41,470 Each of those is a point where I have a little basis function. 159 00:07:41,470 --> 00:07:45,179 And each one of those has a magnitude of how much weight 160 00:07:45,179 --> 00:07:46,720 I have on that basis function, things 161 00:07:46,720 --> 00:07:49,199 they called the d's before in the basis functions. 162 00:07:49,199 --> 00:07:50,240 And I have a lot of them. 163 00:07:50,240 --> 00:07:52,080 Now, it's going to be, if I have 100 points this way 164 00:07:52,080 --> 00:07:54,850 and 100 points this way, I have 10,000 of these little points. 165 00:07:54,850 --> 00:07:57,910 That means I have 10,000 unknowns I have to solve for. 166 00:07:57,910 --> 00:08:00,612 10,000 unknowns are a lot harder to solve for than 100 unknowns. 167 00:08:00,612 --> 00:08:03,070 and if you saw, on some of the problems we did with the ODE 168 00:08:03,070 --> 00:08:04,759 BBP problems, even with 100 unknowns, 169 00:08:04,759 --> 00:08:06,550 we were having a lot of trouble, sometimes, 170 00:08:06,550 --> 00:08:08,380 to get the real solution. 171 00:08:08,380 --> 00:08:09,520 All right? 172 00:08:09,520 --> 00:08:11,680 So you can just see how it just gets to be, like, 173 00:08:11,680 --> 00:08:12,850 a big numerical problem. 174 00:08:12,850 --> 00:08:13,891 So that's two dimensions. 175 00:08:13,891 --> 00:08:16,030 The third dimension, you'll get another 100 times 176 00:08:16,030 --> 00:08:16,750 as many points. 177 00:08:16,750 --> 00:08:18,860 You might have a million unknowns. 178 00:08:18,860 --> 00:08:20,610 And once you get up to a million unknowns, 179 00:08:20,610 --> 00:08:22,960 we're starting to talk serious math, now. 180 00:08:22,960 --> 00:08:24,670 OK? 181 00:08:24,670 --> 00:08:27,720 Even for your good computers that you have. 182 00:08:27,720 --> 00:08:31,470 In addition, it's not just you have the number of mesh points, 183 00:08:31,470 --> 00:08:34,600 but it's the number of mesh points 184 00:08:34,600 --> 00:08:39,390 times the number of variables at each mesh point. 185 00:08:39,390 --> 00:08:44,039 So in the Navier-Stokes equations, what variables 186 00:08:44,039 --> 00:08:44,858 do you have? 187 00:08:44,858 --> 00:08:49,250 AUDIENCE: [INAUDIBLE] 188 00:08:49,250 --> 00:08:50,130 PROFESSOR: Louder. 189 00:08:50,130 --> 00:08:51,050 AUDIENCE: [INAUDIBLE] 190 00:08:51,050 --> 00:08:52,070 PROFESSOR: Velocity and pressure, OK. 191 00:08:52,070 --> 00:08:53,330 So you have pressure. 192 00:08:53,330 --> 00:08:57,277 We have the different components of velocity. 193 00:08:57,277 --> 00:08:58,275 What else we got? 194 00:09:02,766 --> 00:09:04,263 AUDIENCE: [INAUDIBLE] 195 00:09:04,263 --> 00:09:06,259 PROFESSOR: Sorry, what? 196 00:09:06,259 --> 00:09:07,840 AUDIENCE: [INAUDIBLE] 197 00:09:07,840 --> 00:09:09,140 PROFESSOR: Yeah, though, yeah. 198 00:09:09,140 --> 00:09:11,340 Density or temperature. 199 00:09:11,340 --> 00:09:12,550 Yeah, temperature maybe, so. 200 00:09:12,550 --> 00:09:14,830 So maybe temperature would be what you might use. 201 00:09:14,830 --> 00:09:17,980 And then you'd have a state function maybe with these guys. 202 00:09:17,980 --> 00:09:20,826 What else would you have? 203 00:09:20,826 --> 00:09:22,700 So then you'd have all these species, right,? 204 00:09:22,700 --> 00:09:24,330 The mesh fractions are all the species. 205 00:09:24,330 --> 00:09:26,400 So you have, like, y, species one. 206 00:09:26,400 --> 00:09:27,925 Mesh fraction, species two. 207 00:09:27,925 --> 00:09:28,850 However many you got. 208 00:09:31,490 --> 00:09:33,680 So that's a pretty big set of variables 209 00:09:33,680 --> 00:09:37,280 and you need to know all these numbers at each point. 210 00:09:37,280 --> 00:09:38,560 So this might be-- 211 00:09:38,560 --> 00:09:39,560 how many have we've got? 212 00:09:39,560 --> 00:09:42,571 One, two, three, four, five, plus however many species 213 00:09:42,571 --> 00:09:43,070 we have? 214 00:09:43,070 --> 00:09:49,780 So this is n species plus 5 times the mesh. 215 00:09:49,780 --> 00:09:53,310 And I told you the mesh was what, a million? 216 00:09:53,310 --> 00:09:56,640 So maybe I have 10 million unknowns. 217 00:09:59,340 --> 00:10:05,930 So the issue, to a large extent, has to do with just, 218 00:10:05,930 --> 00:10:07,160 can I solve it? 219 00:10:07,160 --> 00:10:11,050 Can I solve a system of 10 million unknowns? 220 00:10:11,050 --> 00:10:13,360 10 million equations, 10 million unknowns. 221 00:10:13,360 --> 00:10:15,652 All right, so that's going to be one big issue. 222 00:10:15,652 --> 00:10:17,860 Now the other issues, we still have the same problems 223 00:10:17,860 --> 00:10:20,680 we had in ordinary differential equations, 224 00:10:20,680 --> 00:10:22,702 that we have to worry about numerical stability. 225 00:10:22,702 --> 00:10:24,910 We have to worry about actually, is the whole problem 226 00:10:24,910 --> 00:10:26,410 even stable to begin with? 227 00:10:26,410 --> 00:10:27,860 Physically? 228 00:10:27,860 --> 00:10:30,040 Because if it's not stable physically, then probably 229 00:10:30,040 --> 00:10:33,100 we're not going to get a stable solution. 230 00:10:33,100 --> 00:10:35,800 So you have all those kinds of problems 231 00:10:35,800 --> 00:10:37,780 that we had before in the ODE case, 232 00:10:37,780 --> 00:10:40,452 and now it's just sort of amplified by the fact, 233 00:10:40,452 --> 00:10:42,910 now we have boundary conditions in more than one dimension. 234 00:10:42,910 --> 00:10:45,311 So we have multiple chances to screw it up. 235 00:10:45,311 --> 00:10:47,560 And we're trying to figure out how to set the boundary 236 00:10:47,560 --> 00:10:49,240 conditions correctly. 237 00:10:49,240 --> 00:10:51,310 And we have to worry about stability in all 238 00:10:51,310 --> 00:10:54,030 the different directions. 239 00:10:54,030 --> 00:10:55,780 So those are the different kinds of things 240 00:10:55,780 --> 00:10:56,960 we have to worry about. 241 00:10:56,960 --> 00:11:00,008 But it's all the same problems you had before, just amplified. 242 00:11:03,840 --> 00:11:08,280 Now, because of this problem of the very large number 243 00:11:08,280 --> 00:11:12,900 of unknowns, currently, if you have 244 00:11:12,900 --> 00:11:16,680 a problem that has two dimensions, 245 00:11:16,680 --> 00:11:18,690 you can kind of very routinely solve it. 246 00:11:18,690 --> 00:11:20,880 And things like COMSOL on your laptop, 247 00:11:20,880 --> 00:11:22,410 you can usually do those. 248 00:11:22,410 --> 00:11:24,540 When you get up to three dimensions, 249 00:11:24,540 --> 00:11:28,290 you might be able to solve it or might not be able to solve it. 250 00:11:28,290 --> 00:11:30,930 And it could be, like, a really big challenge sometimes 251 00:11:30,930 --> 00:11:31,930 to solve it. 252 00:11:31,930 --> 00:11:34,430 And then in reality, we have problems like the Navier-Stokes 253 00:11:34,430 --> 00:11:35,513 is four dimensions, right? 254 00:11:35,513 --> 00:11:36,940 There's time, as well. 255 00:11:36,940 --> 00:11:38,646 So then that gets to be really hard. 256 00:11:38,646 --> 00:11:40,020 And then the Schrodinger equation 257 00:11:40,020 --> 00:11:43,590 has three times the number of electrons dimensions in it. 258 00:11:43,590 --> 00:11:46,470 So that gets impossibly hard. 259 00:11:46,470 --> 00:11:49,050 And so, specialized methods are developed 260 00:11:49,050 --> 00:11:52,105 to handle those systems with a lot of dimensionality. 261 00:11:52,105 --> 00:11:53,730 So, like, for the Schrodinger equation, 262 00:11:53,730 --> 00:11:56,100 people almost always do it with basis set methods. 263 00:11:56,100 --> 00:11:58,320 And there's been a huge effort, over decades, 264 00:11:58,320 --> 00:12:00,930 to figure out really clever basis sets that 265 00:12:00,930 --> 00:12:02,940 are really close to the real solutions, 266 00:12:02,940 --> 00:12:06,300 so that you don't need to use too many basis functions to get 267 00:12:06,300 --> 00:12:07,590 the good solution. 268 00:12:07,590 --> 00:12:09,881 So that's the way to keep the number of variables down. 269 00:12:12,054 --> 00:12:15,060 In the Navier-Stokes world, I'll show you 270 00:12:15,060 --> 00:12:17,010 what people do to try to deal with that. 271 00:12:19,650 --> 00:12:20,959 There's many other problems. 272 00:12:20,959 --> 00:12:22,750 I don't know all the problems in the world. 273 00:12:22,750 --> 00:12:26,490 But in general, when you have three or more dimensions 274 00:12:26,490 --> 00:12:30,004 in your PDE system, you're looking up special tricks. 275 00:12:30,004 --> 00:12:31,920 You're looking up, how do people in this field 276 00:12:31,920 --> 00:12:33,840 deal with this particular PDE? 277 00:12:33,840 --> 00:12:36,400 Special ways that are good for that kind of PDE. 278 00:12:39,265 --> 00:12:40,640 All right, so I said, one problem 279 00:12:40,640 --> 00:12:42,848 was that you could have, that the PDE system, itself, 280 00:12:42,848 --> 00:12:44,250 can be intrinsically unstable. 281 00:12:44,250 --> 00:12:48,440 So one example I know like that is, detonation problems. 282 00:12:48,440 --> 00:12:52,820 So if I have a system that has some explosive in it, 283 00:12:52,820 --> 00:12:55,130 or a mixture of hydrogen and oxygen even, 284 00:12:55,130 --> 00:12:59,720 and I have it in a tube or something, 285 00:12:59,720 --> 00:13:01,137 if I change my initial conditions, 286 00:13:01,137 --> 00:13:03,386 I can go from the case where it's stable, that it just 287 00:13:03,386 --> 00:13:05,100 sits there, a hydrogen and oxygen. 288 00:13:05,100 --> 00:13:08,239 Or it could turn into a regular burning. 289 00:13:08,239 --> 00:13:10,530 Or if I change the conditions a little bit differently, 290 00:13:10,530 --> 00:13:12,470 it could change into a detonation, where it actually 291 00:13:12,470 --> 00:13:14,100 sends a shockwave down the tube faster 292 00:13:14,100 --> 00:13:16,917 than the speed of sound, totally different situation. 293 00:13:16,917 --> 00:13:18,500 So those kind of systems can be really 294 00:13:18,500 --> 00:13:22,670 sensitive to the physical initial conditions. 295 00:13:22,670 --> 00:13:25,130 A very small spark can make a really big difference 296 00:13:25,130 --> 00:13:26,600 in a system like that. 297 00:13:26,600 --> 00:13:29,660 And also, in your numerical solution, 298 00:13:29,660 --> 00:13:31,820 if the thing is physically sensitive, 299 00:13:31,820 --> 00:13:34,134 it'll typically be sensitive also to numerical noise. 300 00:13:34,134 --> 00:13:36,050 If you get numerical noise at some point, that 301 00:13:36,050 --> 00:13:39,080 might be enough to kick it over from, say, not 302 00:13:39,080 --> 00:13:41,642 igniting at all to having a detonation, 303 00:13:41,642 --> 00:13:43,350 which is a completely different solution. 304 00:13:48,390 --> 00:13:51,270 But another case, like in biology, 305 00:13:51,270 --> 00:13:54,120 if you're trying to study the growth of a cell culture 306 00:13:54,120 --> 00:13:58,380 or the growth of a tumor, as you know, 307 00:13:58,380 --> 00:14:02,521 if somebody has a tumor and one of the cancer cells mutates 308 00:14:02,521 --> 00:14:04,020 and does something different, it can 309 00:14:04,020 --> 00:14:05,760 have a really different outcome of what 310 00:14:05,760 --> 00:14:07,920 happens to the patient, what happens to the tumor. 311 00:14:07,920 --> 00:14:10,210 The whole morphology can change completely. 312 00:14:10,210 --> 00:14:11,920 Everything about it can change. 313 00:14:11,920 --> 00:14:13,586 So that's the kind of system, also, that 314 00:14:13,586 --> 00:14:16,510 can be very sensitive to small details. 315 00:14:16,510 --> 00:14:18,450 A small change in one cell, for example, 316 00:14:18,450 --> 00:14:20,230 could be a really big difference. 317 00:14:20,230 --> 00:14:22,087 So this is not just in combustion problems 318 00:14:22,087 --> 00:14:23,670 you have these kinds of instabilities. 319 00:14:23,670 --> 00:14:26,087 All kinds of problems have this kind of thing. 320 00:14:26,087 --> 00:14:27,920 You can have just regular chemical reactors, 321 00:14:27,920 --> 00:14:30,650 like in multiple steady states, little tiny fluctuations 322 00:14:30,650 --> 00:14:32,960 can make it jump from one steady state to another one. 323 00:14:32,960 --> 00:14:34,751 That's another very famous kind of problem. 324 00:14:38,160 --> 00:14:40,640 There's another class of PDEs where 325 00:14:40,640 --> 00:14:45,830 it's stable in some sense, but it causes a problem. 326 00:14:45,830 --> 00:14:48,740 We call those hyperbolic. 327 00:14:48,740 --> 00:14:51,620 And those are like wave equations. 328 00:14:51,620 --> 00:14:53,120 The solutions are propagating waves. 329 00:14:53,120 --> 00:14:54,953 So the equations of acoustics, the equations 330 00:14:54,953 --> 00:14:57,560 of electromagnetism, are wave equations 331 00:14:57,560 --> 00:15:00,590 where, basically, a wave propagates more or less 332 00:15:00,590 --> 00:15:03,000 without any damping. 333 00:15:03,000 --> 00:15:05,540 So what happens, then, is you have some numerical noise that 334 00:15:05,540 --> 00:15:07,970 will make a little wavelet that will propagate, more or less, 335 00:15:07,970 --> 00:15:09,080 without any damping, too. 336 00:15:09,080 --> 00:15:11,660 And it'll keep bouncing back and forth inside your domain 337 00:15:11,660 --> 00:15:13,010 numerically. 338 00:15:13,010 --> 00:15:14,870 And if you keep introducing numerical noise, 339 00:15:14,870 --> 00:15:16,916 eventually the whole thing is just full of noise. 340 00:15:16,916 --> 00:15:18,290 And it doesn't have much relation 341 00:15:18,290 --> 00:15:20,240 to the real physical solution. 342 00:15:20,240 --> 00:15:23,810 So hyperbolic systems of differential equations 343 00:15:23,810 --> 00:15:25,490 need special solution methods. 344 00:15:25,490 --> 00:15:27,994 That in the numerical method, it's 345 00:15:27,994 --> 00:15:30,410 carefully trying to damp out that kind of propagating wave 346 00:15:30,410 --> 00:15:32,060 that's coming from the noise. 347 00:15:32,060 --> 00:15:33,680 And they set them up a special way. 348 00:15:33,680 --> 00:15:36,221 I'm not going to talk about how you solve them in this class, 349 00:15:36,221 --> 00:15:38,570 but there's a special, whole group of solution methods 350 00:15:38,570 --> 00:15:40,695 that people use who are solving acoustics equations 351 00:15:40,695 --> 00:15:44,540 and solving electromagnetic equations. 352 00:15:44,540 --> 00:15:47,219 If you get into solving shockwave equations-- 353 00:15:47,219 --> 00:15:49,010 I actually work for Professor [? Besant. ?] 354 00:15:49,010 --> 00:15:51,517 He has, like, electrochemical shocks in some of his systems. 355 00:15:51,517 --> 00:15:52,850 You have the same kind of thing. 356 00:15:52,850 --> 00:15:55,141 You have to use special mathematics to correctly handle 357 00:15:55,141 --> 00:15:58,166 that, special numerical tricks. 358 00:15:58,166 --> 00:16:00,290 But we're not going to get into that in this class. 359 00:16:00,290 --> 00:16:01,664 If you get into those things, you 360 00:16:01,664 --> 00:16:03,760 should take the PDE class for hyperbolic equations 361 00:16:03,760 --> 00:16:05,130 and they'll tell you just all about the solvers 362 00:16:05,130 --> 00:16:07,310 for those kinds of things and special tricks. 363 00:16:07,310 --> 00:16:08,976 And there's journal papers all about it. 364 00:16:11,352 --> 00:16:13,060 Here, what we're going to focus primarily 365 00:16:13,060 --> 00:16:14,800 on two kinds of problems: elliptic problems 366 00:16:14,800 --> 00:16:15,758 and parabolic problems. 367 00:16:18,640 --> 00:16:20,260 And those problems are ones where 368 00:16:20,260 --> 00:16:23,230 there's some dissipation that's very significant. 369 00:16:23,230 --> 00:16:26,290 And so, in an elliptic problem, you 370 00:16:26,290 --> 00:16:31,750 introduce some numerical noise and as you solve it, 371 00:16:31,750 --> 00:16:34,000 the numerical noise kind of goes away. 372 00:16:34,000 --> 00:16:36,524 Sort of like intrinsically stable all throughout. 373 00:16:36,524 --> 00:16:38,190 Now the real definition of these things, 374 00:16:38,190 --> 00:16:40,350 parabolic, hyperbolic, elliptic, has 375 00:16:40,350 --> 00:16:42,760 to do with flow of information. 376 00:16:42,760 --> 00:16:48,540 So in a hyperbolic system, information-- 377 00:16:48,540 --> 00:16:51,090 there's regions of the domain that you can make a change here 378 00:16:51,090 --> 00:16:54,460 and it doesn't make any effect on some domain part over there. 379 00:16:54,460 --> 00:16:56,250 And so that causes-- because that 380 00:16:56,250 --> 00:16:59,860 is the real situation, that also propagates 381 00:16:59,860 --> 00:17:02,330 into the numerical solution method. 382 00:17:02,330 --> 00:17:04,450 So for example, if I'm modeling a shockwave that's 383 00:17:04,450 --> 00:17:06,460 moving faster than the speed of sound, 384 00:17:06,460 --> 00:17:08,319 I can introduce any pressure fluctuation 385 00:17:08,319 --> 00:17:10,089 I want behind the shockwave, and it 386 00:17:10,089 --> 00:17:12,230 has no effect on the shockwave. 387 00:17:12,230 --> 00:17:14,920 Because the pressure waves from that pressure fluctuation 388 00:17:14,920 --> 00:17:15,760 are traveling at the speed of sound. 389 00:17:15,760 --> 00:17:17,410 They won't catch up to the shock. 390 00:17:17,410 --> 00:17:19,786 So they will have no effect whatsoever. 391 00:17:19,786 --> 00:17:21,369 So if I was judging what was happening 392 00:17:21,369 --> 00:17:23,202 by looking at what's happening in the shock, 393 00:17:23,202 --> 00:17:28,830 I can introduce any kind of random noise I want back over, 394 00:17:28,830 --> 00:17:30,260 downstream of the shock, I guess. 395 00:17:30,260 --> 00:17:31,050 Yeah? 396 00:17:31,050 --> 00:17:33,212 And it won't make any difference at all. 397 00:17:33,212 --> 00:17:34,836 And so, that's part of the reason why I 398 00:17:34,836 --> 00:17:37,110 need a special kind of solver. 399 00:17:37,110 --> 00:17:39,900 If I have elliptic equations, every point 400 00:17:39,900 --> 00:17:42,090 affects every other point. 401 00:17:42,090 --> 00:17:45,000 A famous case for that is, like, this steady state heat 402 00:17:45,000 --> 00:17:46,500 equation. 403 00:17:46,500 --> 00:17:50,310 You have some temperature source here, a heat source here, 404 00:17:50,310 --> 00:17:53,640 and you have some coolness source over here. 405 00:17:53,640 --> 00:17:57,585 And there's some heat flowing from the hot to the cold. 406 00:17:57,585 --> 00:18:00,104 And the heat flow is kind of slow enough 407 00:18:00,104 --> 00:18:02,520 that everything sort of affects everything else, since you 408 00:18:02,520 --> 00:18:04,890 actually have, the heat is trying to flow by, 409 00:18:04,890 --> 00:18:06,660 sort of, every path. 410 00:18:06,660 --> 00:18:10,282 And there's some insulation that's resisting its flow. 411 00:18:10,282 --> 00:18:11,740 And it has been in the connectivity 412 00:18:11,740 --> 00:18:12,869 of flow a certain way. 413 00:18:12,869 --> 00:18:14,910 And those kinds of equations are called elliptic. 414 00:18:14,910 --> 00:18:16,830 And the methods that we use for the relaxation mesh 415 00:18:16,830 --> 00:18:19,246 that we showed are really good for those kinds of methods, 416 00:18:19,246 --> 00:18:20,550 for those kind of problems. 417 00:18:20,550 --> 00:18:22,841 And then another kind of problem that we get into a lot 418 00:18:22,841 --> 00:18:23,890 is like this one. 419 00:18:23,890 --> 00:18:25,800 This is called parabolic. 420 00:18:25,800 --> 00:18:28,230 And typically, in the cases we see that are parabolic, 421 00:18:28,230 --> 00:18:32,390 we have time as part of the problem. 422 00:18:32,390 --> 00:18:34,374 Always what we think should happen 423 00:18:34,374 --> 00:18:35,790 is that what happens in the future 424 00:18:35,790 --> 00:18:40,170 depends on what happened in the past, but not vice versa. 425 00:18:40,170 --> 00:18:45,060 So you would expect that when you're 426 00:18:45,060 --> 00:18:48,662 computing what's happening at the next time point, 427 00:18:48,662 --> 00:18:50,870 it's going to depend on what happened at earlier time 428 00:18:50,870 --> 00:18:51,590 points. 429 00:18:51,590 --> 00:18:53,846 But if you try to do it the other way around, 430 00:18:53,846 --> 00:18:55,220 it would be kind of weird, right? 431 00:18:55,220 --> 00:18:57,011 You wouldn't expect that what in the future 432 00:18:57,011 --> 00:18:59,210 really affected stuff in the past. 433 00:18:59,210 --> 00:19:04,370 So it naturally leads us to sort of pose it as an ODE IVP, 434 00:19:04,370 --> 00:19:05,510 or IVP kind of problem. 435 00:19:05,510 --> 00:19:07,844 And we know initially, some time zero something, 436 00:19:07,844 --> 00:19:09,260 and we're trying to calculate what 437 00:19:09,260 --> 00:19:11,660 happens at some later time. 438 00:19:11,660 --> 00:19:14,210 And so those kinds of problems we call parabolic. 439 00:19:14,210 --> 00:19:15,700 And the methods that we're going to talk about mostly 440 00:19:15,700 --> 00:19:17,908 are ones for solving parabolic and elliptic problems. 441 00:19:21,120 --> 00:19:24,894 And you can read, in the middle of chapter 6, 442 00:19:24,894 --> 00:19:25,810 has a nice discussion. 443 00:19:25,810 --> 00:19:27,580 Just has the mathematical definitions 444 00:19:27,580 --> 00:19:28,621 of what these things are. 445 00:19:32,660 --> 00:19:36,710 Now, even in a steady state problem, 446 00:19:36,710 --> 00:19:39,860 the problem can have a sort of directionality to it. 447 00:19:39,860 --> 00:19:42,170 So if you have a problem that's very convective. 448 00:19:42,170 --> 00:19:43,820 So you have a flow in a pipe. 449 00:19:43,820 --> 00:19:46,362 And you have a nice, fast velocity flow down the pipe. 450 00:19:46,362 --> 00:19:48,320 What's happening downwind is very much affected 451 00:19:48,320 --> 00:19:51,750 by what happened upwind, but not vice versa. 452 00:19:51,750 --> 00:19:54,710 So maybe you're doing some reaction upwind. 453 00:19:54,710 --> 00:19:56,639 Say you're burning some hydrogen. 454 00:19:56,639 --> 00:19:58,430 You're going to get water vapor made upwind 455 00:19:58,430 --> 00:19:59,763 and it's going to flow downwind. 456 00:19:59,763 --> 00:20:03,440 And you're going to have steam downwind. 457 00:20:03,440 --> 00:20:06,117 If I squirt a little methane in downwind, 458 00:20:06,117 --> 00:20:08,450 it doesn't actually do anything to what happened upwind. 459 00:20:08,450 --> 00:20:10,491 I still have a hydrogen and oxygen flame up there 460 00:20:10,491 --> 00:20:12,260 and it's been burning and making steam. 461 00:20:12,260 --> 00:20:13,700 Do you see what I mean? 462 00:20:13,700 --> 00:20:15,470 So it has a directionality, even though I 463 00:20:15,470 --> 00:20:18,011 might do it as a steady state problem and not have time in it 464 00:20:18,011 --> 00:20:18,704 at all. 465 00:20:18,704 --> 00:20:20,870 And you've probably already done this, where you can 466 00:20:20,870 --> 00:20:23,300 convert, like, a flow reactor. 467 00:20:23,300 --> 00:20:26,464 You can write it, sort of, as a time thing or as a space thing. 468 00:20:26,464 --> 00:20:28,880 Time and space are so related by the velocity of the flow. 469 00:20:31,710 --> 00:20:34,050 We'll have similar kinds of phenomenon then. 470 00:20:34,050 --> 00:20:36,510 It's almost like a parabolic time system 471 00:20:36,510 --> 00:20:39,110 that, if the flow is fast enough, 472 00:20:39,110 --> 00:20:42,480 the fusion backup, anything moving upstream is negligible. 473 00:20:42,480 --> 00:20:45,640 It's all moving from upstream to downstream. 474 00:20:45,640 --> 00:20:49,220 And so then, you'll have the same kind of issues. 475 00:20:49,220 --> 00:20:51,060 In fact, you may even try to solve it 476 00:20:51,060 --> 00:20:54,750 by just starting it at the upstream end and competing 477 00:20:54,750 --> 00:20:57,872 stuff and then propagating what that output from that first one 478 00:20:57,872 --> 00:20:59,580 does to the next one downstream, and then 479 00:20:59,580 --> 00:21:00,538 do them, one at a time. 480 00:21:00,538 --> 00:21:02,130 That would be one possible way to try 481 00:21:02,130 --> 00:21:05,430 to solve those kinds of problems that are highly convective. 482 00:21:05,430 --> 00:21:08,572 As you slow the velocity down, then the diffusion 483 00:21:08,572 --> 00:21:10,530 will start to fight the velocity more and more. 484 00:21:10,530 --> 00:21:13,000 If you make the velocity very slow, 485 00:21:13,000 --> 00:21:15,030 the Peclet number very low, then you really 486 00:21:15,030 --> 00:21:17,520 have to worry about things diffusing back upstream. 487 00:21:17,520 --> 00:21:19,410 And then it's a different situation. 488 00:21:19,410 --> 00:21:23,260 And that would be more like the elliptic problems we have. 489 00:21:23,260 --> 00:21:25,740 Because things downstream would actually affect upstream, 490 00:21:25,740 --> 00:21:27,090 if the flow is really small. 491 00:21:27,090 --> 00:21:28,400 But typically for our problems, the velocity 492 00:21:28,400 --> 00:21:29,580 is always really high. 493 00:21:29,580 --> 00:21:32,434 Peclet numbers are really large and so not much goes 494 00:21:32,434 --> 00:21:33,600 from downstream to upstream. 495 00:21:33,600 --> 00:21:38,030 It's almost all upstream to downstream. 496 00:21:38,030 --> 00:21:41,235 This leads to a famous situation. 497 00:21:41,235 --> 00:21:50,627 If you look in the textbook at figures 6.7, 6.8, 498 00:21:50,627 --> 00:21:52,085 it's kind of a very famous problem. 499 00:21:54,730 --> 00:22:03,640 That if you just do this problem. 500 00:22:23,619 --> 00:22:24,910 So this is just an ODE problem. 501 00:22:29,760 --> 00:22:32,730 So I just have convection and diffusion, no reaction, 502 00:22:32,730 --> 00:22:33,900 very simple. 503 00:22:33,900 --> 00:22:35,070 This problem. 504 00:22:35,070 --> 00:22:38,100 You'd think that I could just use my finite differences, 505 00:22:38,100 --> 00:22:40,540 where I would put, I could use the center difference 506 00:22:40,540 --> 00:22:47,860 for here and here and change this into this obvious looking 507 00:22:47,860 --> 00:22:48,360 thing. 508 00:23:05,669 --> 00:23:08,210 You think that that would be a reasonable thing to do, right? 509 00:23:08,210 --> 00:23:09,380 And you can do the same one here. 510 00:23:09,380 --> 00:23:10,004 You can write-- 511 00:23:17,556 --> 00:23:18,930 I'm terrible with minuses, sorry. 512 00:23:27,620 --> 00:23:30,900 So these, the finite difference formulas. 513 00:23:30,900 --> 00:23:34,640 And so, you can write these like this. 514 00:23:34,640 --> 00:23:36,430 This stays equal to zero. 515 00:23:36,430 --> 00:23:39,640 And now it's just an algebraic problem to solve. 516 00:23:39,640 --> 00:23:41,530 A system of equations, you have equations 517 00:23:41,530 --> 00:23:43,430 like this for every value of n. 518 00:23:43,430 --> 00:23:46,290 And there's the mesh points. 519 00:23:46,290 --> 00:23:48,420 You guys have done this before, yes? 520 00:23:48,420 --> 00:23:49,170 Yes. 521 00:23:49,170 --> 00:23:51,870 OK, so famous problem is, if you actually 522 00:23:51,870 --> 00:23:57,740 try to solve it this way, unless you make delta z really tiny, 523 00:23:57,740 --> 00:23:58,880 you get crazy stuff. 524 00:23:58,880 --> 00:24:02,210 So on figure 6.7, they actually show 525 00:24:02,210 --> 00:24:08,780 what happens for different values of delta z, which 526 00:24:08,780 --> 00:24:11,930 is different values of what they call the local Peclet number. 527 00:24:11,930 --> 00:24:18,430 So there's a thing called local Peclet number, which is-- 528 00:24:23,680 --> 00:24:26,090 OK? 529 00:24:26,090 --> 00:24:29,120 If you make the local Peclet number too large, actually 530 00:24:29,120 --> 00:24:30,950 anywhere bigger than 2, and you try 531 00:24:30,950 --> 00:24:32,450 to solve this equation, what you get 532 00:24:32,450 --> 00:24:35,210 is oscillations, unphysical oscillations. 533 00:24:35,210 --> 00:24:37,800 They'll make the phi go negative. 534 00:24:37,800 --> 00:24:38,310 It's crazy. 535 00:24:41,107 --> 00:24:43,440 It looks like the simplest equation in the world, right? 536 00:24:43,440 --> 00:24:47,320 It's a linear equation, so what's the problem with this? 537 00:24:47,320 --> 00:24:49,250 But it doesn't work. 538 00:24:49,250 --> 00:24:53,590 And so, then you could think, well, why doesn't it work? 539 00:24:53,590 --> 00:24:55,660 And there's, like, multiple ways to explain this 540 00:24:55,660 --> 00:24:58,249 about why this doesn't work. 541 00:24:58,249 --> 00:25:00,040 I think that the best way to think about it 542 00:25:00,040 --> 00:25:01,400 is about the information flow. 543 00:25:01,400 --> 00:25:04,400 So if the local Peclet number is large, 544 00:25:04,400 --> 00:25:08,650 that means that phi is just flowing downwind 545 00:25:08,650 --> 00:25:10,660 and diffusion is not doing much. 546 00:25:10,660 --> 00:25:12,920 Because the Peclet number is big. 547 00:25:12,920 --> 00:25:18,140 And so, you really shouldn't use this formula 548 00:25:18,140 --> 00:25:19,970 to calculate the flow. 549 00:25:19,970 --> 00:25:24,769 Because what happens downwind, at point n plus one, I guess. 550 00:25:24,769 --> 00:25:26,060 Here we have the flow this way. 551 00:25:29,240 --> 00:25:32,520 So here's n minus one. 552 00:25:32,520 --> 00:25:34,770 Here's my point n, n plus 1. 553 00:25:38,470 --> 00:25:42,280 I really should not include, is this right, 554 00:25:42,280 --> 00:25:44,991 n plus 1 in this formula. 555 00:25:44,991 --> 00:25:46,240 Hope I have those signs right. 556 00:25:46,240 --> 00:25:49,266 I apologize if I had them backwards. 557 00:25:49,266 --> 00:25:50,890 Because I don't think that what happens 558 00:25:50,890 --> 00:25:54,320 downstream should really affect this at all. 559 00:25:54,320 --> 00:25:57,070 So really, I would do better to change from this center 560 00:25:57,070 --> 00:26:01,910 difference formula to what's called the upwind difference 561 00:26:01,910 --> 00:26:02,410 formula. 562 00:26:09,280 --> 00:26:10,930 Where they would say, OK, instead, 563 00:26:10,930 --> 00:26:23,810 let's use phi of n minus, the phi of n minus 1 over delta z. 564 00:26:23,810 --> 00:26:27,346 Just use that instead of using this formula. 565 00:26:27,346 --> 00:26:28,720 Now, you wouldn't think that that 566 00:26:28,720 --> 00:26:30,742 would make much difference. 567 00:26:30,742 --> 00:26:32,200 But it makes a gigantic difference. 568 00:26:32,200 --> 00:26:35,290 And so if you look at figure 6.8 in the textbook, 569 00:26:35,290 --> 00:26:40,180 you see you get perfectly stable solutions when you do that. 570 00:26:40,180 --> 00:26:41,740 But where you do this one, you get 571 00:26:41,740 --> 00:26:43,480 crazy, oscillatory, unphysical solutions, 572 00:26:43,480 --> 00:26:47,030 unless you choose delta z really tiny. 573 00:26:47,030 --> 00:26:50,487 Now, there's other ways to look at this. 574 00:26:50,487 --> 00:26:52,320 So I think this is the correct way, the best 575 00:26:52,320 --> 00:26:55,210 way to think about it is, it has to do with information flow. 576 00:26:55,210 --> 00:26:57,210 If you try to make your solution depend on stuff 577 00:26:57,210 --> 00:26:58,666 that it doesn't really depend on. 578 00:26:58,666 --> 00:27:00,540 So upwind does not really depend on downwind. 579 00:27:00,540 --> 00:27:03,360 If you force it to, by choosing to include 580 00:27:03,360 --> 00:27:07,050 that information here, you end up with crazy stuff. 581 00:27:07,050 --> 00:27:09,610 So that's one way to look at it. 582 00:27:09,610 --> 00:27:12,000 Another way to look at it is that, if you really 583 00:27:12,000 --> 00:27:14,260 want to compute these derivatives, 584 00:27:14,260 --> 00:27:15,960 you've got to keep delta z small. 585 00:27:15,960 --> 00:27:17,640 Because we're taking a differential. 586 00:27:17,640 --> 00:27:20,010 We're turning it into a finite difference. 587 00:27:20,010 --> 00:27:22,640 And if you take your delta z too big, 588 00:27:22,640 --> 00:27:27,680 then it's a terrible approximation to do this. 589 00:27:27,680 --> 00:27:28,639 And you can look at it. 590 00:27:28,639 --> 00:27:30,013 You can see that what we're doing 591 00:27:30,013 --> 00:27:32,219 is, when we have this problem where the local Peclet 592 00:27:32,219 --> 00:27:34,510 number is getting too large, is we're actually choosing 593 00:27:34,510 --> 00:27:35,885 delta z bigger than the, sort of, 594 00:27:35,885 --> 00:27:38,340 the thickness of the solution. 595 00:27:38,340 --> 00:27:41,270 So if you look at the solution, the analytical solution 596 00:27:41,270 --> 00:27:44,590 of this problem, it's like this. 597 00:27:52,190 --> 00:27:54,930 Where this has a very, very thin layer at the end. 598 00:27:54,930 --> 00:27:58,720 It's like the flow has pushed all your stuff to one side. 599 00:27:58,720 --> 00:28:03,520 And if you choose your delta z to be from here to here, 600 00:28:03,520 --> 00:28:05,530 that's how big your delta z is, and then you're 601 00:28:05,530 --> 00:28:11,680 computing the derivative of this point halfway along from here 602 00:28:11,680 --> 00:28:13,395 to where this is. 603 00:28:13,395 --> 00:28:15,180 You compute your derivative like that. 604 00:28:15,180 --> 00:28:18,520 It's not a very good, not a very accurate representation of what 605 00:28:18,520 --> 00:28:19,912 the real derivative is here. 606 00:28:19,912 --> 00:28:21,370 And in fact, you really should have 607 00:28:21,370 --> 00:28:23,411 a bunch of points here and here and here and here 608 00:28:23,411 --> 00:28:26,789 and here to really resolve the sharp edge. 609 00:28:26,789 --> 00:28:28,330 And so you did something really crazy 610 00:28:28,330 --> 00:28:31,260 to use a big value of delta z to start with. 611 00:28:31,260 --> 00:28:34,070 So that's another way to look at this problem. 612 00:28:34,070 --> 00:28:38,062 Now, the fix, by doing this upwind differencing, 613 00:28:38,062 --> 00:28:40,020 I think the best way to look at this is saying, 614 00:28:40,020 --> 00:28:42,020 well, I'm just going to make sure I only depend, 615 00:28:42,020 --> 00:28:45,069 I only make the equations depend on what's upwind, 616 00:28:45,069 --> 00:28:46,110 with the convective term. 617 00:28:46,110 --> 00:28:47,850 Because there's nothing convecting from downwind. 618 00:28:47,850 --> 00:28:48,930 So just leave that out. 619 00:28:48,930 --> 00:28:51,270 So that's a conceptual way to think about it. 620 00:28:51,270 --> 00:28:53,760 Another way to think about it, which is in the text, 621 00:28:53,760 --> 00:28:57,300 is that by doing this, you're introducing 622 00:28:57,300 --> 00:28:59,550 another thing called numerical diffusion, where 623 00:28:59,550 --> 00:29:03,120 you're actually, effectively increasing the diffusivity. 624 00:29:03,120 --> 00:29:05,317 Because this is a very poor approximation 625 00:29:05,317 --> 00:29:06,150 to the differential. 626 00:29:06,150 --> 00:29:10,300 This is a asymmetrical finite difference formula. 627 00:29:10,300 --> 00:29:13,260 It's not a very good value estimate of the derivative, 628 00:29:13,260 --> 00:29:15,930 if delta z is big. 629 00:29:15,930 --> 00:29:18,829 But you can write it out carefully and write this out 630 00:29:18,829 --> 00:29:20,370 and say, oh, this actually is sort of 631 00:29:20,370 --> 00:29:24,690 like this plus an effective diffusivity 632 00:29:24,690 --> 00:29:27,960 chosen very carefully, so the terms cancel out just right. 633 00:29:27,960 --> 00:29:29,490 And it turns out, if you look at it 634 00:29:29,490 --> 00:29:31,650 that way, the effective diffusion you've 635 00:29:31,650 --> 00:29:34,320 added, by using this instead of that, 636 00:29:34,320 --> 00:29:38,920 is just enough to make the local Peclet number, 637 00:29:38,920 --> 00:29:41,440 if you compute it with the effects of this d effective, 638 00:29:41,440 --> 00:29:42,940 to stay less than 2. 639 00:29:42,940 --> 00:29:45,130 Which is what you need as a condition to make 640 00:29:45,130 --> 00:29:46,130 this numerically stable. 641 00:29:49,040 --> 00:29:51,022 So I strongly suggest you read that. 642 00:29:51,022 --> 00:29:52,730 It's only two pages long in the textbook. 643 00:29:52,730 --> 00:29:53,180 It's very clear. 644 00:29:53,180 --> 00:29:55,638 Just read, just look at it if you haven't seen this before. 645 00:30:01,347 --> 00:30:03,930 There's a similar thing and it might be relevant to the COMSOL 646 00:30:03,930 --> 00:30:09,780 problem that Kristyn showed you on Monday, is you can-- 647 00:30:12,450 --> 00:30:17,430 In that case, there was a region of the boundary here 648 00:30:17,430 --> 00:30:20,650 that had some concentration c naught, 649 00:30:20,650 --> 00:30:24,790 and then there was, over here, the concentration was zero. 650 00:30:24,790 --> 00:30:25,890 The boundary is zero. 651 00:30:25,890 --> 00:30:31,420 Here, it's some number, 17, whatever number it was. 652 00:30:31,420 --> 00:30:39,470 And if you think of how to model what's going on, 653 00:30:39,470 --> 00:30:44,120 one way to look at it is, I have a diffusive flux coming in 654 00:30:44,120 --> 00:30:48,230 from the drug patch, diffusing the drug into the flow. 655 00:30:48,230 --> 00:30:52,870 And I have this big, giant flow, flowing along here. 656 00:30:52,870 --> 00:31:00,900 And the flow gets slower as I get closer to the wall because 657 00:31:00,900 --> 00:31:03,250 of the friction with the wall. 658 00:31:03,250 --> 00:31:06,290 And so if I looked somewhere right in here, 659 00:31:06,290 --> 00:31:08,170 I could just draw a little control volume, 660 00:31:08,170 --> 00:31:09,280 and look what's happening. 661 00:31:09,280 --> 00:31:11,262 I have a diffusive flux coming in here. 662 00:31:11,262 --> 00:31:12,970 I don't have any drug up here, so there's 663 00:31:12,970 --> 00:31:13,970 no drug coming in here. 664 00:31:13,970 --> 00:31:15,770 But I have drug leaving. 665 00:31:15,770 --> 00:31:19,950 So basically, it's coming up and it's making a right-hand turn. 666 00:31:19,950 --> 00:31:21,927 And so then I could think, well, if I 667 00:31:21,927 --> 00:31:23,510 was going to try to figure out, what's 668 00:31:23,510 --> 00:31:26,250 the steady state concentration of drug in there, I could say, 669 00:31:26,250 --> 00:31:27,120 well, there's a certain amount of drug 670 00:31:27,120 --> 00:31:28,744 entering from the diffusive flux, which 671 00:31:28,744 --> 00:31:30,885 would be just this times dcdy. 672 00:31:34,320 --> 00:31:35,741 This is the y direction. 673 00:31:38,970 --> 00:31:41,280 And then this part over here would 674 00:31:41,280 --> 00:31:46,720 be the velocity, actually, just times the concentration 675 00:31:46,720 --> 00:31:47,900 I have in here, right? 676 00:31:47,900 --> 00:31:52,160 So how much is flowing out. 677 00:31:52,160 --> 00:31:54,570 And it's probably, my units are screwed up so there's-- 678 00:31:54,570 --> 00:31:55,700 Well, maybe not. 679 00:31:55,700 --> 00:31:56,200 That's OK. 680 00:31:56,200 --> 00:31:57,610 Is that all right? 681 00:32:00,430 --> 00:32:01,725 Units are screwed up. 682 00:32:01,725 --> 00:32:02,850 You guys will get it right. 683 00:32:05,360 --> 00:32:06,946 It's probably an area. 684 00:32:06,946 --> 00:32:08,770 Something to do with this size right here. 685 00:32:08,770 --> 00:32:09,561 This one here, too. 686 00:32:14,640 --> 00:32:17,570 Is that right? 687 00:32:17,570 --> 00:32:20,747 So you could compute what the steady state 688 00:32:20,747 --> 00:32:23,330 concentration would be if it's just coming in and flowing out. 689 00:32:23,330 --> 00:32:25,320 And so, this is the finite volume view 690 00:32:25,320 --> 00:32:27,620 of this kind of problem. 691 00:32:27,620 --> 00:32:30,464 And so you can write the equations that way, instead. 692 00:32:30,464 --> 00:32:32,880 And what you really care about are what the velocity flows 693 00:32:32,880 --> 00:32:34,879 and the diffusive fluxes are on these boundaries 694 00:32:34,879 --> 00:32:36,040 around the point. 695 00:32:36,040 --> 00:32:38,550 So you don't try to compute things right at the point. 696 00:32:38,550 --> 00:32:39,940 You compute around it. 697 00:32:39,940 --> 00:32:42,200 And one way to look at this is, we just 698 00:32:42,200 --> 00:32:43,310 did finite differencing. 699 00:32:43,310 --> 00:32:45,750 We were looking at what was happening coming in. 700 00:32:49,390 --> 00:32:52,160 And then, this is actually a center difference 701 00:32:52,160 --> 00:32:56,360 around this interface, which is halfway between the two points. 702 00:32:56,360 --> 00:32:59,450 And so actually, this is a good approximation 703 00:32:59,450 --> 00:33:01,925 for the finite volume point of view. 704 00:33:01,925 --> 00:33:04,090 This is all right? 705 00:33:04,090 --> 00:33:05,810 Now, these are all different ways 706 00:33:05,810 --> 00:33:07,220 to write down the equations. 707 00:33:07,220 --> 00:33:10,100 All of them are correct in the limit 708 00:33:10,100 --> 00:33:11,694 that delta z goes to zero. 709 00:33:11,694 --> 00:33:13,610 If you have an infinite number of mesh points, 710 00:33:13,610 --> 00:33:15,350 all these are exactly the same. 711 00:33:15,350 --> 00:33:18,710 But they lead to really different numerical properties 712 00:33:18,710 --> 00:33:21,360 when delta z is large. 713 00:33:21,360 --> 00:33:23,910 And they're all inaccurate when delta z is large, 714 00:33:23,910 --> 00:33:26,390 so it's all about what happens with the inaccuracies. 715 00:33:30,910 --> 00:33:35,650 Now, we can go ahead and take any problem, even 716 00:33:35,650 --> 00:33:38,980 really complicated 3-D problems or 4-D problems like this, 717 00:33:38,980 --> 00:33:42,164 and we can write them with relaxation methods. 718 00:33:42,164 --> 00:33:43,955 And what we're basically doing in a problem 719 00:33:43,955 --> 00:33:48,430 like this is turning it into some function of y 720 00:33:48,430 --> 00:33:53,460 is equal to zero, where y is the value of all the unknowns. 721 00:33:57,320 --> 00:34:01,230 And there's a lot of equations. 722 00:34:01,230 --> 00:34:04,110 So I might have 100 million equations 723 00:34:04,110 --> 00:34:08,130 and 100 million unknowns that are the values, 724 00:34:08,130 --> 00:34:12,179 or every state variable at every state point at all times. 725 00:34:12,179 --> 00:34:14,389 And I can write it down, by finding differences, 726 00:34:14,389 --> 00:34:16,770 by finding elements, by finding volumes. 727 00:34:16,770 --> 00:34:20,239 Any way I want, I can write down an expression like this. 728 00:34:20,239 --> 00:34:21,929 And I know in the limit, as I make 729 00:34:21,929 --> 00:34:24,370 the number of elements of y go to infinity, 730 00:34:24,370 --> 00:34:27,567 it will actually be the real solution. 731 00:34:27,567 --> 00:34:29,400 And the problem is, my computer can't handle 732 00:34:29,400 --> 00:34:30,870 infinite number of unknowns. 733 00:34:30,870 --> 00:34:32,536 In fact, it's even going to have trouble 734 00:34:32,536 --> 00:34:33,719 with a 100 million unknowns. 735 00:34:33,719 --> 00:34:35,622 And so, I have to figure out what to do. 736 00:34:35,622 --> 00:34:37,080 So how would I solve this normally, 737 00:34:37,080 --> 00:34:39,570 is I would take the Jacobian of this and I would do, 738 00:34:39,570 --> 00:34:42,310 I'd take an initial guess and I'd update it and have 739 00:34:42,310 --> 00:34:48,889 the Jacobian times my change from the initial guess 740 00:34:48,889 --> 00:34:52,424 to the equal negative of f, right? 741 00:34:52,424 --> 00:34:54,090 You guys have done this a million times. 742 00:34:54,090 --> 00:34:57,030 As Newton-Raphson is how you'd find 743 00:34:57,030 --> 00:34:58,675 improved from an initial guess. 744 00:34:58,675 --> 00:35:00,300 You have to provide some initial guess. 745 00:35:00,300 --> 00:35:02,110 This is actually pretty hard if you have 100 million unknowns. 746 00:35:02,110 --> 00:35:05,000 You have to provide the value of the initial guess for every one 747 00:35:05,000 --> 00:35:06,300 of the 100 million. 748 00:35:06,300 --> 00:35:08,309 But somehow you did it. 749 00:35:08,309 --> 00:35:10,350 And now you want to refine that guess and this is 750 00:35:10,350 --> 00:35:12,960 the formula for it. 751 00:35:12,960 --> 00:35:15,000 And so, you just use backslash, right? 752 00:35:15,000 --> 00:35:19,320 But the problem here is that now f is a 100 million long vector 753 00:35:19,320 --> 00:35:21,970 and this is a 100 million squared matrix. 754 00:35:21,970 --> 00:35:23,820 And a 100 million squared is pretty big, 755 00:35:23,820 --> 00:35:26,050 10 to 16th elements in it. 756 00:35:26,050 --> 00:35:28,050 So I'm not going to just write that matrix 757 00:35:28,050 --> 00:35:30,060 down directly like this. 758 00:35:30,060 --> 00:35:33,630 Now, we cleverly chose local basis functions, 759 00:35:33,630 --> 00:35:37,480 so this is going to be sparse, this Jacobian. 760 00:35:37,480 --> 00:35:39,950 So most of the numbers, are those 10 to the 16th numbers 761 00:35:39,950 --> 00:35:41,942 in this matrix are zero. 762 00:35:41,942 --> 00:35:43,400 So we don't have to represent them. 763 00:35:43,400 --> 00:35:45,370 But there still might be quite a few. 764 00:35:45,370 --> 00:35:47,120 Probably the diagonal events are non-zero, 765 00:35:47,120 --> 00:35:49,530 so that's like 10 to the 8th of them right there. 766 00:35:49,530 --> 00:35:50,930 And we've got, probably, a couple of other bands. 767 00:35:50,930 --> 00:35:53,390 So I don't know, might be 10 to the 9th non-zero elements 768 00:35:53,390 --> 00:35:56,240 inside this thing, which is pretty nutty. 769 00:35:56,240 --> 00:35:58,790 And depending on how much memory they provided in the laptop 770 00:35:58,790 --> 00:36:01,370 they gave you, you might be filling up the memory 771 00:36:01,370 --> 00:36:02,961 already just to write down j. 772 00:36:02,961 --> 00:36:05,210 And certainly, if you tried to do Gaussian elimination 773 00:36:05,210 --> 00:36:08,477 to solve this, you're going to have a problem because you 774 00:36:08,477 --> 00:36:09,560 get what's called fill in. 775 00:36:09,560 --> 00:36:12,510 Do you guys remember this? 776 00:36:12,510 --> 00:36:14,900 Even if you start from a very sparse Jacobian, 777 00:36:14,900 --> 00:36:17,360 as you run the Gaussian elimination steps 778 00:36:17,360 --> 00:36:19,850 the number of non-zero entries in the intermediate matrices 779 00:36:19,850 --> 00:36:21,180 gets larger and larger. 780 00:36:21,180 --> 00:36:22,359 Remember this? 781 00:36:22,359 --> 00:36:24,900 And so, even if, initially, you can write down j and store it 782 00:36:24,900 --> 00:36:26,710 in your computer, you could easily 783 00:36:26,710 --> 00:36:31,330 overwhelm it with the third iteration or something. 784 00:36:31,330 --> 00:36:33,560 So you're not going to just use Gaussian elimination. 785 00:36:33,560 --> 00:36:34,580 So what do you have to do? 786 00:36:34,580 --> 00:36:35,663 You want to find a method. 787 00:36:35,663 --> 00:36:38,518 They're called direct. 788 00:36:38,518 --> 00:36:43,090 Direct methods to solve this kind of problem, 789 00:36:43,090 --> 00:36:44,730 in order to get your increment, which 790 00:36:44,730 --> 00:36:47,260 is going to be your update, to get a better estimate of what 791 00:36:47,260 --> 00:36:49,050 the solution is. 792 00:36:49,050 --> 00:36:51,870 And the direct methods are ones that you don't actually 793 00:36:51,870 --> 00:36:53,940 store the gigantic matrix. 794 00:36:53,940 --> 00:36:57,260 So you're trading off CPU time for memory. 795 00:36:57,260 --> 00:37:00,450 So you don't have enough memory to store the whole thing. 796 00:37:00,450 --> 00:37:03,810 Which we, normally you would do this with Gaussian elimination 797 00:37:03,810 --> 00:37:04,760 and huge steps. 798 00:37:04,760 --> 00:37:06,712 It would solve it perfectly, right? 799 00:37:06,712 --> 00:37:08,670 Right, Gaussian elimination's great, backslash. 800 00:37:08,670 --> 00:37:10,335 You guys have used it a few times. 801 00:37:10,335 --> 00:37:11,460 It's a really nice program. 802 00:37:11,460 --> 00:37:12,660 It always works. 803 00:37:12,660 --> 00:37:14,660 It's really nice. 804 00:37:14,660 --> 00:37:17,160 But you're going to give up that because you can't afford it 805 00:37:17,160 --> 00:37:19,260 because the memory, it's consuming too much memory 806 00:37:19,260 --> 00:37:21,093 and your cheapo department didn't buy enough 807 00:37:21,093 --> 00:37:23,940 RAM in your computer for you. 808 00:37:23,940 --> 00:37:27,540 So you're going to instead try a direct method which 809 00:37:27,540 --> 00:37:37,830 is, so it's trading off CPU time versus RAM. 810 00:37:37,830 --> 00:37:40,080 So you're going to reduce how much memory you actually 811 00:37:40,080 --> 00:37:42,754 actively have, but you're going to expect that it's 812 00:37:42,754 --> 00:37:43,920 going to cost more CPU time. 813 00:37:43,920 --> 00:37:45,590 Because otherwise, everybody would do something else besides 814 00:37:45,590 --> 00:37:48,360 Gaussian elimination all the time, and they don't. 815 00:37:48,360 --> 00:37:52,080 All right, so what's the method we're going to use? 816 00:37:52,080 --> 00:37:55,260 It turns out that variance on the conjugate gradient method 817 00:37:55,260 --> 00:37:56,310 turned out really well. 818 00:37:56,310 --> 00:38:00,510 So we mentioned this earlier, briefly. 819 00:38:00,510 --> 00:38:05,920 What you do is, instead of trying to solve this, 820 00:38:05,920 --> 00:38:08,140 you try to solve this. 821 00:38:26,060 --> 00:38:28,660 So you just try to minimize that, right? 822 00:38:28,660 --> 00:38:32,310 You know at the solution this is going to be zero. 823 00:38:32,310 --> 00:38:35,800 Right, jy plus [INAUDIBLE] zero when this is solved. 824 00:38:35,800 --> 00:38:38,740 So you try and vary and find the delta y 825 00:38:38,740 --> 00:38:42,390 that makes this go to zero by minimizing it. 826 00:38:42,390 --> 00:38:44,070 And then this method, it turns out 827 00:38:44,070 --> 00:38:47,640 that the conjugate gradient, if everything works great, 828 00:38:47,640 --> 00:38:50,220 this is guaranteed in n iterations 829 00:38:50,220 --> 00:38:53,820 to go directly to the solution of this kind of problem. 830 00:38:53,820 --> 00:38:56,320 Now n is large, now, because we have 10 to the 8th unknowns. 831 00:38:56,320 --> 00:38:58,190 So that's 10 to the 8th iterations. 832 00:38:58,190 --> 00:38:59,330 You do anything in 10 to the 8th iterations, 833 00:38:59,330 --> 00:39:01,300 you're not really sure it's really going to work. 834 00:39:01,300 --> 00:39:02,133 So that's a warning. 835 00:39:02,133 --> 00:39:05,950 But at least it should get, even after a few iterations, 836 00:39:05,950 --> 00:39:10,270 you should get a smaller value of this than you had before. 837 00:39:10,270 --> 00:39:11,380 And so this is the method. 838 00:39:11,380 --> 00:39:12,520 And what's really nice about this method, 839 00:39:12,520 --> 00:39:14,470 if you look into the details of how it works, 840 00:39:14,470 --> 00:39:17,740 it never has to store any intermediate matrices. 841 00:39:17,740 --> 00:39:20,650 So all it needs is the capability 842 00:39:20,650 --> 00:39:24,730 to multiply your matrix times some vector. 843 00:39:24,730 --> 00:39:29,970 And from that multiplication, it figures out a step direction. 844 00:39:29,970 --> 00:39:33,370 That's the trick of this method. 845 00:39:33,370 --> 00:39:36,422 And because it only has to do a forward multiplication, 846 00:39:36,422 --> 00:39:37,880 you don't have to actually store j. 847 00:39:37,880 --> 00:39:40,171 You can just compute elements of j, as needed, in order 848 00:39:40,171 --> 00:39:43,240 to evaluate this top product. 849 00:39:43,240 --> 00:39:46,730 This j times v. And you don't have to ever store 850 00:39:46,730 --> 00:39:48,390 the whole j at once. 851 00:39:48,390 --> 00:39:50,620 So it's really great for RAM to do this method. 852 00:39:53,330 --> 00:39:56,350 Now, as I said, when you do 10 to the 8th iterations, 853 00:39:56,350 --> 00:39:58,400 things don't always work out so well. 854 00:39:58,400 --> 00:40:00,220 And so people have found that when you go more than maybe 855 00:40:00,220 --> 00:40:01,844 10 or 20 iterations like this, you tend 856 00:40:01,844 --> 00:40:04,300 to pick up numerical problems. 857 00:40:04,300 --> 00:40:06,527 And so, they've worked out better methods. 858 00:40:06,527 --> 00:40:08,860 And there's one that people use a lot for these problems 859 00:40:08,860 --> 00:40:12,820 called bi cg stab. 860 00:40:12,820 --> 00:40:16,800 Which is biconjugate gradient stabilized. 861 00:40:16,800 --> 00:40:19,540 The applied mathematicians went crazy trying to figure out 862 00:40:19,540 --> 00:40:20,617 better methods like this. 863 00:40:20,617 --> 00:40:22,450 And this is the one that, I think currently, 864 00:40:22,450 --> 00:40:23,570 people think is the best. 865 00:40:23,570 --> 00:40:25,040 Though probably, there's a lot of research on this, 866 00:40:25,040 --> 00:40:26,623 so maybe there's even better ones now. 867 00:40:26,623 --> 00:40:29,870 But inside Matlab, I think this is the best on they have. 868 00:40:29,870 --> 00:40:34,770 And this is a method for doing this. 869 00:40:34,770 --> 00:40:35,830 Now, it's iterative. 870 00:40:35,830 --> 00:40:37,371 Now, let's remember what we're doing. 871 00:40:37,371 --> 00:40:38,910 We're trying to solve this problem. 872 00:40:38,910 --> 00:40:39,910 We started with a guess. 873 00:40:39,910 --> 00:40:42,535 We're going to break it up into an iterative problem like this, 874 00:40:42,535 --> 00:40:44,520 where we're going to do some steps, delta y. 875 00:40:44,520 --> 00:40:47,910 This is now doing iterative procedure 876 00:40:47,910 --> 00:40:49,710 in order to compute delta y. 877 00:40:49,710 --> 00:40:52,470 So we have an interim procedure and another iterative procedure 878 00:40:52,470 --> 00:40:54,499 inside the first iterative procedure. 879 00:40:54,499 --> 00:40:56,290 So this is going to be a lot more CPU time. 880 00:40:56,290 --> 00:40:57,940 So just to warn you. 881 00:41:01,024 --> 00:41:02,440 And it's also, it's not guaranteed 882 00:41:02,440 --> 00:41:04,356 like Gaussian elimination, to beautifully come 883 00:41:04,356 --> 00:41:06,097 to machine precision solution. 884 00:41:06,097 --> 00:41:07,430 It's going to get you something. 885 00:41:10,282 --> 00:41:11,740 But anyway, this is what people do. 886 00:41:11,740 --> 00:41:15,760 So that's one good thing to try. 887 00:41:15,760 --> 00:41:20,860 Now, this whole approach is sort of based on Newton-Raphson. 888 00:41:20,860 --> 00:41:24,310 We know that doesn't have the greatest radius of convergence. 889 00:41:24,310 --> 00:41:28,520 So you need to provide a pretty darn good initial guess. 890 00:41:28,520 --> 00:41:30,770 So how are we going to get a good initial guess if you 891 00:41:30,770 --> 00:41:32,210 have a 100 million unknowns? 892 00:41:35,250 --> 00:41:36,410 So what methods do we know? 893 00:41:36,410 --> 00:41:37,826 So one idea is you could do things 894 00:41:37,826 --> 00:41:40,210 like homotopy, stuff like that. 895 00:41:40,210 --> 00:41:42,830 If you could somehow take your PDE system, 896 00:41:42,830 --> 00:41:44,300 get rid of the nonlinear terms. 897 00:41:44,300 --> 00:41:47,741 You're really good at solving linear systems of PDEs. 898 00:41:47,741 --> 00:41:49,490 You start with that and then you gradually 899 00:41:49,490 --> 00:41:51,162 turn the non-linear term on. 900 00:41:51,162 --> 00:41:52,370 Maybe you could coax it over. 901 00:41:52,370 --> 00:41:56,180 So that's one possibility, if you're really good at PDEs. 902 00:41:56,180 --> 00:41:58,400 Another possibility that I've run 903 00:41:58,400 --> 00:42:03,140 into a lot is the problem you're trying to solve, for example, 904 00:42:03,140 --> 00:42:04,890 might be a steady state problem like this. 905 00:42:04,890 --> 00:42:07,040 Where this is zero and you try to figure out what 906 00:42:07,040 --> 00:42:11,900 the steady state situation is in a flow reactor, for example. 907 00:42:11,900 --> 00:42:14,240 And so, it's actually, your steady state problem 908 00:42:14,240 --> 00:42:17,010 came from a time dependent problem. 909 00:42:17,010 --> 00:42:19,640 And if you think that your time dependent problem really 910 00:42:19,640 --> 00:42:22,970 converges to the steady state problem, 911 00:42:22,970 --> 00:42:24,851 sort of without much concern about what 912 00:42:24,851 --> 00:42:26,350 the initial guess is, then you could 913 00:42:26,350 --> 00:42:28,550 put any random initial guess in and march it along 914 00:42:28,550 --> 00:42:32,300 in time enough, and eventually, it should end up 915 00:42:32,300 --> 00:42:35,420 at the solution you want. 916 00:42:35,420 --> 00:42:36,320 So that's the idea. 917 00:42:36,320 --> 00:42:37,903 So that's called a time marching idea. 918 00:42:37,903 --> 00:42:40,790 So there's one idea, sort of the homotopy continuation 919 00:42:40,790 --> 00:42:41,420 kind of idea. 920 00:42:41,420 --> 00:42:44,150 That's one idea about what to do for initial guesses. 921 00:42:44,150 --> 00:42:46,405 Another idea is the time marching idea. 922 00:42:46,405 --> 00:42:48,530 And the key about it is, you have to really believe 923 00:42:48,530 --> 00:42:50,900 that, physically, the system does march 924 00:42:50,900 --> 00:42:52,510 to the steady state you want. 925 00:42:52,510 --> 00:42:54,440 If it does, then this is a good idea. 926 00:42:54,440 --> 00:42:57,470 If it doesn't, then you're not going to end up where you want, 927 00:42:57,470 --> 00:42:59,990 because it's not going where you want. 928 00:42:59,990 --> 00:43:02,120 And so, you have to watch out. 929 00:43:02,120 --> 00:43:04,980 But in some situations, this is true, that things work. 930 00:43:04,980 --> 00:43:07,220 So I do a lot of combustion problems. 931 00:43:07,220 --> 00:43:11,390 A lot of these things, you light a spark on a stove burner. 932 00:43:11,390 --> 00:43:13,610 No matter how you light it, it ends up 933 00:43:13,610 --> 00:43:16,156 basically the same flame. 934 00:43:16,156 --> 00:43:17,030 So that's a good one. 935 00:43:17,030 --> 00:43:20,840 But I have some other ones where I light a candle 936 00:43:20,840 --> 00:43:23,840 and I'm in a convective field, a turbulent field, 937 00:43:23,840 --> 00:43:25,130 and it flickers. 938 00:43:25,130 --> 00:43:26,755 And no matter what, it always flickers. 939 00:43:26,755 --> 00:43:28,588 And I never get a steady state, so I'm never 940 00:43:28,588 --> 00:43:30,890 going to find a steady state solution with that one. 941 00:43:30,890 --> 00:43:34,717 And so I could time march until I retire. 942 00:43:34,717 --> 00:43:36,800 My computer's still burning CPU time all that time 943 00:43:36,800 --> 00:43:38,980 and it'll never get to a steady state solution. 944 00:43:38,980 --> 00:43:41,539 So you really have to know what the real situation is. 945 00:43:41,539 --> 00:43:43,580 But if you have one where it's going to converge, 946 00:43:43,580 --> 00:43:45,170 then this is a reasonable idea. 947 00:43:45,170 --> 00:43:48,900 Now let's talk about the time marching. 948 00:43:48,900 --> 00:43:53,436 The time marching, you have, say, dy, 949 00:43:53,436 --> 00:43:59,820 dt is equal to some f of y where this 950 00:43:59,820 --> 00:44:03,750 is the discretized spatial vergence. 951 00:44:03,750 --> 00:44:05,977 I took all the space dimensions and replaced them 952 00:44:05,977 --> 00:44:07,560 with the finite difference expressions 953 00:44:07,560 --> 00:44:08,820 for the derivatives. 954 00:44:08,820 --> 00:44:12,570 Or I did colocation or something to convert this 955 00:44:12,570 --> 00:44:15,912 into algebraic equations. 956 00:44:15,912 --> 00:44:18,120 So the right-hand side is now algebraic equations, no 957 00:44:18,120 --> 00:44:21,090 longer differentials. 958 00:44:21,090 --> 00:44:25,200 And now, I'm going to time march this and I just-- 959 00:44:25,200 --> 00:44:25,950 This is very long. 960 00:44:25,950 --> 00:44:28,470 It's a 100 million. 961 00:44:28,470 --> 00:44:30,720 So how are we going to time march that? 962 00:44:30,720 --> 00:44:33,690 Well there's kind of two schools of thought for this. 963 00:44:33,690 --> 00:44:38,090 One school of thought is, if I can use an explicit method, 964 00:44:38,090 --> 00:44:41,220 an explicit ODE solver, those are pretty nice. 965 00:44:41,220 --> 00:44:44,070 Because all I have to do is evaluate f and then 966 00:44:44,070 --> 00:44:47,610 I can march along and I never have to save any matrices. 967 00:44:47,610 --> 00:44:51,630 Because remember, look at those formulas, y, 968 00:44:51,630 --> 00:44:56,340 y new, by the explicit method, is 969 00:44:56,340 --> 00:45:00,320 equal to some function g of y old. 970 00:45:04,469 --> 00:45:05,760 And this is the update formula. 971 00:45:05,760 --> 00:45:06,760 It can be forward Euler. 972 00:45:06,760 --> 00:45:08,240 It could be RK 2. 973 00:45:08,240 --> 00:45:09,346 It could be whatever. 974 00:45:09,346 --> 00:45:10,970 But there's just some explicit formula. 975 00:45:10,970 --> 00:45:12,500 I plug in the old vector. 976 00:45:12,500 --> 00:45:13,730 I compute a new vector. 977 00:45:13,730 --> 00:45:18,050 I never had to do any inversions of matrices or anything. 978 00:45:18,050 --> 00:45:18,980 So that's good. 979 00:45:18,980 --> 00:45:21,530 And this is, indeed, the way that the Navier-Stokes 980 00:45:21,530 --> 00:45:25,310 equations are solved currently by the best possible solver. 981 00:45:25,310 --> 00:45:30,050 So the best solver in the world discretizes in space 982 00:45:30,050 --> 00:45:34,970 and then does ODE IVP using explicit method 983 00:45:34,970 --> 00:45:35,899 for marching it. 984 00:45:35,899 --> 00:45:37,440 They do problems that are pretty big. 985 00:45:37,440 --> 00:45:40,190 So they'll have like 100 million mesh points and then 100 986 00:45:40,190 --> 00:45:41,450 species. 987 00:45:41,450 --> 00:45:44,870 So they have 10 to the 10th unknowns. 988 00:45:44,870 --> 00:45:47,120 And at each time step, you're computing another vector 989 00:45:47,120 --> 00:45:49,310 that's 10 to the 10th long. 990 00:45:49,310 --> 00:45:51,270 So it's a pretty big problem. 991 00:45:51,270 --> 00:45:55,460 And they use, like, 40% of one of the national supercomputer 992 00:45:55,460 --> 00:45:58,310 centers will be running at one time for one job like this. 993 00:45:58,310 --> 00:45:59,701 But it works. 994 00:45:59,701 --> 00:46:01,700 You never have to store anything that's too big. 995 00:46:01,700 --> 00:46:04,070 The biggest thing you have to store are the vectors. 996 00:46:04,070 --> 00:46:05,030 Which, you're always going to have 997 00:46:05,030 --> 00:46:06,988 to store some object the size of your solution, 998 00:46:06,988 --> 00:46:09,230 right, if you're going to get your solution. 999 00:46:09,230 --> 00:46:10,880 So that's one possibility. 1000 00:46:10,880 --> 00:46:16,240 Now, this is really good if an explicit solver can solve it. 1001 00:46:16,240 --> 00:46:20,510 So that means it's just not stiff, your problem. 1002 00:46:20,510 --> 00:46:24,140 And also, if you want a time accurate solution. 1003 00:46:24,140 --> 00:46:26,810 So in that case, they actually want the time accurate 1004 00:46:26,810 --> 00:46:27,405 solution. 1005 00:46:27,405 --> 00:46:29,780 But suppose we're trying to really solve the steady state 1006 00:46:29,780 --> 00:46:31,547 problem over here. 1007 00:46:31,547 --> 00:46:33,630 Then we really don't care about the time accuracy. 1008 00:46:33,630 --> 00:46:36,720 We're not going to report anything about our time steps. 1009 00:46:36,720 --> 00:46:38,180 In the end, we're only going to report our steady state 1010 00:46:38,180 --> 00:46:38,679 solution. 1011 00:46:38,679 --> 00:46:40,044 That's all we cared about. 1012 00:46:40,044 --> 00:46:41,460 So in that case, there's no reason 1013 00:46:41,460 --> 00:46:44,580 to try to have a really high accuracy, explicit formula 1014 00:46:44,580 --> 00:46:47,250 and try to make sure we all our time steps are exactly right. 1015 00:46:47,250 --> 00:46:48,780 We're just going to throw away all those time points anyway. 1016 00:46:48,780 --> 00:46:51,360 We'll just keep the final one as our initial guess 1017 00:46:51,360 --> 00:46:54,630 for a Newton-Raphson solve to find the steady state. 1018 00:46:54,630 --> 00:46:57,840 So in those cases, you might instead want to use-- 1019 00:46:57,840 --> 00:46:59,850 So this is time accurate. 1020 00:47:05,312 --> 00:47:06,270 There's another method. 1021 00:47:06,270 --> 00:47:10,260 I don't know if I should call it time inaccurate. 1022 00:47:10,260 --> 00:47:14,250 If you don't care, if the solution is really, 1023 00:47:14,250 --> 00:47:16,320 what the solution yft is, all you 1024 00:47:16,320 --> 00:47:19,770 care about is where you get to, then you can do other methods. 1025 00:47:19,770 --> 00:47:22,230 And one of them is the backward Euler. 1026 00:47:22,230 --> 00:47:30,380 So you can say, well, y new is equal to y old 1027 00:47:30,380 --> 00:47:35,660 plus delta t times f of y new. 1028 00:47:35,660 --> 00:47:39,150 Now, the advantage of this is, I can make delta t pretty large. 1029 00:47:39,150 --> 00:47:41,250 This is guaranteed stable as long 1030 00:47:41,250 --> 00:47:42,640 as the ODE system is stable. 1031 00:47:42,640 --> 00:47:43,890 Remember? 1032 00:47:43,890 --> 00:47:49,110 And so, this will go to this, go to the solution pretty well, 1033 00:47:49,110 --> 00:47:50,745 to the steady state solution. 1034 00:47:50,745 --> 00:47:52,620 It won't go there in the right amount of time 1035 00:47:52,620 --> 00:47:54,760 because this is a very low order of formula. 1036 00:47:54,760 --> 00:47:57,190 It's not really accurate. 1037 00:47:57,190 --> 00:47:59,620 But it will end up to the real steady state solution. 1038 00:47:59,620 --> 00:48:01,220 This is what's actually used a lot. 1039 00:48:01,220 --> 00:48:04,210 However, the problem with this is, this is now 1040 00:48:04,210 --> 00:48:07,530 an implicit equation. 1041 00:48:07,530 --> 00:48:09,810 So we may have to solve it with some method 1042 00:48:09,810 --> 00:48:11,904 like Newton-Raphson. 1043 00:48:11,904 --> 00:48:13,570 But we got into this because we couldn't 1044 00:48:13,570 --> 00:48:15,569 solve our Newton-Raphson steps because we didn't 1045 00:48:15,569 --> 00:48:17,230 have a good initial guess. 1046 00:48:17,230 --> 00:48:19,890 So then the question is, can we find a good initial guess 1047 00:48:19,890 --> 00:48:21,839 for this problem? 1048 00:48:21,839 --> 00:48:23,380 And in this case, you can because you 1049 00:48:23,380 --> 00:48:27,970 can use the explicit formulas to again provide the initial guess 1050 00:48:27,970 --> 00:48:30,110 for this implicit formula. 1051 00:48:30,110 --> 00:48:33,220 So that's what they do, also. 1052 00:48:33,220 --> 00:48:36,204 And so now, if we're doing an iterative procedure in order 1053 00:48:36,204 --> 00:48:38,620 to get the initial guess for our other iterative procedure 1054 00:48:38,620 --> 00:48:40,090 ever there, which we're going to break up 1055 00:48:40,090 --> 00:48:42,715 into another iterative procedure inside it to solve every step. 1056 00:48:45,910 --> 00:48:47,860 But this is currently what people do. 1057 00:48:47,860 --> 00:48:49,360 So you guys are smart. 1058 00:48:49,360 --> 00:48:51,610 Maybe you can figure a better way to do it. 1059 00:48:51,610 --> 00:48:54,360 The world is waiting. 1060 00:48:54,360 --> 00:48:55,360 I think I'll stop there. 1061 00:48:57,880 --> 00:48:59,890 By the way, actually just names-- 1062 00:48:59,890 --> 00:49:02,980 This kind of thing is called method of lines. 1063 00:49:07,350 --> 00:49:09,630 And if so, if anybody ever says I solved something 1064 00:49:09,630 --> 00:49:12,060 by method of lines, what they meant was, 1065 00:49:12,060 --> 00:49:14,389 they discretized in some of the coordinates 1066 00:49:14,389 --> 00:49:16,680 and they kept one of the coordinates as a differential. 1067 00:49:16,680 --> 00:49:20,660 And then they solved it with an ODE solver at the end. 1068 00:49:20,660 --> 00:49:23,630 And in certain systems like this, it's a smart thing to do. 1069 00:49:23,630 --> 00:49:26,240 You need to be kind of lucky that the set up is right, 1070 00:49:26,240 --> 00:49:27,900 so you can use it. 1071 00:49:27,900 --> 00:49:29,960 But if you can, it can be pretty good. 1072 00:49:29,960 --> 00:49:31,290 All right. 1073 00:49:31,290 --> 00:49:33,190 See you on Friday.