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,320 at ocw.mit.edu. 8 00:00:23,289 --> 00:00:24,830 WILLIAM GREEN JR: All right, so we're 9 00:00:24,830 --> 00:00:33,260 going to start a new topic today about ordinary differential 10 00:00:33,260 --> 00:00:35,555 equations initial value problems. 11 00:00:35,555 --> 00:00:39,020 But before I start to talk about that, 12 00:00:39,020 --> 00:00:41,840 I want to remind you next week's schedule is weird. 13 00:00:41,840 --> 00:00:45,830 So there's no classes on Monday, but instead we'll 14 00:00:45,830 --> 00:00:48,860 have Monday class on Tuesday. 15 00:00:48,860 --> 00:00:51,950 So I'll see you next on Tuesday morning. 16 00:00:51,950 --> 00:00:55,980 And then we're not going to have a class on Wednesday, 17 00:00:55,980 --> 00:00:59,490 but we have the quiz on Wednesday night. 18 00:00:59,490 --> 00:01:02,250 And then we'll get back to normal on Friday. 19 00:01:02,250 --> 00:01:04,709 And I would expect that possibly the TAs might 20 00:01:04,709 --> 00:01:08,770 be interested in giving a help session or review on Wednesday 21 00:01:08,770 --> 00:01:10,137 the class period. 22 00:01:10,137 --> 00:01:11,470 AUDIENCE: We'll be here at 3:00. 23 00:01:11,470 --> 00:01:11,910 WILLIAM GREEN JR: They'll be here. 24 00:01:11,910 --> 00:01:13,390 So if you want to come at class period, 25 00:01:13,390 --> 00:01:14,764 and just ask questions and stuff. 26 00:01:14,764 --> 00:01:18,140 AUDIENCE: And I'll have your PSAT graded back, if you like. 27 00:01:18,140 --> 00:01:19,750 WILLIAM GREEN JR: All right, got that? 28 00:01:19,750 --> 00:01:22,400 PSAT be ready as well. 29 00:01:22,400 --> 00:01:26,850 OK, so I'll talk about today ODE-IVPs. 30 00:01:26,850 --> 00:01:31,080 AUDIENCE: Yeah, we'll still have office hours Monday [INAUDIBLE] 31 00:01:31,080 --> 00:01:33,170 WILLIAM GREEN JR: Can you guys hear this? 32 00:01:33,170 --> 00:01:35,210 Yes, office hours like normal Monday. 33 00:01:38,047 --> 00:01:40,130 And for those of you who want to get your homework 34 00:01:40,130 --> 00:01:41,921 started early, we might even post homework. 35 00:01:54,290 --> 00:02:01,632 And I'm going to focus exclusively on ODE-IVPs 36 00:02:01,632 --> 00:02:02,590 you can write this way. 37 00:02:11,600 --> 00:02:16,807 All right, and I would just comment, 38 00:02:16,807 --> 00:02:17,890 what do these things mean? 39 00:02:17,890 --> 00:02:20,060 ODE means ordinary differential equation. 40 00:02:20,060 --> 00:02:22,856 It means that the only differentials in the problem 41 00:02:22,856 --> 00:02:24,230 are with respect to one variable, 42 00:02:24,230 --> 00:02:27,200 so I'm going to call that t. 43 00:02:27,200 --> 00:02:29,070 It can be z or x or whatever you want, 44 00:02:29,070 --> 00:02:31,190 but I'm going to call it t here. 45 00:02:31,190 --> 00:02:34,220 And so this is not partial differential equations 46 00:02:34,220 --> 00:02:38,030 where you have differentials respect to many variables. 47 00:02:38,030 --> 00:02:44,460 Initial value problems mean that all the initial conditions 48 00:02:44,460 --> 00:02:46,780 are given in at a single t point where 49 00:02:46,780 --> 00:02:50,402 t is that the variable that appears in the differential. 50 00:02:53,100 --> 00:02:55,100 All right so this is the form we're going to do. 51 00:02:55,100 --> 00:02:59,880 This is first order a way of writing the ODEs. 52 00:02:59,880 --> 00:03:03,965 If you recall from homework zero problem number one, 53 00:03:03,965 --> 00:03:06,590 we showed you there how you can write higher order differential 54 00:03:06,590 --> 00:03:07,250 equations-- 55 00:03:10,580 --> 00:03:12,650 convert them into the form of first order. 56 00:03:16,100 --> 00:03:21,134 And I guess that's worthwhile to point down. 57 00:03:21,134 --> 00:03:23,050 I would also comment that I going to just talk 58 00:03:23,050 --> 00:03:28,170 about f of y, however the ODE solvers in Matlab 59 00:03:28,170 --> 00:03:34,600 expect f of the t comma y as their from. 60 00:03:34,600 --> 00:03:42,820 All right, so let's just explain all the details about this. 61 00:03:42,820 --> 00:03:47,900 So suppose I had a differential equation like this. 62 00:03:58,887 --> 00:04:02,080 All right, so this is like a differential equation you'd 63 00:04:02,080 --> 00:04:06,260 have for a body move moving-- 64 00:04:06,260 --> 00:04:09,520 a particle moving with some friction force on it 65 00:04:09,520 --> 00:04:11,760 and some time varying force. 66 00:04:11,760 --> 00:04:12,378 Yeah? 67 00:04:12,378 --> 00:04:14,170 AUDIENCE: Are you going to be posting your notes online? 68 00:04:14,170 --> 00:04:15,310 WILLIAM GREEN JR: No, I don't have any notes. 69 00:04:15,310 --> 00:04:16,268 This is it. 70 00:04:20,579 --> 00:04:26,180 All right, so this is a differential equation. 71 00:04:26,180 --> 00:04:29,420 What you can do is to convert it into the form I have above, 72 00:04:29,420 --> 00:04:35,100 is to write a new variable, say v, that's defined to be dx dt. 73 00:04:35,100 --> 00:04:38,250 That's the velocity. 74 00:04:38,250 --> 00:04:41,344 And has a time dependence, but I told you 75 00:04:41,344 --> 00:04:43,010 I don't like time dependence, so I'm not 76 00:04:43,010 --> 00:04:44,640 going to write it that way. 77 00:04:44,640 --> 00:04:47,340 So we can also define another variable. 78 00:04:47,340 --> 00:04:57,950 So we can define our new y vector to be x, v, and t. 79 00:05:01,110 --> 00:05:02,850 So this is y-- the first component of y, 80 00:05:02,850 --> 00:05:04,266 this is the second component of y, 81 00:05:04,266 --> 00:05:06,210 this is the third component of y. 82 00:05:06,210 --> 00:05:10,110 And then I can write down what f of y is. 83 00:05:22,000 --> 00:05:27,190 So the equation for dx dt is y2. 84 00:05:27,190 --> 00:05:32,260 It's the second component of y, dx dt is equal to v. 85 00:05:32,260 --> 00:05:38,750 The equation for the dv dt is given by this equation. 86 00:05:38,750 --> 00:05:53,442 So that's f of y3 over m minus gamma y2. 87 00:05:53,442 --> 00:05:57,600 And the equation dt dt is just one. 88 00:06:00,390 --> 00:06:03,360 So this is my f of y. 89 00:06:03,360 --> 00:06:05,780 So this is f1, f2, f3. 90 00:06:05,780 --> 00:06:06,980 Is that all right? 91 00:06:10,114 --> 00:06:11,655 So this is just how you can convert-- 92 00:06:11,655 --> 00:06:13,488 if somebody writes you an equation this way, 93 00:06:13,488 --> 00:06:16,520 and you want to get it into the standard form that I'm going 94 00:06:16,520 --> 00:06:19,700 to show you here-- 95 00:06:19,700 --> 00:06:21,980 dy dt equals f of y-- 96 00:06:21,980 --> 00:06:22,980 this is how you convert. 97 00:06:25,710 --> 00:06:29,747 And a skill we haven't really emphasized so far 98 00:06:29,747 --> 00:06:31,830 in the course, but I think a very important skill, 99 00:06:31,830 --> 00:06:35,310 is to be able to get like a chemical engineering problem 100 00:06:35,310 --> 00:06:39,630 and rearrange it so it becomes a standard form 101 00:06:39,630 --> 00:06:43,930 that maps up to one of the solvers you know how to use. 102 00:06:43,930 --> 00:06:46,260 So you've got a problem has some differential in it. 103 00:06:46,260 --> 00:06:48,450 Right away you actually, can I somehow rewrite it 104 00:06:48,450 --> 00:06:51,710 so it looks like this standard form-- or actually for Matlab, 105 00:06:51,710 --> 00:06:52,480 this form. 106 00:06:52,480 --> 00:06:58,800 So Matlab will say f of t, y. 107 00:06:58,800 --> 00:07:02,070 If you can do that, then you can go ahead and use those Matlab 108 00:07:02,070 --> 00:07:04,230 solvers to solve it. 109 00:07:04,230 --> 00:07:08,141 But your variables may not have a letter y appearing anywhere, 110 00:07:08,141 --> 00:07:08,640 right? 111 00:07:08,640 --> 00:07:09,570 You have to figure out how to write it 112 00:07:09,570 --> 00:07:10,980 into the standard form. 113 00:07:10,980 --> 00:07:14,190 And then once you can do that, you can use the solvers. 114 00:07:14,190 --> 00:07:16,110 And the same thing with the non-linear 115 00:07:16,110 --> 00:07:18,700 equations solvers, with the optimization software. 116 00:07:18,700 --> 00:07:20,580 If you can rewrite it into the form that 117 00:07:20,580 --> 00:07:23,310 matches the standard form, then once you're done, you're done. 118 00:07:23,310 --> 00:07:26,220 Then you can just go and call the canned programs 119 00:07:26,220 --> 00:07:27,780 that work well for those forms. 120 00:07:32,130 --> 00:07:34,160 All right, now what do we want as the output? 121 00:07:34,160 --> 00:07:36,900 So this is the problem. 122 00:07:36,900 --> 00:07:38,840 This is the problem as proposed. 123 00:07:38,840 --> 00:07:41,150 What do we really want out of this? 124 00:07:41,150 --> 00:07:47,190 What we're going to get out of it is a bunch of points y-- 125 00:07:47,190 --> 00:07:48,260 actually a matrix, right? 126 00:07:48,260 --> 00:07:50,240 It would be a vector of y-values. 127 00:07:50,240 --> 00:07:51,710 Typically y is a vector. 128 00:07:51,710 --> 00:07:55,220 And we want to know it a bunch of time points. 129 00:07:55,220 --> 00:07:56,700 And we would a lot of time points 130 00:07:56,700 --> 00:07:59,110 so we're getting pretty plots of how the components of y 131 00:07:59,110 --> 00:08:01,404 vary with time. 132 00:08:01,404 --> 00:08:02,570 That's what we usually want. 133 00:08:02,570 --> 00:08:07,480 Sometimes we just want the final value, y time value. 134 00:08:07,480 --> 00:08:09,770 We integrate to some point and then stop, and just 135 00:08:09,770 --> 00:08:10,550 run through the final value. 136 00:08:10,550 --> 00:08:12,440 But a lot of times, we actually want to make pretty plots, 137 00:08:12,440 --> 00:08:15,250 so we really want a lot of points here and making plots. 138 00:08:15,250 --> 00:08:16,910 Now how many points do you need to make a nice plot? 139 00:08:16,910 --> 00:08:18,701 Maybe, I don't know, 100 points, 50 points, 140 00:08:18,701 --> 00:08:20,000 you can make a nice plot. 141 00:08:20,000 --> 00:08:22,520 So that's the kind numbers you'd like. 142 00:08:22,520 --> 00:08:26,270 And so the output of these programs 143 00:08:26,270 --> 00:08:30,530 is going to be something like a vector of the time points, 144 00:08:30,530 --> 00:08:40,380 and then a matrix of the y-values corresponding 145 00:08:40,380 --> 00:08:41,750 to each time point. 146 00:08:41,750 --> 00:08:44,510 And the way it does in that lab, it's reasonable. 147 00:08:44,510 --> 00:08:47,870 Each row of y are the y-values that 148 00:08:47,870 --> 00:08:52,050 correspond to the same time point, the first time point. 149 00:08:52,050 --> 00:08:54,470 So the first row of y corresponds to the first number 150 00:08:54,470 --> 00:08:56,000 in the time vector. 151 00:08:56,000 --> 00:08:57,830 The second run is the second. 152 00:08:57,830 --> 00:09:01,400 Last one will be at time final. 153 00:09:01,400 --> 00:09:03,692 And then you can plot them. 154 00:09:03,692 --> 00:09:04,790 OK so far? 155 00:09:08,940 --> 00:09:16,940 And ideally, we would like these y-values that we compute, 156 00:09:16,940 --> 00:09:19,280 these y's we get out, ideally one 157 00:09:19,280 --> 00:09:21,740 that really close to what the true solution 158 00:09:21,740 --> 00:09:23,914 of the true differential equation is. 159 00:09:23,914 --> 00:09:25,830 But we're doing numerical stuff, so it's never 160 00:09:25,830 --> 00:09:27,810 going to be exactly the same. 161 00:09:27,810 --> 00:09:30,340 And so a big part of the effort is trying to figure out, 162 00:09:30,340 --> 00:09:33,155 how can we get the thing to compute y calculated 163 00:09:33,155 --> 00:09:34,530 at each time point to be as close 164 00:09:34,530 --> 00:09:36,750 as possible to the truth, of what the true solution 165 00:09:36,750 --> 00:09:39,390 of the equations are? 166 00:09:39,390 --> 00:09:41,580 And that turns out to be difficult, 167 00:09:41,580 --> 00:09:44,010 so that's a very important thing to worry about. 168 00:09:47,237 --> 00:09:48,570 What else can we say about this? 169 00:09:52,100 --> 00:09:57,280 We specify this from t naught, and you typically have to input 170 00:09:57,280 --> 00:10:00,610 your t final, how far you want to go away from your t naught 171 00:10:00,610 --> 00:10:04,000 as another input. 172 00:10:04,000 --> 00:10:07,870 Normally, people write it where t naught is less than t final, 173 00:10:07,870 --> 00:10:10,540 and so you're like integrating from left to right. 174 00:10:10,540 --> 00:10:12,970 But you could actually write them, put a negative sign 175 00:10:12,970 --> 00:10:14,950 equation which corresponds to a negative t, 176 00:10:14,950 --> 00:10:17,140 change your variable to of t prime that 177 00:10:17,140 --> 00:10:20,404 was equal to t final minus t if you wanted to, 178 00:10:20,404 --> 00:10:21,820 and integrate the other direction. 179 00:10:21,820 --> 00:10:26,230 And in fact in homework zero problem one, we did that too. 180 00:10:26,230 --> 00:10:27,580 So you can integrate frontwards. 181 00:10:27,580 --> 00:10:28,867 You can integrate backwards. 182 00:10:28,867 --> 00:10:31,200 If somebody tells you the initial condition or condition 183 00:10:31,200 --> 00:10:33,810 in the center of the range, you integrate forward 184 00:10:33,810 --> 00:10:36,270 for part of the range, and backwards the other way, 185 00:10:36,270 --> 00:10:37,224 so anything like this. 186 00:10:37,224 --> 00:10:38,640 So what's important for this to be 187 00:10:38,640 --> 00:10:40,140 what's called an initial value problem is just 188 00:10:40,140 --> 00:10:41,940 that all of the specifications of y 189 00:10:41,940 --> 00:10:43,780 are given at the same value of t. 190 00:10:43,780 --> 00:10:46,480 It doesn't matter exactly what value of t it is. 191 00:10:46,480 --> 00:10:47,980 Later, we'll talk about what happens 192 00:10:47,980 --> 00:10:50,005 if you don't know all the specifications at one 193 00:10:50,005 --> 00:10:52,470 value of time. 194 00:10:52,470 --> 00:11:01,250 All right, now how would you approach this? 195 00:11:01,250 --> 00:11:04,240 Probably, a bunch of you've done this already as undergraduates, 196 00:11:04,240 --> 00:11:05,530 or maybe even high school. 197 00:11:05,530 --> 00:11:08,690 You can write like a Taylor expansion, 198 00:11:08,690 --> 00:11:12,070 y of t plus delta t-- 199 00:11:12,070 --> 00:11:16,550 is y of t plus delta t-- 200 00:11:16,550 --> 00:11:20,261 times dy dt. 201 00:11:23,560 --> 00:11:25,570 You can write that. 202 00:11:25,570 --> 00:11:27,040 And then I say, wow, I know dy dt. 203 00:11:27,040 --> 00:11:30,045 That's actually f, so this is actually equal to-- 204 00:11:47,750 --> 00:11:48,910 right? 205 00:11:48,910 --> 00:11:51,140 And then if I truncate the Taylor expansion, 206 00:11:51,140 --> 00:11:52,390 now I have an update formula. 207 00:11:52,390 --> 00:11:55,600 So if I just forget this last bit, 208 00:11:55,600 --> 00:11:58,880 I can put in my current value of y, like the initial condition. 209 00:11:58,880 --> 00:12:01,000 The initial condition here, evaluate f, 210 00:12:01,000 --> 00:12:03,800 multiply by delta t, add them together, 211 00:12:03,800 --> 00:12:06,460 and I'll get a new value of y. 212 00:12:06,460 --> 00:12:09,300 And then I can repeat over and over again. 213 00:12:09,300 --> 00:12:11,120 And in fact, this is a very famous method. 214 00:12:11,120 --> 00:12:12,661 It's called the Forward Euler method. 215 00:12:12,661 --> 00:12:15,330 Euler was a pretty smart guy, so it's not completely ridiculous, 216 00:12:15,330 --> 00:12:17,121 but as I'll show you, it has some problems. 217 00:12:21,000 --> 00:12:25,850 So the Forward Euler, another way to write it is y new 218 00:12:25,850 --> 00:12:37,026 is equal to y old plus delta t times f of y old. 219 00:12:42,510 --> 00:12:44,010 It's called Forward Euler. 220 00:12:53,090 --> 00:13:01,030 And we can write a little implementation of this 221 00:13:01,030 --> 00:13:01,700 pretty easily. 222 00:13:01,700 --> 00:13:05,707 So you can write y is equal to y naught, 223 00:13:05,707 --> 00:13:13,625 t is equal to t naught for i equals 1 224 00:13:13,625 --> 00:13:22,540 to the number of steps, y is equal to y 225 00:13:22,540 --> 00:13:30,082 plus delta t times f of y, t is equal to t plus delta t. 226 00:13:30,082 --> 00:13:31,873 Somewhere up here, I need to write delta t. 227 00:13:41,540 --> 00:13:44,765 OK, so there's your code for doing Forward Euler. 228 00:13:44,765 --> 00:13:47,700 How many of you have done this before? 229 00:13:47,700 --> 00:13:49,630 OK, so I can skip over this very fast. 230 00:13:52,630 --> 00:13:55,236 How many did this, and it didn't work for a problem? 231 00:13:58,110 --> 00:14:00,170 OK, so you're not doing hard enough problems. 232 00:14:00,170 --> 00:14:03,710 Well, we'll correct that. 233 00:14:03,710 --> 00:14:06,800 All right, so I noticed that this closely 234 00:14:06,800 --> 00:14:11,350 resembles the rectangle rule. 235 00:14:11,350 --> 00:14:16,000 So the rectangle rule would look a lot the same. 236 00:14:16,000 --> 00:14:20,830 However, if you're integrating a function-- 237 00:14:20,830 --> 00:14:29,110 so I have I is equal to the integral of f of t dt, 238 00:14:29,110 --> 00:14:31,060 then I can notice that the derivative dI 239 00:14:31,060 --> 00:14:36,100 dt is equal to f of t. 240 00:14:39,900 --> 00:14:42,995 And so if I was going to do the rectangle rule, 241 00:14:42,995 --> 00:14:45,578 it would look just the same as this, except this would be a t. 242 00:14:48,330 --> 00:14:51,410 Maybe I would pick y0 to be 0. 243 00:14:51,410 --> 00:14:55,250 And the final value of y I would get to 244 00:14:55,250 --> 00:14:59,640 would be my integral, my i value that I want. 245 00:14:59,640 --> 00:15:02,490 So how many of you have done rectangle rule? 246 00:15:02,490 --> 00:15:04,020 Yes? 247 00:15:04,020 --> 00:15:06,600 All right, so you did this already. 248 00:15:06,600 --> 00:15:09,180 Now you remember in high school, when they took the rectangle 249 00:15:09,180 --> 00:15:14,639 rule, they probably mentioned that it's not very accurate. 250 00:15:14,639 --> 00:15:16,680 And so then they also taught you some other ones. 251 00:15:16,680 --> 00:15:19,496 You guys, how many of you learned the trapezoid rule? 252 00:15:19,496 --> 00:15:23,040 How about the midpoint rule? 253 00:15:23,040 --> 00:15:27,050 How about Simpson's rule? 254 00:15:27,050 --> 00:15:30,380 All right, so at least a high school math teachers 255 00:15:30,380 --> 00:15:36,661 thought that this is inadequate for doing real work, 256 00:15:36,661 --> 00:15:38,660 so they tell you all these other things instead. 257 00:15:38,660 --> 00:15:42,020 This was like your first baby thing. 258 00:15:42,020 --> 00:15:43,760 So let's talk about why was that. 259 00:15:59,780 --> 00:16:03,870 So let's do the Taylor expansion a little bit further. 260 00:16:03,870 --> 00:16:08,650 So let's do it for the numeric [INAUDIBLE] case. 261 00:16:08,650 --> 00:16:17,670 So y t plus delta t is equal to y of t plus delta t times 262 00:16:17,670 --> 00:16:26,910 f of t plus 1/2 delta t squared times the second derivative. 263 00:16:32,320 --> 00:16:34,510 And so when we made the approximation 264 00:16:34,510 --> 00:16:37,000 to do the rectangle rule, we just threw this term away. 265 00:16:37,000 --> 00:16:39,400 So that's a leading term that we threw away. 266 00:16:39,400 --> 00:16:43,280 So we made an error order or delta t squared. 267 00:16:43,280 --> 00:16:46,330 Because normally the second derivative is not exactly zero. 268 00:16:46,330 --> 00:16:57,330 So we have an error in each step that's 269 00:16:57,330 --> 00:17:00,941 the order of delta t squared. 270 00:17:00,941 --> 00:17:02,190 But how many steps do we make? 271 00:17:02,190 --> 00:17:07,225 The number of steps is equal to t final minus t 272 00:17:07,225 --> 00:17:08,550 naught over delta t. 273 00:17:11,270 --> 00:17:16,660 So therefore, we're making one over delta t order steps. 274 00:17:16,660 --> 00:17:19,089 So the total error is going to be approximately the number 275 00:17:19,089 --> 00:17:21,700 of steps times how much error you making each step. 276 00:17:21,700 --> 00:17:24,099 So it's going to be something the second power in delta t 277 00:17:24,099 --> 00:17:26,223 divided by something to the first power of delta t. 278 00:17:26,223 --> 00:17:33,730 So the total error in the integral I 279 00:17:33,730 --> 00:17:39,080 is going to be older of delta t if you do the rectangle rule. 280 00:17:39,080 --> 00:17:41,730 And so order delta t is not very good. 281 00:17:41,730 --> 00:17:45,870 So you want to increase your accuracy-- 282 00:17:45,870 --> 00:17:47,550 if you cut your step size in half, 283 00:17:47,550 --> 00:17:49,050 that means you double your CPU time, 284 00:17:49,050 --> 00:17:51,870 because you have to do twice as many function evaluations. 285 00:17:51,870 --> 00:17:55,710 But you only gain-- you reduce your error by a factor of two. 286 00:17:55,710 --> 00:17:59,510 So if you want to, say, gain three significant figures 287 00:17:59,510 --> 00:18:01,200 of accuracy in your number, then you 288 00:18:01,200 --> 00:18:06,770 have to do 1,000 times much work as you did before. 289 00:18:06,770 --> 00:18:08,640 So I have some amount of some precision 290 00:18:08,640 --> 00:18:09,979 with some number of steps. 291 00:18:09,979 --> 00:18:12,270 I want to increase, get three more significant figures. 292 00:18:12,270 --> 00:18:14,520 I'll have to do 1,000 times as much work. 293 00:18:14,520 --> 00:18:16,504 That's kind of bad scaling. 294 00:18:16,504 --> 00:18:18,420 If I want to get six more significant figures, 295 00:18:18,420 --> 00:18:20,630 I've got to do a million times more work. 296 00:18:20,630 --> 00:18:25,000 So at some point, this going to be kind of expensive. 297 00:18:25,000 --> 00:18:31,730 So if you contrast it, you can do things like trapezoid rule. 298 00:18:31,730 --> 00:18:35,580 And as you're probably-- 299 00:18:35,580 --> 00:18:38,193 well let's talk about trapezoid rule for a bit. 300 00:18:48,030 --> 00:18:52,140 So the idea of this is you have your function, 301 00:18:52,140 --> 00:18:56,340 you have some point t, and t plus delta t. 302 00:18:56,340 --> 00:18:59,420 You're trying to integrate between them. 303 00:18:59,420 --> 00:19:03,384 Let's say t naught, t naught plus delta t. 304 00:19:03,384 --> 00:19:04,960 Try to integrate. 305 00:19:04,960 --> 00:19:07,510 Here's the value of the function at t. 306 00:19:07,510 --> 00:19:13,301 Here is the value of the function that t plus delta t. 307 00:19:13,301 --> 00:19:15,050 If you do the rectangle rule, what you say 308 00:19:15,050 --> 00:19:20,125 is I just draw a line here and integrate that rectangle. 309 00:19:20,125 --> 00:19:24,480 But say the real function really looks like this. 310 00:19:24,480 --> 00:19:27,560 Then I missing that area there, so that's the error. 311 00:19:27,560 --> 00:19:30,280 That's why the error's so big. 312 00:19:30,280 --> 00:19:37,140 If I do that trapezoid rule, then the error-- 313 00:19:37,140 --> 00:19:40,570 I actually overestimate the value here, 314 00:19:40,570 --> 00:19:42,130 but the integral of the difference 315 00:19:42,130 --> 00:19:43,930 between this dotted line and the solid line 316 00:19:43,930 --> 00:19:45,990 is less than the integral between the solid line 317 00:19:45,990 --> 00:19:48,637 and that dotted line, so its more accurate. 318 00:19:48,637 --> 00:19:50,470 So that's why people like the trapezoid rule 319 00:19:50,470 --> 00:19:52,886 better than the rectangle rule. 320 00:19:52,886 --> 00:19:55,270 If you did the midpoint rule, you'd 321 00:19:55,270 --> 00:20:00,320 evaluate this function halfway between these two points 322 00:20:00,320 --> 00:20:02,815 and then draw a dotted line here and a dotted line there, 323 00:20:02,815 --> 00:20:05,020 and integrate that guy. 324 00:20:05,020 --> 00:20:07,144 And it would overestimate for part, 325 00:20:07,144 --> 00:20:09,310 and underestimate for part, and they would partially 326 00:20:09,310 --> 00:20:10,940 cancel each other out. 327 00:20:10,940 --> 00:20:13,589 And so that also would be more accurate than doing 328 00:20:13,589 --> 00:20:14,380 the rectangle rule. 329 00:20:17,745 --> 00:20:18,620 So you did it before. 330 00:20:23,740 --> 00:20:25,660 And so what we're going to try to do 331 00:20:25,660 --> 00:20:27,430 is, instead of using this formula 332 00:20:27,430 --> 00:20:31,870 we did before, we're going to replace this with something 333 00:20:31,870 --> 00:20:35,300 else that I'm going to call g. 334 00:20:35,300 --> 00:20:37,040 And g is something that's supposed 335 00:20:37,040 --> 00:20:41,510 to be the average value of f over the step. 336 00:20:41,510 --> 00:20:45,410 Instead of just taking the value f at the starting point, 337 00:20:45,410 --> 00:20:47,239 I want to get sort of the average. 338 00:20:47,239 --> 00:20:48,780 Right, because an integral is related 339 00:20:48,780 --> 00:20:51,170 to an average over the delta t. 340 00:20:51,170 --> 00:20:54,500 So g is going to be my way to estimate the average. 341 00:20:54,500 --> 00:20:56,429 And there's going to be a zillion update 342 00:20:56,429 --> 00:20:58,220 formulas that different mathematicians have 343 00:20:58,220 --> 00:20:59,626 proposed over the years. 344 00:20:59,626 --> 00:21:02,250 They give you different g's that tell you how to do the update. 345 00:21:05,240 --> 00:21:07,500 And the key is try to find a g that's really accurate. 346 00:21:07,500 --> 00:21:09,420 You want to have one that's the most accurate you can, 347 00:21:09,420 --> 00:21:11,060 because then your errors will be less. 348 00:21:11,060 --> 00:21:13,910 And if your errors are less, than your delta t can 349 00:21:13,910 --> 00:21:17,250 be bigger, which means your number of steps will be less, 350 00:21:17,250 --> 00:21:19,700 which means your CPU time will be less. 351 00:21:19,700 --> 00:21:21,920 So you might be able to get more accuracy 352 00:21:21,920 --> 00:21:24,367 and use less CPU time both if you 353 00:21:24,367 --> 00:21:26,575 can find a really good formula for g, for the update. 354 00:21:29,830 --> 00:21:32,670 So what's the g for the trapezoid rule? 355 00:21:37,300 --> 00:21:46,190 It's f of t plus f of t plus delta t over 2. 356 00:21:46,190 --> 00:21:48,690 That was the g you used when you did the trapezoid rule when 357 00:21:48,690 --> 00:21:51,440 you were kids. 358 00:21:51,440 --> 00:21:53,010 What is the g for the midpoint rule? 359 00:21:57,280 --> 00:22:02,382 It's f of t plus delta t over 2. 360 00:22:06,410 --> 00:22:08,640 And so the different g's are called different update 361 00:22:08,640 --> 00:22:10,960 formulas. 362 00:22:10,960 --> 00:22:14,740 And these particular ones reduce the error 363 00:22:14,740 --> 00:22:20,770 from delta t squared for [INAUDIBLE] step to, I think, 364 00:22:20,770 --> 00:22:23,260 delta t to the cubed power. 365 00:22:23,260 --> 00:22:27,670 And then when you multiply by this 1 every delta t, 366 00:22:27,670 --> 00:22:30,980 it turns out you go from being order of delta t 367 00:22:30,980 --> 00:22:35,840 to order of delta t squared, and that's a really big change. 368 00:22:35,840 --> 00:22:40,280 So now if I cut the step size by a factor of 10, 369 00:22:40,280 --> 00:22:43,904 I gain a factor of 100 in the accuracy. 370 00:22:43,904 --> 00:22:45,820 If I want to get six more significant figures, 371 00:22:45,820 --> 00:22:48,910 I do 1,000 times more work, not a million times more work. 372 00:22:48,910 --> 00:22:52,220 So it's a pretty significant improvement. 373 00:22:52,220 --> 00:22:54,640 And then people have pushed this on further and further. 374 00:22:54,640 --> 00:22:56,320 And so actually very common integrators 375 00:22:56,320 --> 00:22:58,270 that you might use nowadays would go out 376 00:22:58,270 --> 00:23:02,160 to fourth and fifth order methods. 377 00:23:02,160 --> 00:23:04,242 And so they have complicated update formulas 378 00:23:04,242 --> 00:23:05,950 that are carefully designed to cancel out 379 00:23:05,950 --> 00:23:09,280 all the low order delta t errors and just leave very high order 380 00:23:09,280 --> 00:23:10,390 ones. 381 00:23:10,390 --> 00:23:11,906 And that's kind of the-- 382 00:23:11,906 --> 00:23:13,030 typical ones are like that. 383 00:23:13,030 --> 00:23:15,446 And then some fancy ones might go even up to eighth order, 384 00:23:15,446 --> 00:23:18,337 I've seen, if you're really pushing it. 385 00:23:18,337 --> 00:23:20,170 You have really complicated g formulas then. 386 00:23:27,260 --> 00:23:31,980 Now, this is all for the integration, 387 00:23:31,980 --> 00:23:35,950 but really our problem was that the ODEs [INAUDIBLE].. 388 00:23:35,950 --> 00:23:41,780 So we have to do f of y not f of t. 389 00:23:41,780 --> 00:23:45,070 Let's see if I can show you that. 390 00:23:45,070 --> 00:23:52,140 Yeah, so we did the rectangle rule, we had f of t. 391 00:23:52,140 --> 00:23:55,320 But the real problem is we have to do f of y. 392 00:23:55,320 --> 00:23:57,690 Now the problem is, we don't know 393 00:23:57,690 --> 00:24:00,710 what y is, y is our unknown. 394 00:24:00,710 --> 00:24:03,880 So we generally have a problem here. 395 00:24:03,880 --> 00:24:06,530 In Forward Euler, we can get away with it because we just 396 00:24:06,530 --> 00:24:08,930 put the y old value in. 397 00:24:08,930 --> 00:24:14,450 Like for example, if we wanted to do midpoint here, 398 00:24:14,450 --> 00:24:17,690 we'd have to know f at a different time point 399 00:24:17,690 --> 00:24:20,260 that we haven't computed yet, because it's forward in time. 400 00:24:20,260 --> 00:24:21,635 Same thing with trapezoidal rule. 401 00:24:21,635 --> 00:24:24,200 We have to know f at some future timepoint. 402 00:24:24,200 --> 00:24:26,540 So it's not so easy, actually, to go directly 403 00:24:26,540 --> 00:24:30,860 from these formulas to some explicit formula like this 404 00:24:30,860 --> 00:24:33,810 with a g in here, where the things inside g, 405 00:24:33,810 --> 00:24:37,810 the y's inside g, are all y-values we actually know. 406 00:24:37,810 --> 00:24:39,729 So this is a serious problem. 407 00:24:45,600 --> 00:24:48,570 But this didn't stop people from doing it. 408 00:24:48,570 --> 00:24:54,240 So one idea is you could go back to the Taylor expansion, 409 00:24:54,240 --> 00:25:01,450 and now do the Taylor expansion but in terms of f of y. 410 00:25:01,450 --> 00:25:03,700 So this is f of y of t. 411 00:25:06,240 --> 00:25:12,140 And then over here, this is d squared y dt squared, 412 00:25:12,140 --> 00:25:20,470 but dy dt is f, so this is really df dt. 413 00:25:20,470 --> 00:25:22,970 But f is actually not a function of t, it's a function of y, 414 00:25:22,970 --> 00:25:23,540 right? 415 00:25:23,540 --> 00:25:24,980 But we can do chain rule. 416 00:25:24,980 --> 00:25:41,440 So df dt is equal to the sum df dyn dyn dt. 417 00:25:45,630 --> 00:25:48,680 But dyn dt is f. 418 00:25:48,680 --> 00:25:50,960 And this is what we call the Jacobian [INAUDIBLE].. 419 00:25:50,960 --> 00:25:56,558 So this is really saying that this is equal to J times f. 420 00:25:59,550 --> 00:26:00,930 So that would be one possibility, 421 00:26:00,930 --> 00:26:02,846 is you could go and write the Taylor expansion 422 00:26:02,846 --> 00:26:05,190 out to the next higher order explicitly, 423 00:26:05,190 --> 00:26:08,580 putting in J times f here, and then just 424 00:26:08,580 --> 00:26:11,920 being in the cubic terms. 425 00:26:11,920 --> 00:26:13,370 Now, people don't do this usually. 426 00:26:13,370 --> 00:26:15,720 Any idea why this is not a popular thing to do? 427 00:26:18,996 --> 00:26:19,496 Yeah? 428 00:26:19,496 --> 00:26:22,059 AUDIENCE: It's computationally expensive to [INAUDIBLE].. 429 00:26:22,059 --> 00:26:23,350 WILLIAM GREEN JR: That's right. 430 00:26:23,350 --> 00:26:23,891 That's right. 431 00:26:23,891 --> 00:26:25,570 So how many function evaluations does it 432 00:26:25,570 --> 00:26:26,730 take to get the Jacobian? 433 00:26:32,270 --> 00:26:37,390 Right, so each function is n values, right, 434 00:26:37,390 --> 00:26:38,930 because it's a vector. 435 00:26:38,930 --> 00:26:43,530 And you have to do at least, say, forward differences 436 00:26:43,530 --> 00:26:46,700 would be n function evaluations. 437 00:26:46,700 --> 00:26:48,700 So instead of just doing one function evaluation 438 00:26:48,700 --> 00:26:50,710 like we were doing with the Forward Euler, now 439 00:26:50,710 --> 00:26:52,668 we have to do n plus 1, or something like that, 440 00:26:52,668 --> 00:26:54,190 function evaluations. 441 00:26:54,190 --> 00:26:56,606 And actually, you probably want to use center differences, 442 00:26:56,606 --> 00:26:59,260 so you need 2 times n plus 1. 443 00:26:59,260 --> 00:27:02,421 And if you do them analytically, unless it's sparse, 444 00:27:02,421 --> 00:27:04,670 it's still going to be really expensive to compute it. 445 00:27:04,670 --> 00:27:06,667 So it's a lot more effort to do it. 446 00:27:06,667 --> 00:27:08,750 So you need to think, well, is it really worth it? 447 00:27:08,750 --> 00:27:11,410 And so people have found other formulas that basically 448 00:27:11,410 --> 00:27:13,510 get rid of the delta t squared term that 449 00:27:13,510 --> 00:27:15,710 don't require so much effort, and so that's 450 00:27:15,710 --> 00:27:17,224 what people do instead. 451 00:27:17,224 --> 00:27:18,140 But this is an option. 452 00:27:18,140 --> 00:27:18,806 You could do it. 453 00:27:25,580 --> 00:27:27,620 All right, so instead, what people 454 00:27:27,620 --> 00:27:32,810 do is they've decided to generalize the midpoint rule 455 00:27:32,810 --> 00:27:34,790 and try to figure out, how can we 456 00:27:34,790 --> 00:27:39,380 get an estimate of the midpoint value that's cheap to evaluate? 457 00:27:39,380 --> 00:27:40,970 And so what they say is, well, this 458 00:27:40,970 --> 00:27:44,356 is g midpoint for numerical integration. 459 00:27:44,356 --> 00:27:45,980 Let's make an new g midpoint that works 460 00:27:45,980 --> 00:27:48,480 for when f is a function of y. 461 00:27:48,480 --> 00:27:52,700 So let's make it function of y evaluated 462 00:27:52,700 --> 00:27:53,930 at t plus delta t over 2. 463 00:27:59,787 --> 00:28:01,370 So that's what the formula looks like. 464 00:28:01,370 --> 00:28:02,870 But then you say, well, I don't know 465 00:28:02,870 --> 00:28:05,510 y is at t plus delta t over 2, because I only 466 00:28:05,510 --> 00:28:07,360 know it at y of t. 467 00:28:07,360 --> 00:28:10,140 So then they say, well, let's do Forward Euler for that. 468 00:28:10,140 --> 00:28:14,720 So then they say, well, let's say it's approximately f of-- 469 00:28:18,580 --> 00:28:20,380 what's the Forward Euler formula? 470 00:28:20,380 --> 00:28:28,442 It's y of t plus delta t over 2 times f of y of t. 471 00:28:31,154 --> 00:28:33,018 How many parenthesis there? 472 00:28:38,410 --> 00:28:43,680 So this is the formula that's actually 473 00:28:43,680 --> 00:28:45,360 used in practice a lot. 474 00:28:45,360 --> 00:28:48,120 And this one is called Runge-Kutta two, second order 475 00:28:48,120 --> 00:28:50,030 Runge-Kutta formula. 476 00:28:50,030 --> 00:28:52,080 And it's just a-- 477 00:28:52,080 --> 00:28:56,400 it's the midpoint rule where we estimate the value of y of t 478 00:28:56,400 --> 00:28:59,074 plus delta t using Forward Euler. 479 00:28:59,074 --> 00:29:01,240 And the advantage of this is it's very cheap, right? 480 00:29:01,240 --> 00:29:05,152 Just evaluate-- I just have to evaluate one function here, 481 00:29:05,152 --> 00:29:06,860 one function here, so two function calls. 482 00:29:10,269 --> 00:29:16,860 [INAUDIBLE] and that way I can get an update formula. 483 00:29:16,860 --> 00:29:19,810 And this one, turns out that it has the nice delta t 484 00:29:19,810 --> 00:29:21,360 scaling that you might like. 485 00:29:24,640 --> 00:29:27,950 So that's the kind of thing people do. 486 00:29:27,950 --> 00:29:30,600 And the Runge-Kutta guys went crazy, 487 00:29:30,600 --> 00:29:33,260 and so one of the most popular solvers 488 00:29:33,260 --> 00:29:36,920 is called Runge-Kutta four five, and that's the ODE four five 489 00:29:36,920 --> 00:29:38,690 program in Matlab. 490 00:29:38,690 --> 00:29:41,257 And what that does is the fourth order extension of this, 491 00:29:41,257 --> 00:29:42,840 and the fifth order extension of this. 492 00:29:42,840 --> 00:29:45,050 And it uses them both, and it compares them, 493 00:29:45,050 --> 00:29:47,900 and uses the difference between an estimate the error 494 00:29:47,900 --> 00:29:48,860 in the calculation. 495 00:29:48,860 --> 00:29:52,434 And then uses that to decide what size delta t you need, 496 00:29:52,434 --> 00:29:54,100 depending on what tolerances you demand. 497 00:29:58,016 --> 00:29:59,352 And the formulas for all those-- 498 00:29:59,352 --> 00:30:01,310 [INAUDIBLE] formulas are given in the textbook. 499 00:30:06,810 --> 00:30:09,960 Now, this starts to lead the problem. 500 00:30:09,960 --> 00:30:14,230 So we weren't able to actually evaluate the true g midpoint, 501 00:30:14,230 --> 00:30:15,400 which we'd like to be here. 502 00:30:15,400 --> 00:30:18,690 Instead, we have to use some approximation in order 503 00:30:18,690 --> 00:30:21,480 to extrapolate forwards to the y. 504 00:30:21,480 --> 00:30:24,230 So this is hinting at the whole problem. 505 00:30:24,230 --> 00:30:26,610 This whole problem is actually we're just extrapolating. 506 00:30:26,610 --> 00:30:31,260 We know the true value of y only at t naught, 507 00:30:31,260 --> 00:30:34,950 and all the rest we're doing is extrapolating forwards. 508 00:30:34,950 --> 00:30:36,450 And we're do it over and over again, 509 00:30:36,450 --> 00:30:37,949 so we're going to do it for n steps. 510 00:30:37,949 --> 00:30:40,570 Maybe we'll do 1,000 steps. 511 00:30:40,570 --> 00:30:41,769 So you're extrapolating. 512 00:30:41,769 --> 00:30:43,310 And then based on that extrapolation, 513 00:30:43,310 --> 00:30:44,760 you going to extrapolate again. 514 00:30:44,760 --> 00:30:47,610 And then based on extrapolation, you going to extrapolate again. 515 00:30:47,610 --> 00:30:50,520 And you can see that this is not really the most robust thing 516 00:30:50,520 --> 00:30:52,050 to do, right? 517 00:30:52,050 --> 00:30:54,090 Because if you have any experience, 518 00:30:54,090 --> 00:30:57,480 you're not that comfortable extrapolating at all. 519 00:30:57,480 --> 00:31:00,350 And then you have to extrapolate 1,000 times. 520 00:31:00,350 --> 00:31:03,220 You should be a little bit worried about what can happen. 521 00:31:03,220 --> 00:31:05,500 So let's think about how the error will propagate. 522 00:31:05,500 --> 00:31:09,180 So we're going to have some errors while we 523 00:31:09,180 --> 00:31:10,421 do these extrapolations. 524 00:31:18,780 --> 00:31:23,520 So we have here a t naught, t naught, we're happy. 525 00:31:23,520 --> 00:31:24,810 We know y naught. 526 00:31:24,810 --> 00:31:26,110 This is exact. 527 00:31:26,110 --> 00:31:27,570 We're like woohoo. 528 00:31:27,570 --> 00:31:29,940 We know one value of y, really good. 529 00:31:29,940 --> 00:31:34,716 All right, so now the question is the first step. 530 00:31:34,716 --> 00:31:36,360 Now we're going to get to value y1 531 00:31:36,360 --> 00:31:38,510 that we compute with one of the formulas, 532 00:31:38,510 --> 00:31:40,260 and it's going to be good to some accuracy 533 00:31:40,260 --> 00:31:43,470 depending on how good or g is, our update formula. 534 00:31:43,470 --> 00:31:44,930 And so it's going to have an error. 535 00:31:47,750 --> 00:31:50,290 I should comment that this y naught, 536 00:31:50,290 --> 00:31:52,980 although we treat it as exact, it could have an error too. 537 00:31:52,980 --> 00:31:55,188 But it will probably be more like a machine precision 538 00:31:55,188 --> 00:31:57,010 kind of error because we don't really 539 00:31:57,010 --> 00:31:58,080 know the initial conditions that well. 540 00:31:58,080 --> 00:31:59,290 Or might be an experimental error 541 00:31:59,290 --> 00:32:00,940 about how well we know what the initial conditions are. 542 00:32:00,940 --> 00:32:02,590 But I'm not going to worry about that one for now. 543 00:32:02,590 --> 00:32:04,600 But be aware, this is true, we don't really 544 00:32:04,600 --> 00:32:07,700 know initial conditions perfectly either. 545 00:32:07,700 --> 00:32:11,075 But certainly, even if we treat this as being known perfectly, 546 00:32:11,075 --> 00:32:13,450 we're going to have some error by the time we extrapolate 547 00:32:13,450 --> 00:32:15,660 to y1. 548 00:32:15,660 --> 00:32:17,660 So now what happens at y2? 549 00:32:17,660 --> 00:32:21,780 So now we've added some more, so now we're at t2. 550 00:32:21,780 --> 00:32:28,350 And we see started from y1 true plus some error. 551 00:32:28,350 --> 00:32:29,760 This is y-- y naught was true. 552 00:32:32,600 --> 00:32:35,520 Now we're starting from y1 true with some error in it, 553 00:32:35,520 --> 00:32:38,290 because of a mistake we made the first step. 554 00:32:38,290 --> 00:32:42,240 And now we see how that step propagates further. 555 00:32:42,240 --> 00:32:45,170 So let's see what that's going to say. 556 00:32:45,170 --> 00:32:53,400 Well say y2 is going to equal y1 plus delta t times 557 00:32:53,400 --> 00:32:55,710 some update formula which would depend on y1. 558 00:33:00,770 --> 00:33:03,820 But y1 was not true anymore, so this is actually 559 00:33:03,820 --> 00:33:05,990 equal to what should have been the true value of y1 560 00:33:05,990 --> 00:33:09,176 if we only knew what it was, but we don't know. 561 00:33:09,176 --> 00:33:11,480 And we have our error that we don't know exactly how 562 00:33:11,480 --> 00:33:14,600 big it is that we added on. 563 00:33:14,600 --> 00:33:18,980 And then we have another term from this guy, 564 00:33:18,980 --> 00:33:22,370 so it's going to be delta t times 565 00:33:22,370 --> 00:33:28,520 g of y1 true plus delta 1. 566 00:33:31,240 --> 00:33:34,240 And then we know our update formula's not right anyway, 567 00:33:34,240 --> 00:33:37,621 so there's actually another error, delta 2, 568 00:33:37,621 --> 00:33:39,621 just because of the error in the update formula. 569 00:33:45,030 --> 00:33:48,089 Now if we had a great update formula, 570 00:33:48,089 --> 00:33:49,380 these deltas are kind of small. 571 00:33:52,240 --> 00:33:57,844 And if somehow, by luck, our delta 1 was 0, 572 00:33:57,844 --> 00:33:59,260 and we had a great update formula, 573 00:33:59,260 --> 00:34:01,926 then we think that this error in y2 is going to be pretty small. 574 00:34:04,440 --> 00:34:05,580 But it's not true. 575 00:34:05,580 --> 00:34:07,080 So really the errors are pretty big. 576 00:34:10,880 --> 00:34:13,790 So let's see if we can figure out-- 577 00:34:13,790 --> 00:34:17,010 we're going to do a Taylor expansion here, 578 00:34:17,010 --> 00:34:27,139 so this'll be g of y1 true plus d dgy-- 579 00:34:27,139 --> 00:34:29,800 the Jacobian of g with respect to y-- 580 00:34:29,800 --> 00:34:33,199 times the delta 1. 581 00:34:33,199 --> 00:34:36,884 That's the first order Taylor extension. 582 00:34:36,884 --> 00:34:38,550 All right, so now let's group the terms. 583 00:34:48,679 --> 00:34:58,580 So y2 is equal to y1 true plus the update formula of y1 584 00:34:58,580 --> 00:35:03,710 true plus the error in our update formula 585 00:35:03,710 --> 00:35:08,310 if we had started at the true value. 586 00:35:08,310 --> 00:35:09,780 So this is what we-- 587 00:35:09,780 --> 00:35:12,320 if somebody told us what the true value of y1, this 588 00:35:12,320 --> 00:35:14,020 is the size error we'd expect. 589 00:35:14,020 --> 00:35:17,390 It's to be that little error we calculated before. 590 00:35:17,390 --> 00:35:20,600 But then on top of it, we have these other terms. 591 00:35:20,600 --> 00:35:25,490 We have a delta 1, because y1 was not the true value. 592 00:35:25,490 --> 00:35:26,780 Sorry, I lost a delta t here. 593 00:35:31,490 --> 00:35:46,080 And then we have delta t times dg dy times delta 1. 594 00:35:46,080 --> 00:35:47,860 And I'll emphasize these are all vectors. 595 00:35:53,283 --> 00:35:56,670 I think that's it if we only keep to first order. 596 00:36:00,140 --> 00:36:01,970 So what is this thing? 597 00:36:01,970 --> 00:36:06,833 This is a matrix, right, dg dy. 598 00:36:06,833 --> 00:36:09,600 I could write it this way. 599 00:36:09,600 --> 00:36:13,980 It's the identity matrix plus delta t times the Jacobian dg 600 00:36:13,980 --> 00:36:18,500 dy times the vector d1. 601 00:36:21,990 --> 00:36:24,430 And you can see, at every iteration when I do this, 602 00:36:24,430 --> 00:36:27,500 I'm going to get a similar factor like this. 603 00:36:27,500 --> 00:36:30,820 And it's going to again multiply delta 1. 604 00:36:30,820 --> 00:36:32,385 I keep on multiplying the-- 605 00:36:32,385 --> 00:36:33,760 errors from my previous step keep 606 00:36:33,760 --> 00:36:37,090 getting multiplied by this kind of factor. 607 00:36:37,090 --> 00:36:38,840 So it's going to be really important to us 608 00:36:38,840 --> 00:36:40,540 about what the norm of this matrix is. 609 00:36:52,641 --> 00:36:55,099 What's going to happen if the normal of this matrix is big? 610 00:36:58,452 --> 00:36:59,410 What's going to happen? 611 00:36:59,410 --> 00:37:00,534 AUDIENCE: [INAUDIBLE] 612 00:37:00,534 --> 00:37:02,575 WILLIAM GREEN JR: That's right, the error's going 613 00:37:02,575 --> 00:37:04,210 to get multiplied by some factor, 614 00:37:04,210 --> 00:37:06,452 say bigger than 1, every iteration. 615 00:37:06,452 --> 00:37:08,160 And after we go 1,000 iterations, how big 616 00:37:08,160 --> 00:37:11,770 is the error going to be? 617 00:37:11,770 --> 00:37:12,640 Really big, right? 618 00:37:12,640 --> 00:37:16,060 Even if that thing is quite close to 1, 1 plus something 619 00:37:16,060 --> 00:37:18,760 to the 1,000 power it's gigantic. 620 00:37:18,760 --> 00:37:21,010 So what was originally a pretty small error-- 621 00:37:21,010 --> 00:37:22,880 because we chose a good g method-- 622 00:37:22,880 --> 00:37:26,020 is going to get multiplied by this huge amplification factor. 623 00:37:26,020 --> 00:37:27,850 So the key thing for a method to be 624 00:37:27,850 --> 00:37:30,190 what is called numerically stable 625 00:37:30,190 --> 00:37:35,280 is that the norm of this matrix has got to be less than 1. 626 00:37:35,280 --> 00:37:37,700 Or maybe even equal, equals 1 might be OK. 627 00:37:37,700 --> 00:37:40,860 If this is much bigger than 1, you're doomed. 628 00:37:40,860 --> 00:37:43,390 Your numeric order is just going to blow up 629 00:37:43,390 --> 00:37:46,094 from step to step to step. 630 00:37:46,094 --> 00:37:48,010 And what's bad about it, is actually sometimes 631 00:37:48,010 --> 00:37:50,361 you might not even be aware of this. 632 00:37:50,361 --> 00:37:52,860 So you'll get some solution, because it still will calculate 633 00:37:52,860 --> 00:37:57,070 some numbers for the y's at each type-- y1, y2, y3, y4-- 634 00:37:57,070 --> 00:38:00,740 it's just running through the calculation. 635 00:38:00,740 --> 00:38:04,070 But the numbers may become increasingly meaningless 636 00:38:04,070 --> 00:38:07,970 further away from the y true as time goes on, as the steps go. 637 00:38:14,275 --> 00:38:20,100 AUDIENCE: [INAUDIBLE] 638 00:38:20,100 --> 00:38:21,469 WILLIAM GREEN JR: Yes, it can. 639 00:38:21,469 --> 00:38:24,010 Yeah, so that's what we call a stable numerical method is one 640 00:38:24,010 --> 00:38:26,440 that, if you make an error in an earlier step, 641 00:38:26,440 --> 00:38:29,980 its impact diminishes as the steps go on. 642 00:38:29,980 --> 00:38:31,330 That's what you want. 643 00:38:31,330 --> 00:38:32,580 You can't always achieve this. 644 00:38:32,580 --> 00:38:35,490 But if you can, that's really great. 645 00:38:35,490 --> 00:38:37,930 Then that's you would call a really stable, robust kind 646 00:38:37,930 --> 00:38:39,080 of method. 647 00:38:39,080 --> 00:38:42,070 Even though I give you some stupidity early on, 648 00:38:42,070 --> 00:38:43,570 my stupidity evaporates away, and it 649 00:38:43,570 --> 00:38:45,455 goes to the true solution. 650 00:38:45,455 --> 00:38:46,330 That's what you want. 651 00:38:48,847 --> 00:38:50,680 OK, so let's think of cases where it's going 652 00:38:50,680 --> 00:38:52,660 to be problematic for sure. 653 00:38:52,660 --> 00:38:56,680 So if dg dy-- 654 00:38:56,680 --> 00:39:00,860 say if all it's eigenvalues are positive 655 00:39:00,860 --> 00:39:03,540 and you add it to identity matrix, 656 00:39:03,540 --> 00:39:06,890 the thing is still going to have eigenvalues bigger than 1. 657 00:39:06,890 --> 00:39:09,350 Can you see that? 658 00:39:09,350 --> 00:39:12,357 And actually, what do we say when dg dy is positive, 659 00:39:12,357 --> 00:39:13,440 has a positive eigenvalue? 660 00:39:13,440 --> 00:39:14,600 You guys have a word for this already 661 00:39:14,600 --> 00:39:16,558 you learn in your differential equations class. 662 00:39:16,558 --> 00:39:20,030 What do you call these kind of differential equations? 663 00:39:20,030 --> 00:39:20,870 Maybe not dg dy. 664 00:39:20,870 --> 00:39:22,940 Back up, how about df dy? 665 00:39:22,940 --> 00:39:26,270 If df dy, if that Jacobian has a positive eigenvalue, 666 00:39:26,270 --> 00:39:27,580 what do you say? 667 00:39:34,230 --> 00:39:35,710 You guys remember this? 668 00:39:35,710 --> 00:39:39,700 So this is what we call unstable differential equations. 669 00:39:39,700 --> 00:39:42,970 So if it has any positive eigenvalues, 670 00:39:42,970 --> 00:39:49,900 then two initial conditions that differ by a little differential 671 00:39:49,900 --> 00:39:53,650 exponentially separate as time goes on as long 672 00:39:53,650 --> 00:39:56,050 as the difference as any projection in the direction 673 00:39:56,050 --> 00:39:59,630 of the bad eigenvector. 674 00:39:59,630 --> 00:40:01,970 So those are called numerically unstable. 675 00:40:01,970 --> 00:40:04,980 A lot of those actually exist in reality. 676 00:40:04,980 --> 00:40:08,030 So explosions, that's because you 677 00:40:08,030 --> 00:40:12,200 had some positive eigenvalue somewhere for example. 678 00:40:12,200 --> 00:40:14,090 Amplifiers, like you buy an amplifier 679 00:40:14,090 --> 00:40:17,900 to crank up the music, you get positive feedback when-- 680 00:40:17,900 --> 00:40:20,360 remember Professor Swan walked to some strange place, 681 00:40:20,360 --> 00:40:22,967 and for some reason, he got the feedback from the speakers? 682 00:40:22,967 --> 00:40:24,800 So he had a positive eigenvalue going there. 683 00:40:24,800 --> 00:40:26,716 A little tiny noise and his microphone, and it 684 00:40:26,716 --> 00:40:28,700 amplified to a big squeak. 685 00:40:28,700 --> 00:40:31,790 So there are real systems like this that amplify. 686 00:40:31,790 --> 00:40:33,332 And sometimes we want to model them, 687 00:40:33,332 --> 00:40:34,540 like explosions for examples. 688 00:40:34,540 --> 00:40:37,820 That's a pretty important for chemical engineers. 689 00:40:37,820 --> 00:40:41,810 But it's going to be really tough, because if it's 690 00:40:41,810 --> 00:40:44,990 got positive eigenvalues, then this is probably going 691 00:40:44,990 --> 00:40:46,850 to have a norm bigger than 1. 692 00:40:46,850 --> 00:40:48,320 And then whatever little errors we 693 00:40:48,320 --> 00:40:50,070 may get any step in the whole calculation, 694 00:40:50,070 --> 00:40:52,130 we're going to amplify and amplify and amplify. 695 00:40:52,130 --> 00:40:53,630 And who knows if we're going to have 696 00:40:53,630 --> 00:40:57,110 any reality left in our solution by the time we get to the end. 697 00:41:02,642 --> 00:41:04,350 All right, so that's one kind of problem. 698 00:41:09,520 --> 00:41:13,980 So how about we have a problem that's intrinsically stable. 699 00:41:13,980 --> 00:41:16,140 That means that all the eigenvalues 700 00:41:16,140 --> 00:41:19,320 of the Jacobian of the original problem are negative. 701 00:41:19,320 --> 00:41:22,770 So the thing should be perfectly stable. 702 00:41:22,770 --> 00:41:25,254 From whatever initial condition they started from, 703 00:41:25,254 --> 00:41:26,920 the whole thing should sort of come down 704 00:41:26,920 --> 00:41:28,300 to like an equilibrium point. 705 00:41:28,300 --> 00:41:29,340 So that's a pretty common situation 706 00:41:29,340 --> 00:41:30,506 in chemical engineering too. 707 00:41:30,506 --> 00:41:34,410 You start from some state, and it relaxes down 708 00:41:34,410 --> 00:41:37,040 to like an equilibrium state or a steady state. 709 00:41:37,040 --> 00:41:39,490 Well-behaved systems act like that. 710 00:41:39,490 --> 00:41:44,180 So what I'm going to show you is that, 711 00:41:44,180 --> 00:41:48,470 despite the fact that the df dy is well-behaved, 712 00:41:48,470 --> 00:41:52,460 it's possible that this sum matrix might still 713 00:41:52,460 --> 00:41:56,960 be poorly behaved and still have a norm bigger that one. 714 00:41:56,960 --> 00:41:58,960 And so that means you started for a problem that 715 00:41:58,960 --> 00:42:01,920 was pretty well-behaved, and you broke it by your choice 716 00:42:01,920 --> 00:42:04,030 in numerical method. 717 00:42:04,030 --> 00:42:06,550 So this is a very bad one. 718 00:42:06,550 --> 00:42:09,370 This is the kind that can make you really embarrassed 719 00:42:09,370 --> 00:42:11,440 if you're working in a company. 720 00:42:11,440 --> 00:42:13,000 Your boss gives you a problem. 721 00:42:13,000 --> 00:42:16,060 He says, please tell me the time is going to take our reactor 722 00:42:16,060 --> 00:42:18,070 to relax after a start up. 723 00:42:18,070 --> 00:42:19,774 the system is perfectly well-behaved. 724 00:42:19,774 --> 00:42:21,190 He gives these beautiful equations 725 00:42:21,190 --> 00:42:23,121 that your previous coworker worked out. 726 00:42:23,121 --> 00:42:23,870 Everything's fine. 727 00:42:23,870 --> 00:42:25,000 You put it in there. 728 00:42:25,000 --> 00:42:26,800 You use your stupid numerical solver, 729 00:42:26,800 --> 00:42:30,700 and you get a wild oscillation, and the thing blows up. 730 00:42:30,700 --> 00:42:33,610 OK, and that's going to be because of a poor choice of g 731 00:42:33,610 --> 00:42:39,700 so that it makes this norm of this matrix bigger than one. 732 00:42:39,700 --> 00:42:41,360 It's pretty easy for that to happen. 733 00:42:41,360 --> 00:42:44,792 So you don't need a very complicated case to show it. 734 00:42:52,930 --> 00:42:54,535 So let's just do a scalar case. 735 00:42:58,150 --> 00:43:00,210 Say dy dt is equal to negative 2y. 736 00:43:02,750 --> 00:43:05,370 So this is a well-behaved differential equation. 737 00:43:05,370 --> 00:43:08,120 What's the solution? 738 00:43:08,120 --> 00:43:14,335 So with a y naught y at t equals 0 is equal to 1. 739 00:43:17,800 --> 00:43:19,780 What's the solution? 740 00:43:19,780 --> 00:43:21,770 AUDIENCE: [INAUDIBLE] 741 00:43:21,770 --> 00:43:25,790 WILLIAM GREEN JR: Yeah, OK, so it's 742 00:43:25,790 --> 00:43:27,460 perfectly well-behaved function. 743 00:43:27,460 --> 00:43:31,460 And the plot when you plot it, here's one. 744 00:43:31,460 --> 00:43:33,470 It goes [PLANE DIVING NOISE]. 745 00:43:36,800 --> 00:43:39,890 So now let's solve this with the Forward Euler method. 746 00:43:39,890 --> 00:43:46,690 So I'm going to say that y new is equal to y 747 00:43:46,690 --> 00:43:48,946 plus delta t times f. 748 00:43:52,570 --> 00:43:58,360 So in this case, this is f, right there, negative 2y. 749 00:43:58,360 --> 00:44:00,200 So for Forward Euler, this is just 750 00:44:00,200 --> 00:44:01,926 saying that this is negative 2y. 751 00:44:05,350 --> 00:44:08,700 All right, so let's choose a delta t, 752 00:44:08,700 --> 00:44:12,350 say delta t equals 1 to start with. 753 00:44:12,350 --> 00:44:16,480 So we're going to compute. 754 00:44:16,480 --> 00:44:17,540 We start at t equals 0. 755 00:44:17,540 --> 00:44:21,540 We're computing out to t equals 1. so we're out here. 756 00:44:21,540 --> 00:44:24,400 So we put 1 in here. 757 00:44:24,400 --> 00:44:28,590 Our initial condition was y was 1, so that's 1. 758 00:44:28,590 --> 00:44:33,840 So as 1 times negative 2 times 1 is negative 2. 759 00:44:33,840 --> 00:44:35,410 This is a 1. 760 00:44:35,410 --> 00:44:40,770 1 minus 2 is negative 1. 761 00:44:40,770 --> 00:44:44,898 So at the first step, it goes from here down to there. 762 00:44:49,379 --> 00:44:50,420 This is not looking good. 763 00:44:54,200 --> 00:44:55,750 This physical variable is supposed 764 00:44:55,750 --> 00:44:57,770 to be gradually relaxing towards 0. 765 00:44:57,770 --> 00:45:00,580 Going negative is not what you want. 766 00:45:00,580 --> 00:45:02,990 And then I think we could go the next step, if you want. 767 00:45:02,990 --> 00:45:07,432 So we can put in suppose y is negative 1. 768 00:45:07,432 --> 00:45:09,640 And we're still doing a delta t of one step. 769 00:45:09,640 --> 00:45:14,610 So negative 1 times negative 2 is positive 2. 770 00:45:14,610 --> 00:45:18,731 Plus negative 1, so this goes back to the starting point. 771 00:45:18,731 --> 00:45:24,484 We'll get the sawtooth, is what you get. 772 00:45:24,484 --> 00:45:26,900 And it turns out, actually, if you make delta t any larger 773 00:45:26,900 --> 00:45:30,630 than 1, then this sawtooth amplifies, 774 00:45:30,630 --> 00:45:34,320 so you have a oscillating solution that explodes. 775 00:45:34,320 --> 00:45:37,170 Which is a little bit different than the physical solution, 776 00:45:37,170 --> 00:45:38,100 which gently relaxes. 777 00:45:42,230 --> 00:45:46,220 And so this is just a scalar equation 778 00:45:46,220 --> 00:45:48,260 that you guys could have done in high school. 779 00:45:48,260 --> 00:45:49,710 If you high school math teacher had been brave, 780 00:45:49,710 --> 00:45:51,160 he would have showed it you. 781 00:45:51,160 --> 00:45:54,682 But he don't want to ruin your faith in the Forward Euler 782 00:45:54,682 --> 00:45:57,140 method, so he just showed you a simple one that worked out. 783 00:46:00,320 --> 00:46:01,990 Now you can make this work better 784 00:46:01,990 --> 00:46:04,650 by making the delta t smaller. 785 00:46:04,650 --> 00:46:07,950 If you made the delta t small enough-- 786 00:46:07,950 --> 00:46:11,550 in this case, I think it's delta t is less than half-- 787 00:46:11,550 --> 00:46:15,250 then this actually behaves somewhat reasonably. 788 00:46:15,250 --> 00:46:21,420 And so what's going on there is that the identity matrix 789 00:46:21,420 --> 00:46:23,660 is positive. 790 00:46:23,660 --> 00:46:27,880 My Jacobian is negative. 791 00:46:27,880 --> 00:46:30,260 And if I make delta t small enough, 792 00:46:30,260 --> 00:46:32,755 this big positive number is bigger than this smaller 793 00:46:32,755 --> 00:46:34,880 negative number, and the whole thing stays positive 794 00:46:34,880 --> 00:46:38,910 but less than 1, and so then the normal is good. 795 00:46:38,910 --> 00:46:41,110 But if I make delta t too large, I 796 00:46:41,110 --> 00:46:43,615 make this second term much larger. 797 00:46:43,615 --> 00:46:46,390 If it gets much larger than 1, it overwhelms the first term. 798 00:46:46,390 --> 00:46:51,480 The thing's negative, but I what do the norm it's positive. 799 00:46:51,480 --> 00:46:53,360 And so it has an expansion. 800 00:46:53,360 --> 00:46:55,790 Then the negative is why we get these oscillations. 801 00:46:55,790 --> 00:46:58,290 So if you ever do a numerical differential equation solution 802 00:46:58,290 --> 00:47:00,330 and you start seeing this kind of stuff, 803 00:47:00,330 --> 00:47:06,090 it's almost certainly related to this numerical instability 804 00:47:06,090 --> 00:47:07,590 of the method. 805 00:47:07,590 --> 00:47:09,420 and not really a physical thing. 806 00:47:09,420 --> 00:47:10,140 I mean, there are physical things 807 00:47:10,140 --> 00:47:12,473 that do this, but not too often in chemical engineering, 808 00:47:12,473 --> 00:47:14,355 fortunately. 809 00:47:14,355 --> 00:47:18,230 All right, questions? 810 00:47:18,230 --> 00:47:19,040 No questions, OK. 811 00:47:19,040 --> 00:47:20,956 So in the textbook, they have a big discussion 812 00:47:20,956 --> 00:47:25,190 about how to evaluate different g methods in order 813 00:47:25,190 --> 00:47:27,800 to figure out if they're guaranteed to be stable, 814 00:47:27,800 --> 00:47:30,410 or under what conditions they'll be stable. 815 00:47:30,410 --> 00:47:34,184 And so that's worth a look to see. 816 00:47:34,184 --> 00:47:35,600 And they have detailed derivations 817 00:47:35,600 --> 00:47:36,543 for several of them. 818 00:47:40,000 --> 00:47:42,444 Questions? 819 00:47:42,444 --> 00:47:44,360 All right, so for the homework, homework four, 820 00:47:44,360 --> 00:47:47,870 you're going to have to write your own ODE solver. 821 00:47:47,870 --> 00:47:52,390 So you'll do update formulas and compare them. 822 00:47:52,390 --> 00:47:56,180 When we come back on Tuesday, I'll 823 00:47:56,180 --> 00:47:57,830 show you how this is usually dealt 824 00:47:57,830 --> 00:48:04,530 with by switching to implicit solvers and talk about that. 825 00:48:04,530 --> 00:48:06,980 All right, have a good weekend. 826 00:48:06,980 --> 00:48:09,470 Enjoy the three days off.