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:19,310 at ocw.mit.edu. 8 00:00:23,490 --> 00:00:26,880 WILLIAM GREEN: So welcome to our special class with half the 9 00:00:26,880 --> 00:00:29,550 people gone for tonight's hearing. 10 00:00:29,550 --> 00:00:33,850 You guys are the chosen few, so I'm pleased you're here. 11 00:00:33,850 --> 00:00:38,814 AUDIENCE: [INTERPOSING VOICES] 12 00:00:38,814 --> 00:00:40,980 WILLIAM GREEN: So I'm just going to say a little bit 13 00:00:40,980 --> 00:00:43,770 about models versus data and then 14 00:00:43,770 --> 00:00:47,913 I'm going to go back to undervalue problems and talk 15 00:00:47,913 --> 00:00:49,996 a little bit about how you handle really big ones. 16 00:00:49,996 --> 00:00:50,928 Yeah? 17 00:00:50,928 --> 00:00:52,803 AUDIENCE: Are these going to be [INAUDIBLE].. 18 00:00:52,803 --> 00:00:54,990 WILLIAM GREEN: Yeah. 19 00:00:54,990 --> 00:00:55,940 Yes, sir. 20 00:00:55,940 --> 00:00:56,890 AUDIENCE: All right. 21 00:00:56,890 --> 00:00:57,410 WILLIAM GREEN: And in fact, these 22 00:00:57,410 --> 00:00:58,640 are even going to be on the video, 23 00:00:58,640 --> 00:00:59,780 because the guy just this fixed the video, 24 00:00:59,780 --> 00:01:00,738 so it will be recorded. 25 00:01:03,230 --> 00:01:05,534 All right. 26 00:01:05,534 --> 00:01:06,200 Yeah, I'm sorry. 27 00:01:06,200 --> 00:01:07,670 You said the video was broken last time 28 00:01:07,670 --> 00:01:08,920 so that recording didn't work. 29 00:01:11,870 --> 00:01:13,110 All right. 30 00:01:13,110 --> 00:01:17,794 So when we're talking about models versus data 31 00:01:17,794 --> 00:01:20,580 a really critical question is whether the model 32 00:01:20,580 --> 00:01:23,260 is really consistent with the data or not. 33 00:01:23,260 --> 00:01:29,280 And how do we know this is really true? 34 00:01:29,280 --> 00:01:35,890 So we talked about how we can figure out 35 00:01:35,890 --> 00:01:39,860 some sort of chi-squared maximum that we would tolerate. 36 00:01:39,860 --> 00:01:41,410 So the maximum that we tolerate, we 37 00:01:41,410 --> 00:01:46,180 call it Q, capital Q in a lecture before. 38 00:01:46,180 --> 00:01:51,690 And so if we can find a chi-squared of theta that's 39 00:01:51,690 --> 00:01:55,320 less than this, then it's at least plausible 40 00:01:55,320 --> 00:01:56,740 that our model is true. 41 00:01:56,740 --> 00:01:57,240 Right? 42 00:01:57,240 --> 00:02:03,290 Or, alternatively, if all chi-squared thetas 43 00:02:03,290 --> 00:02:08,030 for any value of thetas are greater than chi-squared max, 44 00:02:08,030 --> 00:02:10,190 then we would say the model is false, 45 00:02:10,190 --> 00:02:14,310 or something else is really wrong with our experiment. 46 00:02:14,310 --> 00:02:15,291 Right? 47 00:02:15,291 --> 00:02:16,790 That's what we say when we're saying 48 00:02:16,790 --> 00:02:20,990 we have a maximum chi-squared that we would tolerate. 49 00:02:20,990 --> 00:02:24,200 So a critical question with us is 50 00:02:24,200 --> 00:02:27,620 whether we really found the very best fit 51 00:02:27,620 --> 00:02:33,800 value of the parameters. 52 00:02:33,800 --> 00:02:37,910 In some simple problems you might have linear dependence 53 00:02:37,910 --> 00:02:41,930 on the parameters, and then there's only one 54 00:02:41,930 --> 00:02:43,880 best fit usually, as Long as you have 55 00:02:43,880 --> 00:02:47,330 a non-singular set of equations and everything 56 00:02:47,330 --> 00:02:49,130 we find there's no trouble. 57 00:02:49,130 --> 00:02:52,070 However, more commonly you have equations 58 00:02:52,070 --> 00:02:54,560 which are singular or close to singular, 59 00:02:54,560 --> 00:02:58,460 and also you normally have non-linear dependence 60 00:02:58,460 --> 00:03:01,710 between the model predictions and the parameters. 61 00:03:01,710 --> 00:03:03,500 So, for example, I do kinetics. 62 00:03:03,500 --> 00:03:06,920 We always care about the ray coefficients, the K's and they 63 00:03:06,920 --> 00:03:09,070 always come in inside a differential equation. 64 00:03:09,070 --> 00:03:09,569 Right? 65 00:03:09,569 --> 00:03:12,005 You always have like d concentrations 66 00:03:12,005 --> 00:03:15,080 dt is equal to some terms, and some of them 67 00:03:15,080 --> 00:03:21,100 are, you know, k concentration I times concentration j. 68 00:03:21,100 --> 00:03:22,250 whatever one it is, ij. 69 00:03:22,250 --> 00:03:22,750 Right? 70 00:03:22,750 --> 00:03:25,350 You have a lot of like terms like this by molecular reaction 71 00:03:25,350 --> 00:03:26,670 terms. 72 00:03:26,670 --> 00:03:30,030 And these K's, some of them I don't know. 73 00:03:30,030 --> 00:03:33,230 And so these would be one of my thetas. 74 00:03:33,230 --> 00:03:34,770 So this would be, you know, theta 17 75 00:03:34,770 --> 00:03:38,350 would be one of these K's that I want to determine. 76 00:03:38,350 --> 00:03:41,130 And when I have a differential equation like this, 77 00:03:41,130 --> 00:03:43,690 when I integrate it, the K is always 78 00:03:43,690 --> 00:03:47,280 going to come in non-linearly in the solution. 79 00:03:47,280 --> 00:03:47,780 Right? 80 00:03:47,780 --> 00:03:49,900 So when I integrate this solution, 81 00:03:49,900 --> 00:03:53,580 even in the simplest case you get things like c of t 82 00:03:53,580 --> 00:03:57,480 is equal to ae to the negative kt. 83 00:03:57,480 --> 00:03:59,800 That might be a really super simple case 84 00:03:59,800 --> 00:04:01,720 if you have a linear kinetics. 85 00:04:01,720 --> 00:04:04,570 The K is still coming in non-linearly 86 00:04:04,570 --> 00:04:05,620 to the predictions. 87 00:04:05,620 --> 00:04:06,640 Right? 88 00:04:06,640 --> 00:04:08,440 And if you have anything complicated 89 00:04:08,440 --> 00:04:11,530 you can't even write it down in a closed form expression, 90 00:04:11,530 --> 00:04:17,200 and you have to integrate with oed45 or ode15s 91 00:04:17,200 --> 00:04:18,940 and so the relationship between K 92 00:04:18,940 --> 00:04:21,820 and the concentrations that you measure 93 00:04:21,820 --> 00:04:23,467 can be really complicated. 94 00:04:23,467 --> 00:04:24,550 So it's always non-linear. 95 00:04:24,550 --> 00:04:29,620 In my world all the models are non-linear. 96 00:04:29,620 --> 00:04:31,930 And I think there's many other situations like this. 97 00:04:31,930 --> 00:04:38,287 So if you try to measure a diffusivity, 98 00:04:38,287 --> 00:04:40,120 it might be possible to set up an experiment 99 00:04:40,120 --> 00:04:41,670 where you'd expect a linear dependency to something 100 00:04:41,670 --> 00:04:44,280 you'd measure in a diffusivity, but most likely not. 101 00:04:44,280 --> 00:04:47,470 So most likely it will come in non-linearly somehow 102 00:04:47,470 --> 00:04:49,787 and you'd have some model for it, 103 00:04:49,787 --> 00:04:51,370 and you'd have a non-linear dependence 104 00:04:51,370 --> 00:04:54,640 between that number and the numbers you would measure. 105 00:04:54,640 --> 00:04:55,420 Right? 106 00:04:55,420 --> 00:04:58,130 And almost all the properties are like that. 107 00:04:58,130 --> 00:05:03,437 So because you have a non-linear dependence, 108 00:05:03,437 --> 00:05:05,270 then you have to worry about how many minima 109 00:05:05,270 --> 00:05:06,520 do you have in your surface? 110 00:05:06,520 --> 00:05:10,300 So chi-squared of theta is a non-linear function 111 00:05:10,300 --> 00:05:11,920 of a lot of variables. 112 00:05:11,920 --> 00:05:14,350 And typically it won't just have one minimum, 113 00:05:14,350 --> 00:05:16,530 you'll have a whole lot of local minima. 114 00:05:16,530 --> 00:05:18,190 Now one of those is the global minimum. 115 00:05:18,190 --> 00:05:19,648 And that's the one you want, that's 116 00:05:19,648 --> 00:05:21,600 the very best possible fit. 117 00:05:21,600 --> 00:05:26,070 But you normally won't know how many minima it has 118 00:05:26,070 --> 00:05:28,170 and so no matter how many minima you 119 00:05:28,170 --> 00:05:31,170 find you don't know if you're done, 120 00:05:31,170 --> 00:05:33,900 and you don't know if you really found the best fit. 121 00:05:33,900 --> 00:05:34,470 OK? 122 00:05:34,470 --> 00:05:37,260 So just because all the ones you found so far bigger 123 00:05:37,260 --> 00:05:39,150 than chi-squared max you might not 124 00:05:39,150 --> 00:05:41,015 be ready to write that article in Nature 125 00:05:41,015 --> 00:05:42,390 to claim that the laws of physics 126 00:05:42,390 --> 00:05:44,389 are incorrect because you're not sure you really 127 00:05:44,389 --> 00:05:45,810 found the very best fit. 128 00:05:45,810 --> 00:05:46,310 OK? 129 00:05:46,310 --> 00:05:49,240 So that's one of the many problems 130 00:05:49,240 --> 00:05:51,380 of trying to figure out if a model of the data 131 00:05:51,380 --> 00:05:53,690 are really consistent. 132 00:05:53,690 --> 00:05:56,530 Now that particular one, you can overcome 133 00:05:56,530 --> 00:06:00,880 if you can find a guaranteed way to find the global optimum, 134 00:06:00,880 --> 00:06:03,070 to find the very best possible that there is. 135 00:06:03,070 --> 00:06:05,320 And some people spend their life working on this, 136 00:06:05,320 --> 00:06:07,120 and one of those people is professor Barton 137 00:06:07,120 --> 00:06:08,290 in this department. 138 00:06:08,290 --> 00:06:10,220 And so he teaches a class on how to do this, 139 00:06:10,220 --> 00:06:12,640 and there are actually methods to do it. 140 00:06:12,640 --> 00:06:15,130 So if you want to see the application of this 141 00:06:15,130 --> 00:06:18,330 to kinetic models there's a paper 142 00:06:18,330 --> 00:06:19,760 by me and professor Barton. 143 00:06:19,760 --> 00:06:22,420 And the lead author is one of his former students, 144 00:06:22,420 --> 00:06:24,200 Adam Singer and J. Phys. 145 00:06:24,200 --> 00:06:29,387 Chem. and it just shows how you can use global optimization 146 00:06:29,387 --> 00:06:31,220 to guarantee you have the very best possible 147 00:06:31,220 --> 00:06:34,430 via the parameters, and then they'll show you what happens. 148 00:06:34,430 --> 00:06:38,460 So this is a case where the experimental data was measured 149 00:06:38,460 --> 00:06:42,150 by one of my students, James Taylor, and that's 150 00:06:42,150 --> 00:06:45,050 the thing with the error bars, and he did a lot of hard work 151 00:06:45,050 --> 00:06:46,800 to make sure he had pretty good error bars 152 00:06:46,800 --> 00:06:48,850 and he really knew what they were. 153 00:06:48,850 --> 00:06:50,020 And so he had those points. 154 00:06:50,020 --> 00:06:52,002 So he measured a lot of points. 155 00:06:52,002 --> 00:06:54,080 This measuring a lot of points is a key thing. 156 00:06:54,080 --> 00:06:57,572 So it means his degrees of freedom is very large. 157 00:06:57,572 --> 00:06:59,280 Because he uses many, many different time 158 00:06:59,280 --> 00:07:01,380 points along the way there and he 159 00:07:01,380 --> 00:07:03,060 has error bars on all of them. 160 00:07:03,060 --> 00:07:06,410 And so because the degrees of freedom are very large 161 00:07:06,410 --> 00:07:11,265 it actually makes it a very good constraint on the model. 162 00:07:11,265 --> 00:07:12,630 It tested very thoroughly. 163 00:07:12,630 --> 00:07:14,440 Where if you only measured one point, 164 00:07:14,440 --> 00:07:17,590 then that's not a very good constraint. 165 00:07:17,590 --> 00:07:20,820 Now what James had done, and we published a paper in 2004 where 166 00:07:20,820 --> 00:07:24,330 he did it, where he said, you know, when I do my fits 167 00:07:24,330 --> 00:07:25,770 this is what I get. 168 00:07:25,770 --> 00:07:29,230 I get these kind of purple and green things, 169 00:07:29,230 --> 00:07:32,770 and therefore they're all outside the error bars, 170 00:07:32,770 --> 00:07:35,500 and therefore I publish in the Journal of Physical Chemistry, 171 00:07:35,500 --> 00:07:39,590 I claim this model is wrong based on my experimental data. 172 00:07:39,590 --> 00:07:40,090 OK? 173 00:07:40,090 --> 00:07:43,270 So we published that in 2004, it was accepted by the referees. 174 00:07:43,270 --> 00:07:45,010 It's in the literature, you can read it. 175 00:07:45,010 --> 00:07:49,180 Don't believe it though, because Adam singer went back and did 176 00:07:49,180 --> 00:07:51,640 global optimisation using the same model 177 00:07:51,640 --> 00:07:56,220 and his red line goes right through all the points really 178 00:07:56,220 --> 00:07:56,720 well. 179 00:07:56,720 --> 00:07:57,530 Do you see that? 180 00:07:57,530 --> 00:07:59,260 You can see the red line there. 181 00:07:59,260 --> 00:08:01,030 See this red curve? 182 00:08:01,030 --> 00:08:02,830 That's the true best fit. 183 00:08:02,830 --> 00:08:05,080 And so all the fits that James had found 184 00:08:05,080 --> 00:08:08,840 were all local minima in the chi-squared squared surface. 185 00:08:08,840 --> 00:08:10,600 And when he actually found the best fit 186 00:08:10,600 --> 00:08:12,436 with Adam's fancy method, then it 187 00:08:12,436 --> 00:08:13,810 turns out that actually the model 188 00:08:13,810 --> 00:08:16,270 matched the data perfectly. 189 00:08:16,270 --> 00:08:18,160 So that's definitely true. 190 00:08:18,160 --> 00:08:20,680 So what we published the first time 191 00:08:20,680 --> 00:08:22,742 was completely wrong, because we said 192 00:08:22,742 --> 00:08:23,950 the model disproved the data. 193 00:08:23,950 --> 00:08:25,666 But in fact the model proved the data. 194 00:08:25,666 --> 00:08:27,040 I mean the data proved the model, 195 00:08:27,040 --> 00:08:28,831 or was completely consistent with the model 196 00:08:28,831 --> 00:08:30,760 in every way in great detail. 197 00:08:30,760 --> 00:08:33,390 So I was a little embarrassed. 198 00:08:33,390 --> 00:08:40,559 So then James, also in the same paper, had found-- 199 00:08:40,559 --> 00:08:42,640 and they measured another condition 200 00:08:42,640 --> 00:08:45,880 in a slightly different case, and actually 201 00:08:45,880 --> 00:08:47,740 was looking at a different solvent. 202 00:08:47,740 --> 00:08:53,530 And in this case he found this green curve as his best fit, 203 00:08:53,530 --> 00:08:58,284 and it looked pretty good by eye, 204 00:08:58,284 --> 00:08:59,950 though you can see it sort of skittering 205 00:08:59,950 --> 00:09:03,280 on the tops of the error bar bands. 206 00:09:03,280 --> 00:09:06,450 And so this is one where he said, you know, 207 00:09:06,450 --> 00:09:10,000 we do the statistical test, there's like a 5% chance 208 00:09:10,000 --> 00:09:12,190 that you would measure data like he 209 00:09:12,190 --> 00:09:17,440 measured if the green model was the truth. 210 00:09:17,440 --> 00:09:21,610 And so we said some wishy-washy words in the paper, 211 00:09:21,610 --> 00:09:23,210 basically because we weren't really 212 00:09:23,210 --> 00:09:27,160 if it was consistent or not consistent with the data. 213 00:09:27,160 --> 00:09:29,722 But in that case Adam used his fancy global optimization 214 00:09:29,722 --> 00:09:31,930 techniques and found that James' fit was actually not 215 00:09:31,930 --> 00:09:34,180 that far from the global optimum, 216 00:09:34,180 --> 00:09:37,159 and he got the red curve there, marked by the red arrow, 217 00:09:37,159 --> 00:09:38,950 and that one looks pretty good, it actually 218 00:09:38,950 --> 00:09:40,960 goes through the data. 219 00:09:40,960 --> 00:09:43,570 But it doesn't really do it right, 220 00:09:43,570 --> 00:09:47,770 and found that the confidence intervals are only 16%. 221 00:09:47,770 --> 00:09:50,800 So it's only a 16% chance, one chance out of six, 222 00:09:50,800 --> 00:09:53,350 that James would have measured the data he measured 223 00:09:53,350 --> 00:09:55,330 if that model was true. 224 00:09:55,330 --> 00:09:57,724 And now we're sure that we have the global best fit, 225 00:09:57,724 --> 00:09:59,140 so now we're, like, thinking well, 226 00:09:59,140 --> 00:10:02,350 you know, maybe we should say that the model is not 227 00:10:02,350 --> 00:10:04,390 true for this case. 228 00:10:04,390 --> 00:10:07,570 Because it's an 84% chance it's not true anyway. 229 00:10:07,570 --> 00:10:10,670 So actually we got it wrong both ways. 230 00:10:10,670 --> 00:10:12,740 So in the first case we said that our data 231 00:10:12,740 --> 00:10:14,400 disproved the model. 232 00:10:14,400 --> 00:10:20,229 And in the second case, we said that our model was pretty good 233 00:10:20,229 --> 00:10:22,270 with the data, but actually they were both wrong. 234 00:10:22,270 --> 00:10:26,267 So one was wrong that actually-- 235 00:10:26,267 --> 00:10:28,600 the data did prove the models correct for the first case 236 00:10:28,600 --> 00:10:34,330 and it actually looks a little fishy for the second case. 237 00:10:34,330 --> 00:10:37,840 So anyway this is just a concrete example 238 00:10:37,840 --> 00:10:41,444 to make you be a little bit more worried about what you do when 239 00:10:41,444 --> 00:10:42,610 you do these fitting things. 240 00:10:45,250 --> 00:10:49,810 Now there's a variety of ways to try to do the least course 241 00:10:49,810 --> 00:10:52,600 fitting to try to find optima. 242 00:10:52,600 --> 00:10:56,639 There was a guy named Carr from University of Minnesota. 243 00:10:56,639 --> 00:10:58,180 And he published this [INAUDIBLE] way 244 00:10:58,180 --> 00:11:00,910 to just try a whole lot of different initial guesses, 245 00:11:00,910 --> 00:11:05,364 and that usually works to find the global minimum. 246 00:11:05,364 --> 00:11:07,780 But if you really want to have a guaranteed global minimum 247 00:11:07,780 --> 00:11:09,488 that you're absolutely sure is definitely 248 00:11:09,488 --> 00:11:11,530 is the global minimum, then you really 249 00:11:11,530 --> 00:11:13,922 need to use global optimization methods 250 00:11:13,922 --> 00:11:15,880 and you really should talk to professor Barton. 251 00:11:15,880 --> 00:11:18,010 And for these kinetic models he actually has 252 00:11:18,010 --> 00:11:21,190 a distributed code called GDOC that will just 253 00:11:21,190 --> 00:11:23,440 find the global optimization as long as you don't have 254 00:11:23,440 --> 00:11:25,450 too many kinetic parameters. 255 00:11:25,450 --> 00:11:27,220 I think if you go to about 6-- 256 00:11:27,220 --> 00:11:29,980 around 6 adjustable parameters it can handle. 257 00:11:29,980 --> 00:11:33,557 If you have more than that then usually typically not. 258 00:11:33,557 --> 00:11:35,723 Because the problem with global optimisation methods 259 00:11:35,723 --> 00:11:38,980 is that they scale exponentially with the number 260 00:11:38,980 --> 00:11:42,370 of adjustable parameters, or the number of things 261 00:11:42,370 --> 00:11:44,385 you're trying to optimize. 262 00:11:44,385 --> 00:11:46,000 Now on the other hand, if you have 263 00:11:46,000 --> 00:11:47,620 a model that you need to use more 264 00:11:47,620 --> 00:11:50,362 than six adjustable parameters to fit your data, 265 00:11:50,362 --> 00:11:51,820 then you might want to try to think 266 00:11:51,820 --> 00:11:53,050 of a different experiment. 267 00:11:53,050 --> 00:11:56,830 Because you know the famous saying is like, you know, 268 00:11:56,830 --> 00:12:00,910 give me N parameters I can fit an elephant. 269 00:12:00,910 --> 00:12:03,220 So if you have more than six parameters 270 00:12:03,220 --> 00:12:05,714 and you're trying to fit some data set, 271 00:12:05,714 --> 00:12:07,630 probably if you do it right you'll find a fit, 272 00:12:07,630 --> 00:12:09,171 but it may not really prove anything. 273 00:12:11,410 --> 00:12:13,270 So from that point of view professor 274 00:12:13,270 --> 00:12:15,759 Barton's code, GDOC, is pretty convenient for doing 275 00:12:15,759 --> 00:12:17,050 these non-linear optimizations. 276 00:12:17,050 --> 00:12:20,500 So that's a code for non-linear optimization with differential 277 00:12:20,500 --> 00:12:22,120 equations embedded. 278 00:12:22,120 --> 00:12:24,060 And they can be stiff and they can be large, 279 00:12:24,060 --> 00:12:24,870 it doesn't matter. 280 00:12:24,870 --> 00:12:26,536 What really matters it's just the number 281 00:12:26,536 --> 00:12:31,550 of adjustable parameters, that's the cost. 282 00:12:31,550 --> 00:12:35,570 All right now there was a really good question, 283 00:12:35,570 --> 00:12:38,770 that you asked I think, in class last time, maybe. 284 00:12:38,770 --> 00:12:39,710 Maybe the time before. 285 00:12:39,710 --> 00:12:43,010 And it was what happens-- 286 00:12:43,010 --> 00:12:45,340 so let's look at different case you get. 287 00:12:45,340 --> 00:12:50,630 So one case you get, you know, if chi-squared square, 288 00:12:50,630 --> 00:12:55,260 theta best, your best fit parameters is a lot bigger, 289 00:12:55,260 --> 00:12:59,620 think chi-squared max, then you're ready to say the model 290 00:12:59,620 --> 00:13:00,608 disproves it. 291 00:13:05,760 --> 00:13:07,760 So that one you're OK. 292 00:13:07,760 --> 00:13:14,350 If chi-squared of your best fit is less than, 293 00:13:14,350 --> 00:13:17,341 or significantly less than chi-squared max-- 294 00:13:17,341 --> 00:13:21,060 it has be much less than it, really less, not equal. 295 00:13:23,610 --> 00:13:25,320 --then you're really happy to say 296 00:13:25,320 --> 00:13:27,510 that the models are consistent. 297 00:13:27,510 --> 00:13:29,685 The model and data are consistent to each other. 298 00:13:33,405 --> 00:13:39,080 And, in this case, what you would normally do next 299 00:13:39,080 --> 00:13:48,700 is determine the confidence intervals on the thetas 300 00:13:48,700 --> 00:13:50,470 that you adjusted. 301 00:13:50,470 --> 00:13:52,510 Right? 302 00:13:52,510 --> 00:13:55,030 And if you read the chapter in Numerical Recipes, 303 00:13:55,030 --> 00:13:56,940 for example, some of the notes online, 304 00:13:56,940 --> 00:13:58,540 it could tell you how to do this. 305 00:13:58,540 --> 00:14:00,060 Right? 306 00:14:00,060 --> 00:14:01,715 And so that was OK? 307 00:14:01,715 --> 00:14:02,890 It's OK. 308 00:14:02,890 --> 00:14:05,140 So the case that's really problematic 309 00:14:05,140 --> 00:14:09,360 is you find chi-squared theta best, 310 00:14:09,360 --> 00:14:10,830 and it's less than chi-squared max, 311 00:14:10,830 --> 00:14:12,890 so you think the model is consistent with data, 312 00:14:12,890 --> 00:14:17,760 but it's pretty darn close to theta-squared max. 313 00:14:17,760 --> 00:14:22,980 And so then you have to figure out, well, the weird aspect 314 00:14:22,980 --> 00:14:27,779 of this is I now have-- 315 00:14:27,779 --> 00:14:29,820 here's my famous plot of theta one and theta two. 316 00:14:29,820 --> 00:14:31,599 So I have two parameters I'm adjusting. 317 00:14:31,599 --> 00:14:33,390 If they want theta two, here's my best fit. 318 00:14:33,390 --> 00:14:34,312 My best fit. 319 00:14:37,540 --> 00:14:41,590 If chi-squared max-- if the curve 320 00:14:41,590 --> 00:14:45,916 of constant chi-squared max is something like this, 321 00:14:45,916 --> 00:14:52,650 this is chi-squared of theta equals chi-squared max. 322 00:14:52,650 --> 00:14:55,444 If it's like that then I have some big region 323 00:14:55,444 --> 00:14:56,860 of a lot of different theta values 324 00:14:56,860 --> 00:14:58,492 that give pretty good fits, and I'm 325 00:14:58,492 --> 00:14:59,950 happy to draw confidence intervals, 326 00:14:59,950 --> 00:15:01,845 like draw lines here and here, for example, and say 327 00:15:01,845 --> 00:15:03,760 that's the confidence interval on theta one, 328 00:15:03,760 --> 00:15:04,870 and I'm all right. 329 00:15:04,870 --> 00:15:05,740 OK? 330 00:15:05,740 --> 00:15:08,247 And it's sloped here so I might want 331 00:15:08,247 --> 00:15:10,330 to tell people about the co-variance between theta 332 00:15:10,330 --> 00:15:11,830 one and theta two. 333 00:15:11,830 --> 00:15:12,350 All right? 334 00:15:12,350 --> 00:15:14,620 But I know what to do, it's fine. 335 00:15:14,620 --> 00:15:17,817 Now suppose this is the case that's like this. 336 00:15:17,817 --> 00:15:19,900 Where chi-squared theta best is significantly less 337 00:15:19,900 --> 00:15:21,550 than chi-squared max. 338 00:15:21,550 --> 00:15:25,620 Now suppose instead chi-squared best 339 00:15:25,620 --> 00:15:27,080 was very close chi-squared max. 340 00:15:27,080 --> 00:15:29,416 So then here is chi-squared max. 341 00:15:29,416 --> 00:15:33,360 This is chi-squared of theta equals chi-squared max, 342 00:15:33,360 --> 00:15:37,497 this is the good case. 343 00:15:37,497 --> 00:15:38,330 Here's the bad case. 344 00:15:41,740 --> 00:15:44,860 So in the bad case, there is some region of thetas 345 00:15:44,860 --> 00:15:47,200 that the model is consistent with the data, 346 00:15:47,200 --> 00:15:49,720 but it's very tiny. 347 00:15:49,720 --> 00:15:51,920 There's a little tiny range of theta values 348 00:15:51,920 --> 00:15:55,130 that would make the models consistent with the data. 349 00:15:55,130 --> 00:15:57,800 Now what's weird about this is that my best fit, it's not 350 00:15:57,800 --> 00:15:59,180 very good. 351 00:15:59,180 --> 00:16:01,470 The chi-squared is almost up to the maximum 352 00:16:01,470 --> 00:16:02,725 chi-squared I would tolerate. 353 00:16:02,725 --> 00:16:04,850 So it means my deviations between my model and data 354 00:16:04,850 --> 00:16:05,490 are pretty bad. 355 00:16:05,490 --> 00:16:07,380 Like this plot actually. 356 00:16:07,380 --> 00:16:08,690 OK? 357 00:16:08,690 --> 00:16:11,090 And so I'm, like, you know, it could-- it's possible. 358 00:16:11,090 --> 00:16:13,990 I was a little bit unlucky and my data set 359 00:16:13,990 --> 00:16:18,070 is pretty far off but still within the possible range. 360 00:16:18,070 --> 00:16:20,940 And so maybe the model data are consistent. 361 00:16:20,940 --> 00:16:24,340 But then when I try to compute the confidence intervals 362 00:16:24,340 --> 00:16:27,550 on my parameters, there's only a really small range 363 00:16:27,550 --> 00:16:30,520 of parameters that are going to give good fits here. 364 00:16:30,520 --> 00:16:32,620 Because the best possible fit is this red line. 365 00:16:32,620 --> 00:16:35,050 There's no choice of the parameters 366 00:16:35,050 --> 00:16:39,690 that's going to go through all the circles with this model. 367 00:16:39,690 --> 00:16:41,990 So now I'm like, what's going on here? 368 00:16:44,937 --> 00:16:46,520 So if I went ahead and do the same way 369 00:16:46,520 --> 00:16:50,190 I did before, I draw my confidence intervals like this 370 00:16:50,190 --> 00:16:51,390 and say, OK. 371 00:16:51,390 --> 00:16:54,180 Say, wow, I've done an amazingly great job 372 00:16:54,180 --> 00:16:55,940 of determining theta one. 373 00:16:55,940 --> 00:16:58,486 I know theta one very, very, very precisely. 374 00:16:58,486 --> 00:16:59,610 And the same for theta two. 375 00:16:59,610 --> 00:17:01,526 Here's theta two, it's got to be in that range 376 00:17:01,526 --> 00:17:04,574 and I've really done a great job. 377 00:17:04,574 --> 00:17:06,270 So I'm awesome. 378 00:17:06,270 --> 00:17:08,640 But I'm only awesome because my original best 379 00:17:08,640 --> 00:17:11,515 fit wasn't very good. 380 00:17:11,515 --> 00:17:13,890 So this might make you feel a little sick in your stomach 381 00:17:13,890 --> 00:17:15,764 if you're getting ready to publish in nature. 382 00:17:15,764 --> 00:17:17,089 Right? 383 00:17:17,089 --> 00:17:21,119 And then if the reviewers are on their game, 384 00:17:21,119 --> 00:17:22,829 you might get something really scathing 385 00:17:22,829 --> 00:17:25,569 back from the reviewers. 386 00:17:25,569 --> 00:17:30,771 So there's kind of two ways to play it at this point. 387 00:17:30,771 --> 00:17:34,410 One way is, you are absolutely sure that your equations are 388 00:17:34,410 --> 00:17:35,230 correct. 389 00:17:35,230 --> 00:17:37,140 You are Sir Isaac Newton, you have just 390 00:17:37,140 --> 00:17:39,540 written the law of gravitation. 391 00:17:39,540 --> 00:17:42,150 You are measuring the orbits of the moons of Jupiter, 392 00:17:42,150 --> 00:17:44,970 and you are sure that your equations are right. 393 00:17:44,970 --> 00:17:49,770 And so, you know, my problem is my antiquated telescopes 394 00:17:49,770 --> 00:17:52,050 of 1600 are not that great, and so I 395 00:17:52,050 --> 00:17:54,300 think I have maybe-- maybe I misestimated my error 396 00:17:54,300 --> 00:17:55,971 bars on my measurements. 397 00:17:55,971 --> 00:17:57,720 Maybe I have some aberration or something, 398 00:17:57,720 --> 00:17:58,960 the angles are off, right? 399 00:17:58,960 --> 00:18:01,930 My twisted telescope looks as if the star moves or something, 400 00:18:01,930 --> 00:18:03,400 right? 401 00:18:03,400 --> 00:18:03,900 OK? 402 00:18:03,900 --> 00:18:08,040 So that's one possible way to go, 403 00:18:08,040 --> 00:18:11,650 say I really believe the model, and so I 404 00:18:11,650 --> 00:18:14,140 believe this is the truth, and I really 405 00:18:14,140 --> 00:18:15,980 believe my chi-squared max to be what it is, 406 00:18:15,980 --> 00:18:18,670 and I double checked the error bars and I think it's true. 407 00:18:18,670 --> 00:18:20,415 So I think have done a great job and I should publish it 408 00:18:20,415 --> 00:18:21,385 in Nature because I have determined 409 00:18:21,385 --> 00:18:23,650 these parameters with 17 more decimal places 410 00:18:23,650 --> 00:18:25,750 than anybody ever could before. 411 00:18:25,750 --> 00:18:26,530 OK? 412 00:18:26,530 --> 00:18:28,259 So that's one way to look at it. 413 00:18:28,259 --> 00:18:29,800 The other way to look at it is what's 414 00:18:29,800 --> 00:18:31,550 recommended in the Numerical Recipes book. 415 00:18:31,550 --> 00:18:34,890 And they say don't use-- 416 00:18:34,890 --> 00:18:40,320 instead of writing the equation to be chi-squared of theta 417 00:18:40,320 --> 00:18:42,967 has to be less than chi-squared max, 418 00:18:42,967 --> 00:18:45,300 they say when you're computing the confidence intervals, 419 00:18:45,300 --> 00:18:49,500 instead write chi-squared of theta 420 00:18:49,500 --> 00:19:00,256 has to be less than chi-squared best fit plus chi-squared max. 421 00:19:00,256 --> 00:19:03,280 So they say don't use a stinky little circle here, 422 00:19:03,280 --> 00:19:07,030 draw a nice big one around there. 423 00:19:10,490 --> 00:19:12,314 And what is this saying? 424 00:19:12,314 --> 00:19:14,230 I guess what this is saying is I don't believe 425 00:19:14,230 --> 00:19:18,080 that my poor agreement between the model and data 426 00:19:18,080 --> 00:19:22,730 justifies me declaring very narrow parameter values. 427 00:19:22,730 --> 00:19:23,230 So 428 00:19:23,230 --> 00:19:25,690 I just use the kinds of parameter value confidence 429 00:19:25,690 --> 00:19:29,200 intervals that it would have actually computed a priori when 430 00:19:29,200 --> 00:19:31,030 I do my experimental design. 431 00:19:31,030 --> 00:19:35,220 I think that's really how well I can determine stuff. 432 00:19:35,220 --> 00:19:37,550 And I think something is wrong that my chi-squared best 433 00:19:37,550 --> 00:19:38,630 fir is so bad. 434 00:19:38,630 --> 00:19:41,140 Maybe I didn't do my global optimization right 435 00:19:41,140 --> 00:19:43,130 and so I have the wrong parameters. 436 00:19:43,130 --> 00:19:47,360 Maybe I mis-estimated my error bars 437 00:19:47,360 --> 00:19:50,270 and so my chi-squared are all off. 438 00:19:50,270 --> 00:19:52,170 Maybe-- I don't know what else happened. 439 00:19:52,170 --> 00:19:56,210 Something happened, and I don't believe I'm that good, 440 00:19:56,210 --> 00:19:58,920 so therefore I'm going to quote much broader confidence 441 00:19:58,920 --> 00:20:02,150 intervals, which means much less precision in my determination 442 00:20:02,150 --> 00:20:03,092 of parameters. 443 00:20:03,092 --> 00:20:04,550 And that's what the Numerial Recipe 444 00:20:04,550 --> 00:20:06,780 guys recommend you should do. 445 00:20:06,780 --> 00:20:08,030 But you know, it's up to you. 446 00:20:08,030 --> 00:20:08,988 And I see it both ways. 447 00:20:08,988 --> 00:20:15,165 So the people who won the Nobel Prize for the dark energy, 448 00:20:15,165 --> 00:20:17,380 do you guys know about this? 449 00:20:17,380 --> 00:20:19,360 So if you ever saw the data there. 450 00:20:19,360 --> 00:20:20,875 So that's one where they said-- 451 00:20:20,875 --> 00:20:22,250 they were using Einstein's theory 452 00:20:22,250 --> 00:20:24,110 of relativity, General relativity. 453 00:20:24,110 --> 00:20:27,790 They say we believe Einstein, we believe that equations true. 454 00:20:27,790 --> 00:20:29,540 And we're not going to, you know, the fact 455 00:20:29,540 --> 00:20:32,210 that the data only matches when you choose 456 00:20:32,210 --> 00:20:34,910 a certain value of the cosmological constant, 457 00:20:34,910 --> 00:20:36,530 we're just going to go with it. 458 00:20:36,530 --> 00:20:38,030 And we're reporting that's the value 459 00:20:38,030 --> 00:20:39,620 of the cosmological constant. 460 00:20:39,620 --> 00:20:42,562 Now some other people might say that's kind of goofy theory. 461 00:20:42,562 --> 00:20:44,270 Einstein changed his mind six times about 462 00:20:44,270 --> 00:20:46,040 whether to put the cosmological constant 463 00:20:46,040 --> 00:20:49,100 in the equation at all, so I'm not so confident 464 00:20:49,100 --> 00:20:51,950 and I would report something different. 465 00:20:51,950 --> 00:20:53,360 But they went with their gut. 466 00:20:53,360 --> 00:20:54,680 They said we believe Einstein. 467 00:20:54,680 --> 00:20:55,665 They went with it. 468 00:20:55,665 --> 00:20:58,040 The Nobel Committee agreed, and they got the Nobel prize. 469 00:20:58,040 --> 00:21:02,060 So you know, you're taking your chances. 470 00:21:02,060 --> 00:21:04,340 The Numerical Recipe guys are more like engineers, 471 00:21:04,340 --> 00:21:06,830 they say, oh, we don't believe models that much anyway. 472 00:21:06,830 --> 00:21:11,220 Let's just use more conservative confidence intervals. 473 00:21:11,220 --> 00:21:12,270 Any questions about this? 474 00:21:15,281 --> 00:21:15,780 OK. 475 00:21:15,780 --> 00:21:16,320 All right. 476 00:21:16,320 --> 00:21:17,445 So now let's change topics. 477 00:21:23,376 --> 00:21:25,510 And we'll change to [INAUDIBLE] splitting. 478 00:21:25,510 --> 00:21:29,670 So if you remember back when we were talking about boundary 479 00:21:29,670 --> 00:21:36,439 value problems, we're often trying 480 00:21:36,439 --> 00:21:37,980 to solve problems that are like this. 481 00:22:00,658 --> 00:22:01,180 All right. 482 00:22:01,180 --> 00:22:03,420 So this is the conservation equation 483 00:22:03,420 --> 00:22:07,710 for a species in a reacting flow situation. 484 00:22:07,710 --> 00:22:11,070 So you have some convection, you have some diffusion, 485 00:22:11,070 --> 00:22:14,190 you have some reactions. 486 00:22:14,190 --> 00:22:16,290 Now in my life I have a lot of models. 487 00:22:16,290 --> 00:22:18,070 I have about 200 species in them. 488 00:22:18,070 --> 00:22:21,342 So this K is 200 different variables. 489 00:22:21,342 --> 00:22:23,550 And then you have these equations, and on top of them 490 00:22:23,550 --> 00:22:28,050 you have [INAUDIBLE] equations for that momentum conservation 491 00:22:28,050 --> 00:22:29,440 and continuity. 492 00:22:29,440 --> 00:22:33,670 And so this is actually a pretty bad set of equations. 493 00:22:33,670 --> 00:22:37,470 And typically the chemistry is stiff. 494 00:22:37,470 --> 00:22:38,970 Because I have chemistry time scales 495 00:22:38,970 --> 00:22:43,770 it'll be, like, in picoseconds, sub-nanosecond time scales. 496 00:22:43,770 --> 00:22:46,260 And I'm caring about simulating-- well in my case 497 00:22:46,260 --> 00:22:48,990 maybe I have an advanced new engine burning some fuel. 498 00:22:48,990 --> 00:22:51,690 And so we're running over a simulation of a piston cycle, 499 00:22:51,690 --> 00:22:55,724 maybe 10 milliseconds, but the chemistry is sub-nanosecond, 500 00:22:55,724 --> 00:22:56,890 so then you multiply it out. 501 00:22:56,890 --> 00:22:59,580 So that's like I'm going to 10 to minus two seconds down to 10 502 00:22:59,580 --> 00:23:00,550 to minus 10 seconds. 503 00:23:00,550 --> 00:23:03,520 That's eight orders of magnitude difference in time scales. 504 00:23:03,520 --> 00:23:07,590 So I'm in a bad way trying to solve this equation. 505 00:23:07,590 --> 00:23:11,520 And I have 200 of these guys coupled to each other. 506 00:23:11,520 --> 00:23:13,920 And then I have to think about what my mesh looks like. 507 00:23:13,920 --> 00:23:16,649 And I'm trying to simulate inside a piston, 508 00:23:16,649 --> 00:23:18,690 actually I have a moving mesh because the pistons 509 00:23:18,690 --> 00:23:21,480 compressing, so the equations are kind of complicated even 510 00:23:21,480 --> 00:23:23,227 with what the mesh is doing. 511 00:23:23,227 --> 00:23:24,810 And then I need a bunch of mesh points 512 00:23:24,810 --> 00:23:27,890 there because the-- do you guys hear the Kolmogorov scale? 513 00:23:27,890 --> 00:23:29,920 You guys learned this yet? 514 00:23:29,920 --> 00:23:32,190 It's a scale like the minimum eddy size 515 00:23:32,190 --> 00:23:33,980 if you have a turbulent flow. 516 00:23:33,980 --> 00:23:36,290 And that's pretty darn small also. 517 00:23:36,290 --> 00:23:38,660 So it might be sub-micron. 518 00:23:38,660 --> 00:23:42,500 And so I have a natural, lens scale 519 00:23:42,500 --> 00:23:44,400 that's like microns maybe. 520 00:23:44,400 --> 00:23:46,820 So 10 to minus six meters, but the piston diameter 521 00:23:46,820 --> 00:23:48,800 is couple centimeters, so again I 522 00:23:48,800 --> 00:23:52,420 have like four orders of magnitude between the mesh 523 00:23:52,420 --> 00:23:56,090 size I need to resolve the physics 524 00:23:56,090 --> 00:23:58,120 and the size of the thing. 525 00:23:58,120 --> 00:24:00,520 So I need about 10,000 points across the cylinder, 526 00:24:00,520 --> 00:24:03,020 but it's a 3-D problem, so I need 10 to the fourth, times 10 527 00:24:03,020 --> 00:24:04,645 to the fourth, times ten to the fourth, 528 00:24:04,645 --> 00:24:07,210 so I need 10 to the 12th mesh points. 529 00:24:07,210 --> 00:24:09,700 And then I have 200 variables, say variables at each mesh 530 00:24:09,700 --> 00:24:15,236 point, so I have 10 to the 14th stayed variables. 531 00:24:15,236 --> 00:24:17,110 So how much memory you have in your computer? 532 00:24:17,110 --> 00:24:17,990 How much ram? 533 00:24:17,990 --> 00:24:20,930 You have probably like 16 gigabytes. 534 00:24:20,930 --> 00:24:24,180 So you have ten to the ninth- tenth to the ninth bytes, 535 00:24:24,180 --> 00:24:25,020 is that right? 536 00:24:25,020 --> 00:24:27,741 No-- what's giga? 537 00:24:27,741 --> 00:24:28,240 Ninth? 538 00:24:28,240 --> 00:24:29,700 OK, so you have ten to the tenth. 539 00:24:29,700 --> 00:24:30,200 Sorry. 540 00:24:30,200 --> 00:24:31,940 Ten to the tenth byes, but I have 10 541 00:24:31,940 --> 00:24:35,300 to the 14th stayed variables. 542 00:24:35,300 --> 00:24:37,670 So I have a problem. 543 00:24:37,670 --> 00:24:39,670 So then you need to figure out how you solve it. 544 00:24:39,670 --> 00:24:42,400 So there's two schools of thought of how 545 00:24:42,400 --> 00:24:43,665 to solve these equations. 546 00:24:43,665 --> 00:24:45,790 The school thought which has actually been the most 547 00:24:45,790 --> 00:24:49,529 successful so far, but is not readily available, 548 00:24:49,529 --> 00:24:51,070 is to get the government to build you 549 00:24:51,070 --> 00:24:55,000 the biggest possible computer in the world with as much memory 550 00:24:55,000 --> 00:24:56,270 as possible. 551 00:24:56,270 --> 00:24:58,930 And then you can actually keep track of your state vector, 552 00:24:58,930 --> 00:25:04,240 all 10 to the 14th elements, and you 553 00:25:04,240 --> 00:25:07,460 reserve this whole computer for yourself for, your know, 554 00:25:07,460 --> 00:25:08,780 a month. 555 00:25:08,780 --> 00:25:09,810 And just run it. 556 00:25:09,810 --> 00:25:11,351 And you figure out how to parallelize 557 00:25:11,351 --> 00:25:16,220 all the calculations, and you just use an explicit solver. 558 00:25:16,220 --> 00:25:18,000 And so that way it's just straightforward. 559 00:25:18,000 --> 00:25:21,110 All you're doing is you're computing, 560 00:25:21,110 --> 00:25:27,146 you know y of T plus delta T is equal to something, 561 00:25:27,146 --> 00:25:28,520 and it's just an explicit formula 562 00:25:28,520 --> 00:25:32,120 ode45, fancy version of that. 563 00:25:32,120 --> 00:25:33,050 OK? 564 00:25:33,050 --> 00:25:36,150 And it's just algebra, and you just add them up 565 00:25:36,150 --> 00:25:38,380 and it's no problem. 566 00:25:38,380 --> 00:25:43,480 But this is a limited feasibility thing. 567 00:25:43,480 --> 00:25:45,660 This only became possible about three years ago 568 00:25:45,660 --> 00:25:48,076 when they built a computer that actually had enough memory 569 00:25:48,076 --> 00:25:49,530 to store the state factor. 570 00:25:49,530 --> 00:25:54,300 It's also-- if the equations are too stiff it still 571 00:25:54,300 --> 00:25:59,040 doesn't work, because the step size and time step you need 572 00:25:59,040 --> 00:26:00,510 to use is so small-- 573 00:26:00,510 --> 00:26:02,520 so say it's 100 picoseconds and you're 574 00:26:02,520 --> 00:26:04,800 trying to get it out to a millisecond, 575 00:26:04,800 --> 00:26:06,790 so you need 10 to the eighth something 576 00:26:06,790 --> 00:26:09,660 time steps, which means you need to double precision numbers all 577 00:26:09,660 --> 00:26:11,040 the way through. 578 00:26:11,040 --> 00:26:15,960 But if the government builds you a 64-bit 128-bit machine 579 00:26:15,960 --> 00:26:17,860 then you can do it. 580 00:26:17,860 --> 00:26:21,640 But anyway this is like a very special kind of approach. 581 00:26:21,640 --> 00:26:23,490 So now most people, including me I 582 00:26:23,490 --> 00:26:26,225 don't have access to that kind of computer power. 583 00:26:26,225 --> 00:26:27,850 Actually if you want to read about that 584 00:26:27,850 --> 00:26:32,280 there's a woman who does that, who is Jacqueline Chen, 585 00:26:32,280 --> 00:26:35,170 and you might want to look on the internet and see about her. 586 00:26:35,170 --> 00:26:37,086 And somehow she became buddies with the people 587 00:26:37,086 --> 00:26:38,650 in charge of the super computers, 588 00:26:38,650 --> 00:26:41,170 and so they're happy to give her a month of, you know, 589 00:26:41,170 --> 00:26:43,540 the best computer in the world every year, 590 00:26:43,540 --> 00:26:46,610 and she runs some really awesome calculation 591 00:26:46,610 --> 00:26:47,694 that's exactly right. 592 00:26:47,694 --> 00:26:49,610 Solves everything, has all the state variables 593 00:26:49,610 --> 00:26:52,300 ten to the 14th, and that's fantastic. 594 00:26:52,300 --> 00:26:54,580 And currently that's the benchmark 595 00:26:54,580 --> 00:26:59,500 for how people test approximate methods for how to solve this. 596 00:26:59,500 --> 00:27:02,340 Because she actually has the full solution numerically 597 00:27:02,340 --> 00:27:05,710 converged. 598 00:27:05,710 --> 00:27:08,710 But this again is a limited applicability. 599 00:27:08,710 --> 00:27:12,010 So a lot of people are working on approximation methods 600 00:27:12,010 --> 00:27:15,280 to try to solve this equation if I don't have access 601 00:27:15,280 --> 00:27:18,480 to such awesome computer resources. 602 00:27:18,480 --> 00:27:25,020 And the general idea of them-- 603 00:27:25,020 --> 00:27:26,020 well the two ideas. 604 00:27:26,020 --> 00:27:29,981 So one idea is somehow get rid of the turbulence. 605 00:27:29,981 --> 00:27:30,480 OK? 606 00:27:30,480 --> 00:27:32,960 And so what they do is they use a mesh coarser than what 607 00:27:32,960 --> 00:27:36,290 you what you would need to resolve the Kolmogorov scale, 608 00:27:36,290 --> 00:27:38,900 and then they use models for what happens 609 00:27:38,900 --> 00:27:40,910 for itty bitty tiny eddies. 610 00:27:40,910 --> 00:27:43,070 And those are called sub-grade scale models, 611 00:27:43,070 --> 00:27:45,715 and this whole approach is called large eddy simulation. 612 00:27:45,715 --> 00:27:47,840 And the idea is you keep track of the large eddies, 613 00:27:47,840 --> 00:27:51,410 the big recirculation zones that's modelled explicitly, 614 00:27:51,410 --> 00:27:53,769 but you don't keep track of all the itty bitty tiny 615 00:27:53,769 --> 00:27:55,310 swirling off little eddies that lasts 616 00:27:55,310 --> 00:27:57,518 for a little second, a little bit of a period of time 617 00:27:57,518 --> 00:27:59,934 before they dissipate because of diffusion. 618 00:27:59,934 --> 00:28:00,860 OK? 619 00:28:00,860 --> 00:28:03,590 So that's the concept. 620 00:28:03,590 --> 00:28:06,320 And there's, I don't know, half of all 621 00:28:06,320 --> 00:28:09,074 of the mechanical engineers who work in large eddy simulations 622 00:28:09,074 --> 00:28:10,740 are trying to figure out how to do this. 623 00:28:10,740 --> 00:28:12,130 Maybe Jim as well, I don't know. 624 00:28:12,130 --> 00:28:12,630 No. 625 00:28:12,630 --> 00:28:15,990 OK, but at least he knows some of them probably. 626 00:28:15,990 --> 00:28:18,442 So that's one whole branch. 627 00:28:18,442 --> 00:28:20,150 It's to try to get rid of the turbulence. 628 00:28:20,150 --> 00:28:23,120 And basically that's trying to reduce 629 00:28:23,120 --> 00:28:24,901 the density of mesh I need. 630 00:28:24,901 --> 00:28:26,150 So trying to get rid of this-- 631 00:28:26,150 --> 00:28:28,559 use the spatial mesh. 632 00:28:28,559 --> 00:28:30,100 But the other issue is the time mesh. 633 00:28:30,100 --> 00:28:31,310 And the time is even worse, right? 634 00:28:31,310 --> 00:28:33,434 Because I said I have ten to the eighth time points 635 00:28:33,434 --> 00:28:35,920 and I only need ten to the fourth spatial points 636 00:28:35,920 --> 00:28:37,050 in the dimension. 637 00:28:37,050 --> 00:28:39,380 So in some ways, if I can do something with the time, 638 00:28:39,380 --> 00:28:40,820 that really can buy me a lot. 639 00:28:40,820 --> 00:28:42,470 I have eight orders of magnitude I'd 640 00:28:42,470 --> 00:28:45,520 like to try to reduce to three orders of magnitude if I could, 641 00:28:45,520 --> 00:28:48,540 that would be really nice. 642 00:28:48,540 --> 00:28:50,600 That idea is to say, well let's separate 643 00:28:50,600 --> 00:28:53,210 the fast timescales from the slow time scales. 644 00:28:53,210 --> 00:28:54,710 And in fact all the fast time scales 645 00:28:54,710 --> 00:28:57,575 are mostly all chemistry, particularly after I changed 646 00:28:57,575 --> 00:29:00,200 the eddy simulation where I get rid of the really small spatial 647 00:29:00,200 --> 00:29:05,000 scales, then the time scale of diffusion over this big mesh 648 00:29:05,000 --> 00:29:06,500 is kind of slow. 649 00:29:06,500 --> 00:29:10,670 And so I can go to a time scale that's completely 650 00:29:10,670 --> 00:29:11,870 controlled by chemistry. 651 00:29:11,870 --> 00:29:14,328 So we already told you, how do you solve stiff differential 652 00:29:14,328 --> 00:29:15,530 equations? 653 00:29:15,530 --> 00:29:17,390 Use a stiff solver, right? 654 00:29:17,390 --> 00:29:21,097 And there are special methods, there are implicit methods, 655 00:29:21,097 --> 00:29:22,430 and there's specialized solvers. 656 00:29:22,430 --> 00:29:25,400 And what they do is they allow you to take larger time steps 657 00:29:25,400 --> 00:29:29,779 than you would be allowed to use with explicit methods 658 00:29:29,779 --> 00:29:31,820 where you would run into numerical instabilities. 659 00:29:31,820 --> 00:29:35,360 If you try to use the explicit method with a big time stamp, 660 00:29:35,360 --> 00:29:37,490 it's much bigger than the real physical time scale, 661 00:29:37,490 --> 00:29:39,140 you'll get some crazy oscillations. 662 00:29:39,140 --> 00:29:41,710 But if we use the stiff solvers you 663 00:29:41,710 --> 00:29:43,225 can get away with time steps that 664 00:29:43,225 --> 00:29:45,350 are much, much larger than the real physical scale. 665 00:29:45,350 --> 00:29:49,010 So that allows you to go from 100 picoseconds up 666 00:29:49,010 --> 00:29:51,019 to maybe 100 nanoseconds. 667 00:29:51,019 --> 00:29:53,060 Get rid of three orders of magnitude of time step 668 00:29:53,060 --> 00:29:54,950 by using a stiff solver. 669 00:29:54,950 --> 00:29:57,480 So that's pretty good, three orders magnitude reduction? 670 00:29:57,480 --> 00:29:58,190 Well go for it. 671 00:29:58,190 --> 00:30:02,716 It's not only reducing the CPU time potentially, 672 00:30:02,716 --> 00:30:04,340 but it's even more importantly reducing 673 00:30:04,340 --> 00:30:06,830 the number of points you have to save, 674 00:30:06,830 --> 00:30:09,740 because when you run this calculation what Jacqueline 675 00:30:09,740 --> 00:30:20,952 Chen ends up with is a list, y at every x, xn, yn, el, times 676 00:30:20,952 --> 00:30:25,184 e, something like this. 677 00:30:25,184 --> 00:30:28,090 Alpha, I don't know. 678 00:30:28,090 --> 00:30:29,630 That's the output of her program. 679 00:30:29,630 --> 00:30:31,330 OK, so how big is this object? 680 00:30:31,330 --> 00:30:35,855 This is 200 times 10 to the fourth, times 10 to the fourth, 681 00:30:35,855 --> 00:30:40,020 times 10 to the fourth, times 10 to the eighth. 682 00:30:40,020 --> 00:30:40,520 Right? 683 00:30:40,520 --> 00:30:42,120 That's how many numbers she's going 684 00:30:42,120 --> 00:30:43,620 to get out of her simulation. 685 00:30:43,620 --> 00:30:45,000 So, actually, her biggest problem 686 00:30:45,000 --> 00:30:47,010 is not actually solving the simulation. 687 00:30:47,010 --> 00:30:49,770 As long as she can get the one month of CPU time she's good. 688 00:30:49,770 --> 00:30:50,610 She's got to solve. 689 00:30:50,610 --> 00:30:52,990 Her Problem is how to analyze her data. 690 00:30:52,990 --> 00:30:54,815 Because now she has this much data 691 00:30:54,815 --> 00:30:56,190 that she has to figure out to say 692 00:30:56,190 --> 00:30:58,096 something rational about it. 693 00:30:58,096 --> 00:30:59,970 And so she has like a whole army of post-docs 694 00:30:59,970 --> 00:31:02,344 whose whole job is trying to invent new day interrogation 695 00:31:02,344 --> 00:31:02,970 methods. 696 00:31:02,970 --> 00:31:05,220 And also this data set is so large it cannot be stored 697 00:31:05,220 --> 00:31:06,697 on a hard disk. 698 00:31:06,697 --> 00:31:09,280 In fact it cannot be stored in the biggest farm of hard disks. 699 00:31:09,280 --> 00:31:11,790 So she has to, as the calculation is running, 700 00:31:11,790 --> 00:31:14,340 send the data to different farms of hard disks 701 00:31:14,340 --> 00:31:17,010 around the world and filling them completely up, 702 00:31:17,010 --> 00:31:18,540 and then sending them to the next one and filling them up. 703 00:31:18,540 --> 00:31:20,040 So she has fragments of this data 704 00:31:20,040 --> 00:31:22,110 set stored all over the place. 705 00:31:22,110 --> 00:31:23,670 And then if she wants to make a graph 706 00:31:23,670 --> 00:31:25,586 she has to have a special graphs software that 707 00:31:25,586 --> 00:31:27,620 looks inside this data farm, get some numbers, 708 00:31:27,620 --> 00:31:29,203 look inside this data farm and puts it 709 00:31:29,203 --> 00:31:30,366 all together in some way. 710 00:31:30,366 --> 00:31:31,740 So this is a pretty hard project. 711 00:31:31,740 --> 00:31:32,410 OK? 712 00:31:32,410 --> 00:31:34,867 Now if I can use a stiff solver I 713 00:31:34,867 --> 00:31:36,450 don't need to use so many time stamps. 714 00:31:36,450 --> 00:31:39,410 So I can get this down, say 10 to the fifth at least 715 00:31:39,410 --> 00:31:42,180 I got rid of three orders of magnitude in my data. 716 00:31:42,180 --> 00:31:42,680 Right? 717 00:31:42,680 --> 00:31:44,305 So now I just need one data farm maybe, 718 00:31:44,305 --> 00:31:46,972 instead of, you know, using all the ones at Amazon and Google 719 00:31:46,972 --> 00:31:47,680 own or something. 720 00:31:47,680 --> 00:31:48,180 OK? 721 00:31:50,690 --> 00:31:53,510 So this is a pretty big advantage to do it this way. 722 00:31:53,510 --> 00:31:55,640 Now the disadvantage, as you know, 723 00:31:55,640 --> 00:31:58,640 is solving implicit differential equations 724 00:31:58,640 --> 00:31:59,930 is way more expensive, right? 725 00:31:59,930 --> 00:32:03,320 Because you have to solve an implicit non-linear algebraic 726 00:32:03,320 --> 00:32:06,535 equation at every time step. 727 00:32:06,535 --> 00:32:08,910 And in this case it can be a pretty complicated equation. 728 00:32:08,910 --> 00:32:12,417 Because if I have this many state variables, so 10 729 00:32:12,417 --> 00:32:15,000 to the fourth, times 10 to the fourth, times 10 to the fourth, 730 00:32:15,000 --> 00:32:19,800 so 10 to 12th times 200, so 10 to the 14th at each time step, 731 00:32:19,800 --> 00:32:21,037 that's my current state. 732 00:32:21,037 --> 00:32:23,370 And I'm trying to solve a non-linear algebraic equation, 733 00:32:23,370 --> 00:32:27,300 10 to the 14th unknowns, you need a pretty good solver. 734 00:32:27,300 --> 00:32:31,200 I mean, I know backslash is good and F [? solved ?] and stuff, 735 00:32:31,200 --> 00:32:33,940 but I really don't think it's going to work. 736 00:32:33,940 --> 00:32:35,620 So that's a really serious problem. 737 00:32:35,620 --> 00:32:40,500 So the next trick is to say, OK, well this stiffness is only 738 00:32:40,500 --> 00:32:42,820 happening because of the chemistry. 739 00:32:42,820 --> 00:32:43,900 So I don't-- 740 00:32:43,900 --> 00:32:49,690 I could solve all of the finite volumes separately. 741 00:32:49,690 --> 00:32:51,730 So the idea is to split the problem up. 742 00:32:51,730 --> 00:32:55,300 I'm going to solve the chemistry using a stiff solver only 743 00:32:55,300 --> 00:32:57,550 at one point a volume at a time. 744 00:32:57,550 --> 00:33:00,410 And then I'm going to do the transport some other way. 745 00:33:00,410 --> 00:33:01,910 So what I'm going to do is I'm going 746 00:33:01,910 --> 00:33:06,330 to rewrite this equation as My unknown is dy to t, 747 00:33:06,330 --> 00:33:10,080 is equal to some kind of transport term 748 00:33:10,080 --> 00:33:12,410 and some kind of reactant term. 749 00:33:12,410 --> 00:33:17,140 And this term is local. 750 00:33:17,140 --> 00:33:20,200 All this in this equation, this is really 751 00:33:20,200 --> 00:33:30,210 same dy at position x depends on this at position x. 752 00:33:30,210 --> 00:33:31,880 Whereas if you compute these guys 753 00:33:31,880 --> 00:33:34,000 I need to do finite differences or something 754 00:33:34,000 --> 00:33:36,392 from adjoining finite volumes. 755 00:33:36,392 --> 00:33:37,122 Right? 756 00:33:37,122 --> 00:33:39,080 I have to do fluxes from the adjoining volumes. 757 00:33:39,080 --> 00:33:41,410 So these guys are sensitive not just 758 00:33:41,410 --> 00:33:43,670 to my local value of y at this finite volume, 759 00:33:43,670 --> 00:33:46,530 but also to the neighboring values of y. 760 00:33:46,530 --> 00:33:49,830 But if I don't worry about them, if I break it up 761 00:33:49,830 --> 00:33:55,110 I have this local thing and this non-local thing. 762 00:33:55,110 --> 00:33:58,470 This one is local, it's stiff. 763 00:33:58,470 --> 00:34:01,440 This one is non-local, but it's not so stiff. 764 00:34:05,900 --> 00:34:10,850 And often it's linear or nearly linear. 765 00:34:10,850 --> 00:34:13,840 Because the transport terms are mostly linear. 766 00:34:13,840 --> 00:34:15,056 And so I can use-- 767 00:34:15,056 --> 00:34:16,639 I have specialized solvers that people 768 00:34:16,639 --> 00:34:18,930 have worked on for many years for each of these things. 769 00:34:18,930 --> 00:34:21,790 So for the stiff thing I have things like ode15s, 770 00:34:21,790 --> 00:34:26,554 I have DAEPACK and [? DASPAK ?] and [? DASL ?] and sundials 771 00:34:26,554 --> 00:34:28,929 and all those programs that people have written and spent 772 00:34:28,929 --> 00:34:31,052 a lot of effort to solve this problem. 773 00:34:31,052 --> 00:34:32,510 And similarly there's a whole bunch 774 00:34:32,510 --> 00:34:33,670 of people who have spent their whole life trying 775 00:34:33,670 --> 00:34:35,980 to solve the convection diffusion problems, 776 00:34:35,980 --> 00:34:38,540 the non-reacting flow problem. 777 00:34:38,540 --> 00:34:40,795 And so I can use specialized solvers that they have. 778 00:34:40,795 --> 00:34:41,920 So I want to split this up. 779 00:34:41,920 --> 00:34:43,461 And what I really would like to solve 780 00:34:43,461 --> 00:34:48,670 is I'd prefer to solve dy/dt equals r and dy/dt equals t 781 00:34:48,670 --> 00:34:50,920 completely separately if I could. 782 00:34:50,920 --> 00:34:52,510 Because I can use specialized methods 783 00:34:52,510 --> 00:34:54,175 that are the best for each. 784 00:34:54,175 --> 00:34:56,050 So not to figure out, well this is not really 785 00:34:56,050 --> 00:34:58,880 the real equation because it has both of them added together. 786 00:34:58,880 --> 00:35:00,590 How can I handle that? 787 00:35:00,590 --> 00:35:06,670 So this operator splitting recipe 788 00:35:06,670 --> 00:35:15,910 is solve dy/dt equals to t for a while, 789 00:35:15,910 --> 00:35:21,250 and at the end of this process I'll have a [? yat ?] final. 790 00:35:21,250 --> 00:35:27,860 And then I'll solve this equation with y equals y0. 791 00:35:27,860 --> 00:35:30,740 That's the equation system of solving, of transport. 792 00:35:30,740 --> 00:35:33,670 I do some transport for a while, I get some values for the y, 793 00:35:33,670 --> 00:35:39,220 then I take these guys and I solve 794 00:35:39,220 --> 00:35:44,810 the reaction equation starting from y equals y final. 795 00:35:44,810 --> 00:35:46,630 How to write this? 796 00:35:46,630 --> 00:35:47,830 One or two 0. 797 00:35:51,190 --> 00:35:52,150 OK? 798 00:35:52,150 --> 00:35:56,050 And I solve this one for a while and I'll get a different, 799 00:35:56,050 --> 00:35:58,621 you know, a somewhat different [? yrt ?] final. 800 00:36:01,160 --> 00:36:01,660 All right? 801 00:36:01,660 --> 00:36:02,770 And the simplest way of putting things 802 00:36:02,770 --> 00:36:04,186 is just keep going back and forth. 803 00:36:04,186 --> 00:36:07,230 Boom, boom, boom, boom, boom, boom, boom, boom. 804 00:36:07,230 --> 00:36:11,600 And in the limit where I make this t final minus 805 00:36:11,600 --> 00:36:17,110 t0 to be very small, then maybe I 806 00:36:17,110 --> 00:36:19,950 might expect that they might converge to the same solution 807 00:36:19,950 --> 00:36:23,680 as the couple of questions. 808 00:36:23,680 --> 00:36:24,190 All right? 809 00:36:24,190 --> 00:36:26,570 So that's the concept. 810 00:36:26,570 --> 00:36:29,550 Now, if you do it this way it doesn't work very well. 811 00:36:29,550 --> 00:36:32,950 You need to make this delta t incredibly small 812 00:36:32,950 --> 00:36:34,920 to get any accuracy. 813 00:36:34,920 --> 00:36:38,400 However you can do a slightly more clever version, which 814 00:36:38,400 --> 00:36:46,775 is I saw dy/dt equal t starting from yt0, so y0, 815 00:36:46,775 --> 00:36:56,100 and I solve this to me y at t plus delta-t over 2. 816 00:36:56,100 --> 00:37:04,330 Then I solve dy/dt is equal to r, starting from y of t0 817 00:37:04,330 --> 00:37:09,820 equal to, call this y star. 818 00:37:09,820 --> 00:37:15,660 y star and this I integrate out to y 819 00:37:15,660 --> 00:37:17,925 of t plus the full delta t. 820 00:37:17,925 --> 00:37:21,690 I'll call that y double star. 821 00:37:21,690 --> 00:37:28,800 And then I solve this one again, and I 822 00:37:28,800 --> 00:37:33,960 solve starting from y of t0 equal to y double star, 823 00:37:33,960 --> 00:37:36,380 and I integrate that out for another delta t over 2. 824 00:37:40,680 --> 00:37:44,680 And this one I called my y final. 825 00:37:44,680 --> 00:37:47,350 And that's the one I use. 826 00:37:47,350 --> 00:37:49,654 So I take a half of a transport step, 827 00:37:49,654 --> 00:37:51,570 sort of like I'm starting to do the transport, 828 00:37:51,570 --> 00:37:52,910 then I suddenly-- oh, hold on. 829 00:37:52,910 --> 00:37:54,020 I got you my reaction. 830 00:37:54,020 --> 00:37:56,950 I do my reaction, then I do the other half of my transport 831 00:37:56,950 --> 00:37:58,570 after that. 832 00:37:58,570 --> 00:38:00,430 And if you can do the analysis of this, 833 00:38:00,430 --> 00:38:02,070 it turns out this thing converges 834 00:38:02,070 --> 00:38:04,460 to the true solution of the couple equations 835 00:38:04,460 --> 00:38:06,030 with second order accuracy. 836 00:38:06,030 --> 00:38:09,390 So you might make this delta t, if I cut it 837 00:38:09,390 --> 00:38:12,477 in half the precision increased by factor of four. 838 00:38:12,477 --> 00:38:15,060 If I cut if by a factor of 10 it increased by a factor of 100. 839 00:38:15,060 --> 00:38:17,950 So this is looking more promising as a way to do it. 840 00:38:17,950 --> 00:38:19,930 This is the most popular way to do it. 841 00:38:19,930 --> 00:38:24,750 So this recipe here for operator splitting 842 00:38:24,750 --> 00:38:27,230 is called Strang splitting. 843 00:38:27,230 --> 00:38:30,780 Named after Gill Strang, he was a professor 844 00:38:30,780 --> 00:38:33,300 in applied mathematics department here. 845 00:38:33,300 --> 00:38:35,970 Former head of CIAM, the international body 846 00:38:35,970 --> 00:38:37,710 for applied mathematics. 847 00:38:37,710 --> 00:38:39,600 And also author of a fantastically great 848 00:38:39,600 --> 00:38:41,766 linear algebra textbook which maybe some of you guys 849 00:38:41,766 --> 00:38:44,320 might have seen some videos of his. 850 00:38:44,320 --> 00:38:48,650 So he invented this 1968, I was five years old. 851 00:38:48,650 --> 00:38:51,120 He didn't know that his invention 852 00:38:51,120 --> 00:38:54,340 was of any practical use, being a mathematician. 853 00:38:54,340 --> 00:38:56,102 So he just published that paper. 854 00:38:56,102 --> 00:38:58,560 Now in the meantime, this paper is like a citation classic. 855 00:38:58,560 --> 00:39:01,150 It's cited like, I don't know, 20,000 times. 856 00:39:01,150 --> 00:39:02,350 He had no idea. 857 00:39:02,350 --> 00:39:06,050 He never looks at citation list, he's a mathematician. 858 00:39:06,050 --> 00:39:09,880 So I was working on this problem and-- 859 00:39:12,680 --> 00:39:14,680 let's see-- let's back up. 860 00:39:14,680 --> 00:39:16,191 So this is how he usually solves it. 861 00:39:16,191 --> 00:39:16,690 Right? 862 00:39:16,690 --> 00:39:18,040 We usually solve the problems, the-- 863 00:39:18,040 --> 00:39:18,748 splitting, sorry. 864 00:39:18,748 --> 00:39:20,129 Just be able to see it. 865 00:39:22,720 --> 00:39:25,154 And we like to solve it this way because I 866 00:39:25,154 --> 00:39:26,320 can use specialized solvers. 867 00:39:26,320 --> 00:39:29,290 I use a transport solver, I use a stiff solver here. 868 00:39:29,290 --> 00:39:31,000 This one I can parallelize because I 869 00:39:31,000 --> 00:39:34,902 can do a different ode solution at each finite volume. 870 00:39:34,902 --> 00:39:37,110 So if I have a big computer, a big parallel computer, 871 00:39:37,110 --> 00:39:39,970 say I have 10,000 nodes, I can solve 10,000 finite volumes 872 00:39:39,970 --> 00:39:43,260 simultaneously in this step, so the CPU time 873 00:39:43,260 --> 00:39:46,060 is reduced a lot because it parallelizes easily. 874 00:39:46,060 --> 00:39:47,230 So this is really popular. 875 00:39:47,230 --> 00:39:50,530 This is the main way that people do things. 876 00:39:50,530 --> 00:39:52,600 And then in the ideal case you do it 877 00:39:52,600 --> 00:39:54,610 with a certain guess of delta t, then 878 00:39:54,610 --> 00:39:56,680 you try delta t divided by five. 879 00:39:56,680 --> 00:39:59,686 And if you get the same result then you publish it. 880 00:39:59,686 --> 00:40:00,342 OK? 881 00:40:00,342 --> 00:40:01,800 And you've converged your solution. 882 00:40:06,420 --> 00:40:11,520 However, what I found doing this is I was trying to find flames. 883 00:40:11,520 --> 00:40:14,500 I was trying to calculate the positions of flames. 884 00:40:14,500 --> 00:40:18,610 And what I found was that with the delta t's that I would try, 885 00:40:18,610 --> 00:40:21,510 sometimes I would get a qualitatively different 886 00:40:21,510 --> 00:40:25,750 solution depending on the delta t I chose. 887 00:40:25,750 --> 00:40:27,190 And so in particular what I found 888 00:40:27,190 --> 00:40:30,710 was that I was solving burner flames. 889 00:40:30,710 --> 00:40:34,270 So I have, like, a burner here, and it had a flame on top of it 890 00:40:34,270 --> 00:40:36,110 like this. 891 00:40:36,110 --> 00:40:38,200 And there's two importantly different cases. 892 00:40:38,200 --> 00:40:41,380 One is where the flame is floating above the burner. 893 00:40:41,380 --> 00:40:43,520 And one is where the flame kind of goes all the way 894 00:40:43,520 --> 00:40:45,525 and touches the metal. 895 00:40:45,525 --> 00:40:47,275 And you've probably seen both those cases, 896 00:40:47,275 --> 00:40:49,460 like a bunsen burner or something like that. 897 00:40:49,460 --> 00:40:51,900 Sometimes if you get the flow rates high enough, 898 00:40:51,900 --> 00:40:54,680 the whole thing lift off and be floating there. 899 00:40:54,680 --> 00:40:57,640 And sometimes if you slow the flow rates down 900 00:40:57,640 --> 00:41:00,920 they say that they will anchor to the metal. 901 00:41:00,920 --> 00:41:03,620 Some piece the flame will actually touch the metal. 902 00:41:03,620 --> 00:41:06,412 Usually it's like if you have a-- 903 00:41:06,412 --> 00:41:08,120 a lot of times you have a case like this. 904 00:41:08,120 --> 00:41:10,940 You have like a metal piece, metal piece, 905 00:41:10,940 --> 00:41:14,990 here's your mixture coming in, out here is the air, 906 00:41:14,990 --> 00:41:17,390 and a lot of times the flame will actually 907 00:41:17,390 --> 00:41:18,797 touch this metal piece here. 908 00:41:18,797 --> 00:41:20,630 You know, the flame sort of looks like that. 909 00:41:20,630 --> 00:41:22,171 Have you ever seen a flame like this? 910 00:41:22,171 --> 00:41:24,230 Anyways, this is a common thing, anchored-flame. 911 00:41:24,230 --> 00:41:26,684 If you just have a slow speed burner that's what you get. 912 00:41:26,684 --> 00:41:28,100 If you go high speed on the burner 913 00:41:28,100 --> 00:41:29,090 you can make the whole thing lift up and sort of 914 00:41:29,090 --> 00:41:31,640 dance around if you look at the flame. 915 00:41:31,640 --> 00:41:32,140 All right? 916 00:41:32,140 --> 00:41:33,931 You could probably see the same if you look 917 00:41:33,931 --> 00:41:35,510 at a piece of wood burning. 918 00:41:35,510 --> 00:41:38,030 If it's burning well the flame will be lifted off the wood, 919 00:41:38,030 --> 00:41:40,312 you'll see a gap between the flame and the wood. 920 00:41:40,312 --> 00:41:42,770 If it's not burning that well you can see bits of the flame 921 00:41:42,770 --> 00:41:44,870 actually touch the pieces of wood. 922 00:41:44,870 --> 00:41:45,681 OK? 923 00:41:45,681 --> 00:41:47,930 So that's a very important thing for a commission guy. 924 00:41:47,930 --> 00:41:50,430 So was really-- to me this is a big deal. 925 00:41:50,430 --> 00:41:52,700 So when I solved it with one choice of delta t 926 00:41:52,700 --> 00:41:53,990 I get the solution. 927 00:41:53,990 --> 00:41:56,000 When I solved with another choice of delta t 928 00:41:56,000 --> 00:42:00,390 I get the solution, with a nice big gap. 929 00:42:00,390 --> 00:42:02,195 Now they both can't be right. 930 00:42:02,195 --> 00:42:03,570 And so I was a little bit alarmed 931 00:42:03,570 --> 00:42:06,990 because it's like, it's really different, right? 932 00:42:06,990 --> 00:42:09,360 So I thought, you know, I didn't really 933 00:42:09,360 --> 00:42:11,937 care that much about the dynamics of my problem. 934 00:42:11,937 --> 00:42:14,520 I just wanted to make sure that the steady state flame I ended 935 00:42:14,520 --> 00:42:16,660 up with was the right one. 936 00:42:16,660 --> 00:42:18,660 So if you're interested in steady state problems 937 00:42:18,660 --> 00:42:19,500 this is a really important thing. 938 00:42:19,500 --> 00:42:21,874 You want to make sure that your solution really converges 939 00:42:21,874 --> 00:42:23,010 to the true steady state. 940 00:42:23,010 --> 00:42:25,634 So what I demonstrated primarily is the Strang splitting method 941 00:42:25,634 --> 00:42:28,410 is not reliable to converge steady state. 942 00:42:28,410 --> 00:42:30,529 Now I believe I made the delta t tiny enough, 943 00:42:30,529 --> 00:42:32,820 probably I would converge the [INAUDIBLE] steady state. 944 00:42:32,820 --> 00:42:34,450 But I don't know how tiny I have to go, 945 00:42:34,450 --> 00:42:36,450 and I can't really judge my convergence criteria 946 00:42:36,450 --> 00:42:38,780 very well, because at some point of delta t 947 00:42:38,780 --> 00:42:41,474 is it suddenly jumps from one solution to another solution 948 00:42:41,474 --> 00:42:43,390 that's qualitatively and completely different. 949 00:42:43,390 --> 00:42:46,277 So the convergence is not smooth and continuous, 950 00:42:46,277 --> 00:42:48,360 it's like it's converging, converging, converging, 951 00:42:48,360 --> 00:42:50,594 jumps to some other solution. 952 00:42:50,594 --> 00:42:52,260 And I don't know where the jump happens. 953 00:42:52,260 --> 00:42:53,460 And I don't know if I keep going further, 954 00:42:53,460 --> 00:42:55,120 will it jump to something else? 955 00:42:55,120 --> 00:42:57,280 So this makes me very loathe to publish the paper. 956 00:42:57,280 --> 00:42:58,260 Right? 957 00:42:58,260 --> 00:42:59,412 So I don't know. 958 00:42:59,412 --> 00:43:00,620 So I was agitated about that. 959 00:43:00,620 --> 00:43:03,547 So I start working on trying to fix this splitting. 960 00:43:03,547 --> 00:43:05,380 The problem was that Strang is really smart. 961 00:43:05,380 --> 00:43:07,005 So his splitting method is really good. 962 00:43:07,005 --> 00:43:09,880 It's stable, it's pretty accurate, 963 00:43:09,880 --> 00:43:11,430 it's really hard to beat it actually. 964 00:43:11,430 --> 00:43:13,638 And so I published the papers and they were that good 965 00:43:13,638 --> 00:43:16,200 because they really didn't really improve it that much. 966 00:43:16,200 --> 00:43:17,783 But then I had a really smart post-doc 967 00:43:17,783 --> 00:43:21,270 named Ray Speth, now works in the AeroAstro department. 968 00:43:21,270 --> 00:43:24,930 And he came up with a different way to look at this. 969 00:43:24,930 --> 00:43:29,450 And he said, you know, if I had a problem that's like this, 970 00:43:29,450 --> 00:43:33,230 I can always add some constant to this 971 00:43:33,230 --> 00:43:34,940 and subtract the same constant from this, 972 00:43:34,940 --> 00:43:38,290 and I'm really getting the same problem. 973 00:43:38,290 --> 00:43:40,050 And so now the thing is, can I cleverly 974 00:43:40,050 --> 00:43:44,670 choose this constant so that, for example, my steady state 975 00:43:44,670 --> 00:43:47,000 solution is going to be good? 976 00:43:47,000 --> 00:43:49,350 So how to figure out what constants to choose? 977 00:43:49,350 --> 00:43:52,710 And so the methods where you add these constants 978 00:43:52,710 --> 00:43:54,520 are called balancing methods. 979 00:43:54,520 --> 00:43:57,700 So it's trying to balance the operator splitting scheme. 980 00:43:57,700 --> 00:44:00,970 And what Ray's method-- 981 00:44:05,850 --> 00:44:07,440 pretty simple balancing method. 982 00:44:07,440 --> 00:44:09,291 So this is called simple balancing. 983 00:44:19,490 --> 00:44:22,120 Let's see if I have a slide that shows it. 984 00:44:22,120 --> 00:44:22,620 There we go. 985 00:44:22,620 --> 00:44:23,400 Here we go. 986 00:44:23,400 --> 00:44:24,450 Simple balancing. 987 00:44:24,450 --> 00:44:30,300 You choose a c to be the average of your right hand 988 00:44:30,300 --> 00:44:34,482 side at the previous time step. 989 00:44:34,482 --> 00:44:35,940 OK, now idea here is that as you're 990 00:44:35,940 --> 00:44:39,000 getting close to a steady state, the previous time step 991 00:44:39,000 --> 00:44:41,375 is not going be that much different than the current time 992 00:44:41,375 --> 00:44:41,920 step. 993 00:44:41,920 --> 00:44:43,378 And so this average number is going 994 00:44:43,378 --> 00:44:46,230 to stay more or less study as you as you get close. 995 00:44:46,230 --> 00:44:55,186 And so what you're really doing is adding the, sort of, 996 00:44:55,186 --> 00:44:57,060 the average of these two guys to each of them 997 00:44:57,060 --> 00:44:59,330 from the previous step, and that has 998 00:44:59,330 --> 00:45:01,440 the consequence of really changing the solution. 999 00:45:01,440 --> 00:45:04,270 So let me just show you. 1000 00:45:04,270 --> 00:45:06,900 So the left hand side is what Strang splitting does. 1001 00:45:14,050 --> 00:45:17,220 The blue chips are the transport, 1002 00:45:17,220 --> 00:45:20,400 and the green steps are the reaction. 1003 00:45:20,400 --> 00:45:23,100 And what's happening is I did a half a transport step, 1004 00:45:23,100 --> 00:45:24,090 then I do a reaction. 1005 00:45:24,090 --> 00:45:29,867 The reaction sort of jumps me off the trajectory, 1006 00:45:29,867 --> 00:45:31,450 then brings me back to the trajectory, 1007 00:45:31,450 --> 00:45:35,830 then overshoots because I do a whole delta t step of reaction. 1008 00:45:35,830 --> 00:45:39,010 And then I have to do a half step of transport back. 1009 00:45:39,010 --> 00:45:41,230 And so I get from the black-- one black dot 1010 00:45:41,230 --> 00:45:43,650 is the first point, and the other black dot 1011 00:45:43,650 --> 00:45:47,170 is what I call y final. 1012 00:45:47,170 --> 00:45:48,537 And that's showing what happens. 1013 00:45:48,537 --> 00:45:50,620 And it just repeats like that over and over again. 1014 00:45:50,620 --> 00:45:51,460 And it makes sense. 1015 00:45:51,460 --> 00:45:53,750 If you're close to a steady state, 1016 00:45:53,750 --> 00:45:58,600 then if I'd just take one of these terms, 1017 00:45:58,600 --> 00:46:00,450 without the balancing theta transport, 1018 00:46:00,450 --> 00:46:02,200 it's going to head away from steady state, 1019 00:46:02,200 --> 00:46:04,750 because at steady state the transport and the reaction 1020 00:46:04,750 --> 00:46:06,730 are balancing against each other. 1021 00:46:06,730 --> 00:46:08,590 Typically you have something transporting in 1022 00:46:08,590 --> 00:46:10,512 and is being consumed by the reaction. 1023 00:46:10,512 --> 00:46:12,220 Or something is being created by reaction 1024 00:46:12,220 --> 00:46:13,472 and is being transported away. 1025 00:46:13,472 --> 00:46:14,680 And they have opposite signs. 1026 00:46:14,680 --> 00:46:16,660 If I only put one in that has a certain sign, 1027 00:46:16,660 --> 00:46:18,760 then I'm going to move away from the steady state. 1028 00:46:18,760 --> 00:46:20,830 And then I need the other one to bring it back. 1029 00:46:20,830 --> 00:46:22,810 And the reason why Strang splitting has trouble 1030 00:46:22,810 --> 00:46:25,450 is it's doing these huge excursions away 1031 00:46:25,450 --> 00:46:29,050 from the steady state at each time step, 1032 00:46:29,050 --> 00:46:31,670 with the hope that it ends up back at the steady state. 1033 00:46:31,670 --> 00:46:32,170 Right? 1034 00:46:32,170 --> 00:46:33,879 That's the way it's supposed to work out. 1035 00:46:33,879 --> 00:46:35,336 And to second order it does, that's 1036 00:46:35,336 --> 00:46:37,370 because it's pretty clever about how it works. 1037 00:46:37,370 --> 00:46:39,780 But it does it's crazy excursions. 1038 00:46:39,780 --> 00:46:41,410 Now if you have a non-linear model, 1039 00:46:41,410 --> 00:46:44,290 the further you do excursions from the real trajectory 1040 00:46:44,290 --> 00:46:45,940 the bigger the area you get. 1041 00:46:45,940 --> 00:46:51,730 So if you do that balance method, 1042 00:46:51,730 --> 00:46:56,537 you can see what happens is that when you get to steady state 1043 00:46:56,537 --> 00:46:58,870 all the steps are just moving you along the steady state 1044 00:46:58,870 --> 00:46:59,462 trajectory. 1045 00:46:59,462 --> 00:47:01,420 They're basically maintaining the steady state, 1046 00:47:01,420 --> 00:47:06,100 and that's because I added-- 1047 00:47:06,100 --> 00:47:14,510 it's like transport plus half of transport 1048 00:47:14,510 --> 00:47:17,710 plus r or something like that. 1049 00:47:17,710 --> 00:47:21,050 And what is it? 1050 00:47:21,050 --> 00:47:25,200 There's probably a minus sign here. 1051 00:47:25,200 --> 00:47:29,420 And so I'm getting rid of half of my transport 1052 00:47:29,420 --> 00:47:32,730 and I'm adding in half of this. 1053 00:47:32,730 --> 00:47:40,550 So I changed-- before I had dy/dt is equal to t, 1054 00:47:40,550 --> 00:47:44,460 but now I have dy/dt is equal to 1/2 of t. 1055 00:47:52,850 --> 00:47:54,460 I feel I need a sign to stick here. 1056 00:47:58,920 --> 00:48:06,020 So half of t minus f of r, something like that. 1057 00:48:06,020 --> 00:48:12,930 And the-- this can't be right. 1058 00:48:12,930 --> 00:48:13,910 I must have a sign. 1059 00:48:13,910 --> 00:48:16,450 Maybe this is not the average it's a half. 1060 00:48:20,680 --> 00:48:21,830 It should be minus. 1061 00:48:21,830 --> 00:48:25,029 It's like the-- these things should be added, 1062 00:48:25,029 --> 00:48:26,570 and these guys are you cancelled out, 1063 00:48:26,570 --> 00:48:28,445 so basically [INAUDIBLE] is going to be zero, 1064 00:48:28,445 --> 00:48:30,290 which it should be as a steady state. 1065 00:48:30,290 --> 00:48:33,760 So I'm sorry, that previous slide had the sign wrong. 1066 00:48:33,760 --> 00:48:36,330 This is what happens as the solution. 1067 00:48:36,330 --> 00:48:39,740 If you do Strang splitting on the left 1068 00:48:39,740 --> 00:48:47,927 you get more accuracy as you cut delta t, but it's only-- 1069 00:48:47,927 --> 00:48:49,630 it's not getting any better. 1070 00:48:49,630 --> 00:48:50,850 It doesn't matter if it's near steady state 1071 00:48:50,850 --> 00:48:53,183 or away from steady state, it treats everything equally, 1072 00:48:53,183 --> 00:48:54,320 the Strang splitting. 1073 00:48:54,320 --> 00:48:57,442 In the balance splittings, as you get to steady state, 1074 00:48:57,442 --> 00:48:59,400 all of a sudden the errors drive exponentially. 1075 00:48:59,400 --> 00:49:02,180 You see this litte logarithm of the error? 1076 00:49:02,180 --> 00:49:03,690 Logarithmic plot. 1077 00:49:03,690 --> 00:49:06,510 So you get tremendously great convergence 1078 00:49:06,510 --> 00:49:09,240 to the real steady state once you get anywhere close 1079 00:49:09,240 --> 00:49:11,620 to the steady state. 1080 00:49:11,620 --> 00:49:17,980 Now the cost of this is that because the formula, which 1081 00:49:17,980 --> 00:49:22,880 has a sign mistake I think, it depends on yn, 1082 00:49:22,880 --> 00:49:26,370 it makes the whole thing more numerically unstable. 1083 00:49:26,370 --> 00:49:27,870 Because it's like an explicit thing. 1084 00:49:27,870 --> 00:49:29,703 We're using what happened at a previous time 1085 00:49:29,703 --> 00:49:32,280 to correct our equation for future times. 1086 00:49:32,280 --> 00:49:35,100 And so you are more likely to get numerical instabilities. 1087 00:49:35,100 --> 00:49:36,870 So then Ray invented an implicit version 1088 00:49:36,870 --> 00:49:38,536 of that, and also a higher order method, 1089 00:49:38,536 --> 00:49:40,280 and that's called re-balance splitting. 1090 00:49:40,280 --> 00:49:41,440 The equation is way too complicated for me 1091 00:49:41,440 --> 00:49:43,939 to write down, since I can't even write down the minus signs 1092 00:49:43,939 --> 00:49:45,380 correctly in this one. 1093 00:49:45,380 --> 00:49:47,810 But you're happy to read the paper, it's in CIAM. 1094 00:49:47,810 --> 00:49:49,340 And he has the equations for that. 1095 00:49:49,340 --> 00:49:51,990 And when you put the second order one in, that's implicit, 1096 00:49:51,990 --> 00:49:55,140 then the thing is just as stable as Strang splitting, 1097 00:49:55,140 --> 00:49:57,360 but it also has this great property 1098 00:49:57,360 --> 00:49:59,370 that as you get close to steady state 1099 00:49:59,370 --> 00:50:02,160 it exponentially converges to the solution, which 1100 00:50:02,160 --> 00:50:05,302 is sort of like how Newton's method is really so good. 1101 00:50:05,302 --> 00:50:06,260 It's sort of like that. 1102 00:50:06,260 --> 00:50:09,120 You get a lot more decimal places really fast 1103 00:50:09,120 --> 00:50:11,740 once you get close to steady state. 1104 00:50:11,740 --> 00:50:13,050 All right, that's enough. 1105 00:50:13,050 --> 00:50:14,730 Enjoy Thanksgiving. 1106 00:50:14,730 --> 00:50:17,050 I posted a bunch of readings and things 1107 00:50:17,050 --> 00:50:20,675 about metropolis Monte Carlo and stochastic stuff 1108 00:50:20,675 --> 00:50:23,050 that we're talking about next week, and also a little bit 1109 00:50:23,050 --> 00:50:25,214 more about models versus data. 1110 00:50:25,214 --> 00:50:26,880 So do the reading when you get a chance. 1111 00:50:26,880 --> 00:50:29,800 Say hi to your families, say hi to your friends. 1112 00:50:29,800 --> 00:50:31,258 Have a nice time.