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:24,488 --> 00:00:27,604 PROFESSOR: You get your quizzes back on Monday. 9 00:00:27,604 --> 00:00:29,520 Maybe you've had some time to think about them 10 00:00:29,520 --> 00:00:31,710 and reflect on your answers and the problems 11 00:00:31,710 --> 00:00:34,950 that were posed, maybe how they compared to the quiz that-- 12 00:00:34,950 --> 00:00:36,283 sample quiz that was put online. 13 00:00:36,283 --> 00:00:38,850 I think it was pretty fair exam overall. 14 00:00:38,850 --> 00:00:40,740 It looked very similar, actually, to the exam 15 00:00:40,740 --> 00:00:43,280 that we've given in the past, but we tried to-- 16 00:00:43,280 --> 00:00:46,860 tried not to put any sort of chemical engineering problems 17 00:00:46,860 --> 00:00:48,720 on top of the numerics, which tends 18 00:00:48,720 --> 00:00:51,160 to lead to confusion basically. 19 00:00:51,160 --> 00:00:52,830 Spent a lot of time reading and trying 20 00:00:52,830 --> 00:00:55,470 to understand what the underlying physical problem is 21 00:00:55,470 --> 00:01:00,450 instead of just testing the numerical methods. 22 00:01:00,450 --> 00:01:02,050 Are there any questions about the quiz 23 00:01:02,050 --> 00:01:02,940 now that you've gotten it back? 24 00:01:02,940 --> 00:01:04,920 You've seen your answers, things that were unclear that 25 00:01:04,920 --> 00:01:05,753 were asked on there? 26 00:01:05,753 --> 00:01:10,780 Now is a good time for reflection if there are. 27 00:01:10,780 --> 00:01:12,200 No. 28 00:01:12,200 --> 00:01:16,090 I heard this last homework assignment was really long, so 29 00:01:16,090 --> 00:01:16,800 sorry about that. 30 00:01:16,800 --> 00:01:17,800 I didn't write this one. 31 00:01:17,800 --> 00:01:19,285 [LAUGHS] 32 00:01:19,285 --> 00:01:21,730 I'm not going to Take credit for it. 33 00:01:21,730 --> 00:01:25,885 It's hard to judge sometimes how long these things take. 34 00:01:25,885 --> 00:01:30,280 A lot of it has to do with how facile you are with coding up 35 00:01:30,280 --> 00:01:31,192 some of the problems. 36 00:01:31,192 --> 00:01:32,650 So conceptually they can be simple, 37 00:01:32,650 --> 00:01:34,233 but the coding can be quite difficult, 38 00:01:34,233 --> 00:01:35,590 so it's hard to estimate. 39 00:01:35,590 --> 00:01:38,080 I think this week's, which should be posted 40 00:01:38,080 --> 00:01:41,020 and is on the topic of DAE, should be easier. 41 00:01:41,020 --> 00:01:43,380 It's only got two problems instead of three. 42 00:01:43,380 --> 00:01:45,130 You're likely going to be able to leverage 43 00:01:45,130 --> 00:01:46,780 some of the code from your solution 44 00:01:46,780 --> 00:01:49,160 last week for this week. 45 00:01:49,160 --> 00:01:52,210 So I don't think it should be quite so challenging. 46 00:01:52,210 --> 00:01:56,530 There is a last part to the first problem 47 00:01:56,530 --> 00:01:59,020 that seems to always give first-year graduate 48 00:01:59,020 --> 00:02:00,130 students trouble. 49 00:02:00,130 --> 00:02:02,470 Somehow you guys don't know how to write down 50 00:02:02,470 --> 00:02:05,050 energy balances on reactors. 51 00:02:05,050 --> 00:02:07,270 So you might ask pointed questions 52 00:02:07,270 --> 00:02:10,210 about this at your office hours, so you don't spend time trying 53 00:02:10,210 --> 00:02:14,200 to solve the wrong equations. 54 00:02:14,200 --> 00:02:14,780 This happens. 55 00:02:14,780 --> 00:02:19,060 I think many of you maybe were taught incorrectly or just 56 00:02:19,060 --> 00:02:22,690 forgot at some point how to write these things down. 57 00:02:22,690 --> 00:02:25,711 So the office hour is a good point to bring that up, 58 00:02:25,711 --> 00:02:27,460 and that way you don't spend a lot of time 59 00:02:27,460 --> 00:02:32,170 troubling yourself with the wrong system of equations. 60 00:02:32,170 --> 00:02:34,051 Any questions? 61 00:02:34,051 --> 00:02:34,550 Yes. 62 00:02:34,550 --> 00:02:39,542 AUDIENCE: [INAUDIBLE] 63 00:02:39,542 --> 00:02:40,500 PROFESSOR: What's that? 64 00:02:40,500 --> 00:02:44,980 I'd-- well, to be fair, almost all the lectures for the rest 65 00:02:44,980 --> 00:02:47,350 of the class are going to be taught by Professor Green. 66 00:02:47,350 --> 00:02:49,940 I'll step in and do some review for two sessions, 67 00:02:49,940 --> 00:02:53,710 review before the second quiz and review before the final. 68 00:02:53,710 --> 00:02:55,480 Since he's leading the lectures, he'll 69 00:02:55,480 --> 00:02:59,200 be responsible for largely generating the homework 70 00:02:59,200 --> 00:02:59,920 assignments. 71 00:02:59,920 --> 00:03:03,760 I'll try to help out with that as best I can, 72 00:03:03,760 --> 00:03:07,420 but this will be the last formal lecture 73 00:03:07,420 --> 00:03:10,210 you get from me this term. 74 00:03:10,210 --> 00:03:11,520 Other questions? 75 00:03:11,520 --> 00:03:13,720 Well, that's very, very sweet of you. 76 00:03:13,720 --> 00:03:17,040 Thank you so much. 77 00:03:17,040 --> 00:03:21,670 Least there wasn't any hushed yes! 78 00:03:21,670 --> 00:03:23,765 Today's lecture is our last one on DAE, 79 00:03:23,765 --> 00:03:26,512 so we're only going to do two. 80 00:03:26,512 --> 00:03:28,720 You're going to see today that differential algebraic 81 00:03:28,720 --> 00:03:31,960 equations are pretty complicated actually. 82 00:03:31,960 --> 00:03:34,870 And when they get sufficiently complicated, 83 00:03:34,870 --> 00:03:39,850 you really need to reach out to existing codes that 84 00:03:39,850 --> 00:03:42,790 are designed to solve particular classes of differential 85 00:03:42,790 --> 00:03:43,840 algebraic equations. 86 00:03:43,840 --> 00:03:46,270 So for you, it's more important to be 87 00:03:46,270 --> 00:03:49,600 able to identify in the models you formulate when 88 00:03:49,600 --> 00:03:52,420 these complications arise and what 89 00:03:52,420 --> 00:03:54,880 are the essential ingredients of the model that 90 00:03:54,880 --> 00:03:56,680 can be put into one of these solvers 91 00:03:56,680 --> 00:03:59,124 that you get a result that isn't nonsense. 92 00:03:59,124 --> 00:04:00,790 And that's what we're going to do today. 93 00:04:00,790 --> 00:04:02,510 There's going to be lots of examples to work through. 94 00:04:02,510 --> 00:04:04,180 So make sure you're sitting next to somebody 95 00:04:04,180 --> 00:04:05,800 you like because I'm going to ask you to try 96 00:04:05,800 --> 00:04:07,570 to think about these things and discuss things 97 00:04:07,570 --> 00:04:08,380 as we go through. 98 00:04:10,857 --> 00:04:12,940 Let me review where these complications come from. 99 00:04:12,940 --> 00:04:15,070 So last time we talked briefly about semi-explicit 100 00:04:15,070 --> 00:04:18,256 and fully-implicit differential algebraic equations. 101 00:04:18,256 --> 00:04:20,589 I told you in principle, you could simulate these things 102 00:04:20,589 --> 00:04:22,620 with backwards difference formulas 103 00:04:22,620 --> 00:04:25,327 in solving nonlinear equations at each time step 104 00:04:25,327 --> 00:04:26,410 and just marching forward. 105 00:04:26,410 --> 00:04:30,320 Did this in your homework for ODE IVP 106 00:04:30,320 --> 00:04:32,950 So you can do the same thing for differential algebraic 107 00:04:32,950 --> 00:04:33,670 equations. 108 00:04:33,670 --> 00:04:37,210 But there was a catch to that, and I illustrated that 109 00:04:37,210 --> 00:04:38,020 at the very end. 110 00:04:38,020 --> 00:04:40,980 So maybe your brain was fried at the end, and you missed this. 111 00:04:40,980 --> 00:04:42,740 It's good to recap. 112 00:04:42,740 --> 00:04:45,070 So I showed you an example of a stirred tank 113 00:04:45,070 --> 00:04:47,260 where you had some-- 114 00:04:47,260 --> 00:04:50,410 you had transport of some solute into the stirred tank, 115 00:04:50,410 --> 00:04:52,970 and then you're pulling it out at a different concentration. 116 00:04:52,970 --> 00:04:55,970 And we're trying to track the dynamics of the system, 117 00:04:55,970 --> 00:04:58,114 the concentration in and the concentration out, 118 00:04:58,114 --> 00:04:59,530 and we imagined a problem in which 119 00:04:59,530 --> 00:05:01,630 we measured and the concentration in 120 00:05:01,630 --> 00:05:04,390 and we tried to predict the concentration out. 121 00:05:04,390 --> 00:05:08,590 So we had a system of differential and algebraic 122 00:05:08,590 --> 00:05:10,540 equations to solve. 123 00:05:10,540 --> 00:05:13,600 And I showed you that if I used the backward Euler method, 124 00:05:13,600 --> 00:05:16,710 the lowest order backwards difference formula that I can-- 125 00:05:16,710 --> 00:05:18,730 of the canonical class of these things 126 00:05:18,730 --> 00:05:22,330 that I can craft, where I approximate the derivative 127 00:05:22,330 --> 00:05:27,250 of an unknown with a relative error or an error proportional 128 00:05:27,250 --> 00:05:31,502 to the time step delta t, that carrying out one time 129 00:05:31,502 --> 00:05:33,460 step with this backward differentiation formula 130 00:05:33,460 --> 00:05:35,990 would determine c1 in principle exactly-- 131 00:05:35,990 --> 00:05:38,765 if I knew gamma exactly, I would know c1 exactly. 132 00:05:38,765 --> 00:05:42,610 And it would determine c2 to within order delta t-squared. 133 00:05:42,610 --> 00:05:46,160 It's the local truncation error here was order delta t-squared. 134 00:05:46,160 --> 00:05:48,190 So you just substitute this formula in, 135 00:05:48,190 --> 00:05:50,860 and you ask does this error term change at all. 136 00:05:50,860 --> 00:05:52,900 At some point, I wind up multiplying by delta t, 137 00:05:52,900 --> 00:05:56,240 and so I go from order delta t to order delta t-squared It's 138 00:05:56,240 --> 00:05:58,240 the numerical error that gets carried around 139 00:05:58,240 --> 00:06:01,790 in this calculation. 140 00:06:01,790 --> 00:06:04,180 I just switch the model a little bit. 141 00:06:04,180 --> 00:06:07,140 It seems like a irrelevant change, 142 00:06:07,140 --> 00:06:09,210 but it turns out to be incredibly significant. 143 00:06:09,210 --> 00:06:11,520 So now I imagine a different problem 144 00:06:11,520 --> 00:06:15,490 in which I'm measuring c2 the outlet, 145 00:06:15,490 --> 00:06:19,740 and I want to predict c1 the inlet. 146 00:06:19,740 --> 00:06:22,290 I still have a system of DAEs that I have to solve, 147 00:06:22,290 --> 00:06:24,690 and if I applied the backwards Euler method, 148 00:06:24,690 --> 00:06:26,820 well, c2 is automatically determined 149 00:06:26,820 --> 00:06:28,350 by this algebraic equation. 150 00:06:28,350 --> 00:06:30,060 So I know that exactly. 151 00:06:30,060 --> 00:06:33,420 I got to go up to this first equation and solve it for c1, 152 00:06:33,420 --> 00:06:36,630 and when I do, I have to have an approximation 153 00:06:36,630 --> 00:06:37,500 for the derivative. 154 00:06:37,500 --> 00:06:40,576 And the derivative is proportional to delta t-- 155 00:06:40,576 --> 00:06:43,200 carries an error around with it that's proportional to delta t. 156 00:06:43,200 --> 00:06:47,879 So I know c1 to within order delta t, not order delta 157 00:06:47,879 --> 00:06:49,420 t-squared like in the previous model. 158 00:06:49,420 --> 00:06:51,300 So there's something fundamentally different 159 00:06:51,300 --> 00:06:53,130 about these two circumstances. 160 00:06:53,130 --> 00:06:57,550 And all I did was go from c1 to c2 in this algebraic equation. 161 00:06:57,550 --> 00:06:59,506 So that's peculiar, and that should 162 00:06:59,506 --> 00:07:01,380 make you really suspicious about your ability 163 00:07:01,380 --> 00:07:04,260 to solve these problems accurately. 164 00:07:04,260 --> 00:07:06,010 And here was the third example I gave you. 165 00:07:08,740 --> 00:07:12,010 So I said imagine this system of DAEs instead. 166 00:07:12,010 --> 00:07:13,750 So here's the differential equation. 167 00:07:13,750 --> 00:07:15,250 Here's an algebraic equation. 168 00:07:15,250 --> 00:07:17,340 Apply the backward Euler method. 169 00:07:17,340 --> 00:07:19,630 Well, c3 is determined automatically 170 00:07:19,630 --> 00:07:21,220 by this algebraic equation. 171 00:07:21,220 --> 00:07:22,680 I know it exactly. 172 00:07:22,680 --> 00:07:25,370 C2 is related to the derivative of c3, 173 00:07:25,370 --> 00:07:28,270 so I need to approximate it with my backward Euler derivative. 174 00:07:28,270 --> 00:07:31,222 And that picks up an error, order delta t. 175 00:07:31,222 --> 00:07:35,110 C1 is equal to the derivative of c2, 176 00:07:35,110 --> 00:07:37,270 so I need an approximation for c2. 177 00:07:37,270 --> 00:07:40,480 The derivative of c2 that's the backwards Euler approximation. 178 00:07:40,480 --> 00:07:44,070 That has an error that's proportional to order delta t. 179 00:07:44,070 --> 00:07:49,270 But c2 itself also has an error proportional to order delta t. 180 00:07:49,270 --> 00:07:53,110 And order delta t divided by delta t is order 1. 181 00:07:53,110 --> 00:07:54,250 So I get c3. 182 00:07:54,250 --> 00:07:55,660 I get c2. 183 00:07:55,660 --> 00:07:57,760 I even know what c1 is. 184 00:07:57,760 --> 00:08:00,280 I solve this problem with my backwards difference formula, 185 00:08:00,280 --> 00:08:06,230 and I have no resolution of c1, no clue what that value is. 186 00:08:06,230 --> 00:08:08,410 So this is some complicated control scheme 187 00:08:08,410 --> 00:08:11,260 that I've set up, and I need to know the value of c1 188 00:08:11,260 --> 00:08:13,210 to figure out how to operate this process. 189 00:08:13,210 --> 00:08:14,590 I'm lost. 190 00:08:14,590 --> 00:08:16,780 This is never going to work. 191 00:08:16,780 --> 00:08:19,330 So these three problems have fundamental differences 192 00:08:19,330 --> 00:08:20,560 between them. 193 00:08:20,560 --> 00:08:23,170 And I'm going to show you how to predict 194 00:08:23,170 --> 00:08:27,746 when these differences are going to occur, how to name them. 195 00:08:27,746 --> 00:08:30,180 So the stirred tank example 1, it 196 00:08:30,180 --> 00:08:33,110 carried a local truncation error, one time step order 197 00:08:33,110 --> 00:08:34,387 delta-t squared. 198 00:08:34,387 --> 00:08:36,720 Stirred tank example 2 carried a local truncation error. 199 00:08:36,720 --> 00:08:38,610 It's order delta-t. 200 00:08:38,610 --> 00:08:41,010 This third DAE example had a local truncation error 201 00:08:41,010 --> 00:08:41,730 that's order 1. 202 00:08:41,730 --> 00:08:45,030 No change in delta-t will improve my solution 203 00:08:45,030 --> 00:08:46,710 to that problem. 204 00:08:46,710 --> 00:08:50,570 It's independent of how I choose to do my time setting, which, 205 00:08:50,570 --> 00:08:54,757 of course, is ridiculous. 206 00:08:54,757 --> 00:08:55,840 So here's another problem. 207 00:08:59,490 --> 00:09:02,942 And I'd like you to try to do this. 208 00:09:02,942 --> 00:09:04,650 We'll see how far you can get through it. 209 00:09:04,650 --> 00:09:05,820 You don't have to get all the way through it, 210 00:09:05,820 --> 00:09:07,861 but see how far you can get through this problem. 211 00:09:07,861 --> 00:09:10,280 So can you apply the backward Euler method 212 00:09:10,280 --> 00:09:14,634 to this system of DAEs and try to predict how 213 00:09:14,634 --> 00:09:15,425 the air propagates. 214 00:09:15,425 --> 00:09:18,834 [INTERPOSING VOICES] 215 00:12:14,760 --> 00:12:17,410 OK, I don't want to break up the conversation too much, 216 00:12:17,410 --> 00:12:21,020 but tell me what you're finding as you try to do this. 217 00:12:21,020 --> 00:12:21,520 Yes. 218 00:12:21,520 --> 00:12:23,795 AUDIENCE: We found that c3 has to be 0. 219 00:12:23,795 --> 00:12:25,510 PROFESSOR: OK. 220 00:12:25,510 --> 00:12:27,610 How did you find that c3 has to be 0? 221 00:12:27,610 --> 00:12:29,466 AUDIENCE: Stick in this [INAUDIBLE] 222 00:12:29,466 --> 00:12:31,466 that's what the difference is of the [INAUDIBLE] 223 00:12:31,466 --> 00:12:32,250 PROFESSOR: Yes. 224 00:12:32,250 --> 00:12:34,400 AUDIENCE: [INAUDIBLE] 225 00:12:34,400 --> 00:12:35,130 PROFESSOR: Right. 226 00:12:35,130 --> 00:12:35,910 AUDIENCE: [INAUDIBLE] 227 00:12:35,910 --> 00:12:36,618 PROFESSOR: Right. 228 00:12:36,618 --> 00:12:38,536 AUDIENCE: [INAUDIBLE] 229 00:12:38,536 --> 00:12:40,655 PROFESSOR: OK, so this is very clever, 230 00:12:40,655 --> 00:12:45,500 and this has to do with how the model itself is formulated. 231 00:12:45,500 --> 00:12:49,040 So you did manipulations with the fundamental equations 232 00:12:49,040 --> 00:12:54,470 that I wrote down to determine that c3 had to be 0. 233 00:12:54,470 --> 00:12:55,465 That's good. 234 00:12:55,465 --> 00:12:56,090 We can do that. 235 00:12:56,090 --> 00:12:56,765 We can look at the equations. 236 00:12:56,765 --> 00:12:59,060 We can figure how to eliminate variables and find out-- 237 00:12:59,060 --> 00:13:01,604 c3 in this case always has to be 0. 238 00:13:01,604 --> 00:13:03,770 That may or may not be obvious, but one way to do it 239 00:13:03,770 --> 00:13:07,264 is take this constraint equation and take its derivative. 240 00:13:07,264 --> 00:13:08,930 This has to hold at every point in time, 241 00:13:08,930 --> 00:13:11,520 so its derivative must hold at every point in time. 242 00:13:11,520 --> 00:13:13,220 Substitute these other two equations 243 00:13:13,220 --> 00:13:18,080 in, c1 dot plus c2 dot eliminates c1 and c2, 244 00:13:18,080 --> 00:13:20,090 and all that will be left is c3 has to be 0. 245 00:13:20,090 --> 00:13:21,950 So if I manipulate these equations, 246 00:13:21,950 --> 00:13:24,369 I know c3 has to be 0. 247 00:13:24,369 --> 00:13:26,660 Let's suppose I don't manipulate the equations, though. 248 00:13:26,660 --> 00:13:30,080 Let's suppose I just put in the backward Euler approximation 249 00:13:30,080 --> 00:13:31,640 over here for the derivative, and I 250 00:13:31,640 --> 00:13:35,930 evaluate c1, c2 and c3 at the current time 251 00:13:35,930 --> 00:13:39,890 and ask what is c1, c2, and c3 at the current time in terms 252 00:13:39,890 --> 00:13:41,369 of what it was at the previous time 253 00:13:41,369 --> 00:13:43,160 And you get an answer that looks like this. 254 00:13:43,160 --> 00:13:44,930 So there's a simple system of equations 255 00:13:44,930 --> 00:13:45,960 that one has to solve. 256 00:13:45,960 --> 00:13:48,500 It's a three-by-three system, and it's easy to eliminate. 257 00:13:48,500 --> 00:13:49,910 I didn't really expect you-- 258 00:13:49,910 --> 00:13:51,230 actually expected this answer to come up 259 00:13:51,230 --> 00:13:52,771 because you guys are all very clever, 260 00:13:52,771 --> 00:13:54,770 so you try to make things easy beforehand. 261 00:13:54,770 --> 00:13:56,900 But it's-- there's a problem here. 262 00:13:56,900 --> 00:13:59,030 If I try to just do the backwards difference 263 00:13:59,030 --> 00:14:03,670 formula with the equations as they're written, 264 00:14:03,670 --> 00:14:05,420 I'll find out that when I solve this thing 265 00:14:05,420 --> 00:14:10,390 that c1 will be c1 minus c2 at the previous time step over 2 266 00:14:10,390 --> 00:14:12,990 plus an order delta t-squared error. 267 00:14:12,990 --> 00:14:15,950 C2 will take on the opposite sign plus an order delta 268 00:14:15,950 --> 00:14:18,590 t-squared error. 269 00:14:18,590 --> 00:14:25,280 C3 will be minus c1 plus c2 over 2 delta t 270 00:14:25,280 --> 00:14:28,580 plus an order delta t error. 271 00:14:28,580 --> 00:14:32,900 And actually if I were to apply successive approximations, 272 00:14:32,900 --> 00:14:34,820 different-- if I step-- 273 00:14:34,820 --> 00:14:37,310 take many steps in a row with this backward Euler method, 274 00:14:37,310 --> 00:14:39,601 there's nothing to guarantee that these constraints are 275 00:14:39,601 --> 00:14:42,260 satisfied exactly. 276 00:14:42,260 --> 00:14:44,330 I can have numerical error in my solution 277 00:14:44,330 --> 00:14:46,700 of this algebraic equation. 278 00:14:46,700 --> 00:14:49,340 So it isn't necessarily true that in my numerical solution 279 00:14:49,340 --> 00:14:52,400 c1 plus c2 is equal to 0. 280 00:14:52,400 --> 00:14:54,312 I can show you right here if I add c1 to c2. 281 00:14:54,312 --> 00:14:55,520 Well, these two are the same. 282 00:14:55,520 --> 00:14:58,370 They'll cancel, but am I guaranteed that these errors 283 00:14:58,370 --> 00:14:59,120 will cancel, too? 284 00:14:59,120 --> 00:15:00,786 Or will there be a small numerical error 285 00:15:00,786 --> 00:15:03,490 that propagates? 286 00:15:03,490 --> 00:15:06,170 So model formulation is key. 287 00:15:06,170 --> 00:15:09,530 This problem is like stirred tank problem 2. 288 00:15:09,530 --> 00:15:12,650 One of our solutions lost in order of accuracy 289 00:15:12,650 --> 00:15:14,279 in the local truncation error. 290 00:15:14,279 --> 00:15:15,320 It's not delta t-squared. 291 00:15:15,320 --> 00:15:18,160 It's order delta t. 292 00:15:18,160 --> 00:15:21,170 And what you want it to do, which is not necessarily 293 00:15:21,170 --> 00:15:23,840 a bad thing, was to try to use this constraint 294 00:15:23,840 --> 00:15:27,801 to somehow eliminate equations and simplify things. 295 00:15:27,801 --> 00:15:29,300 That can be helpful, but it can also 296 00:15:29,300 --> 00:15:31,460 be true that if I use this constraint to simplify 297 00:15:31,460 --> 00:15:34,580 equations and eliminate things, my numerical solutions 298 00:15:34,580 --> 00:15:36,860 may no longer satisfy the constraint at all. 299 00:15:36,860 --> 00:15:41,554 They may drift away from that algebraic equation. 300 00:15:41,554 --> 00:15:42,470 Does that makes sense? 301 00:15:42,470 --> 00:15:44,928 I'm not controlling the error with respect to that equation 302 00:15:44,928 --> 00:15:45,587 anymore. 303 00:15:50,950 --> 00:15:55,810 So how do you give a name to this kind of behavior? 304 00:15:55,810 --> 00:16:00,340 And the thing one talks about is the differential index 305 00:16:00,340 --> 00:16:02,017 of a DAE system. 306 00:16:02,017 --> 00:16:03,850 So let's look at the stirred tank example 1, 307 00:16:03,850 --> 00:16:06,550 and let's ask this question. 308 00:16:06,550 --> 00:16:09,400 How many time derivatives are needed 309 00:16:09,400 --> 00:16:12,910 to convert to a system of independent ODEs having 310 00:16:12,910 --> 00:16:14,770 differentials of all the unknowns. 311 00:16:14,770 --> 00:16:16,960 So there are two unknowns in this system. 312 00:16:16,960 --> 00:16:19,570 One is c1, and the other c2. 313 00:16:19,570 --> 00:16:20,980 Actually have one equation, which 314 00:16:20,980 --> 00:16:23,820 contains a differential of c2. 315 00:16:23,820 --> 00:16:25,960 How many derivatives of these equations with time 316 00:16:25,960 --> 00:16:31,660 do I need to take in order to get an equivalent ODE system? 317 00:16:31,660 --> 00:16:32,700 And it's just one. 318 00:16:32,700 --> 00:16:35,250 If I take the derivative of equation 2, 319 00:16:35,250 --> 00:16:37,140 I'll get a differential equation for c1. 320 00:16:37,140 --> 00:16:39,510 Now, I have a differential equation for c1 321 00:16:39,510 --> 00:16:42,540 and a differential equation for c2. 322 00:16:42,540 --> 00:16:44,280 DAEs of this type, where I only have 323 00:16:44,280 --> 00:16:47,790 to take the time derivative once of the algebraic constraints, 324 00:16:47,790 --> 00:16:51,664 are called index 1 DAEs. 325 00:16:51,664 --> 00:16:53,830 We saw that when we applied the backwards difference 326 00:16:53,830 --> 00:16:56,097 formula to stirred tank example 1, 327 00:16:56,097 --> 00:16:58,180 the local truncation error was the same as what we 328 00:16:58,180 --> 00:17:01,070 would get in an ODE IVP system. 329 00:17:01,070 --> 00:17:02,320 So index 1-- 330 00:17:02,320 --> 00:17:04,720 DAEs are easy to solve. 331 00:17:04,720 --> 00:17:08,619 They behave like ODE IVPs. 332 00:17:08,619 --> 00:17:10,626 You determine whether it's index 1 or not 333 00:17:10,626 --> 00:17:12,000 by asking how many derivatives do 334 00:17:12,000 --> 00:17:16,220 I have to take in order to get a system of independent ODEs. 335 00:17:16,220 --> 00:17:18,770 So I put this together with this, 336 00:17:18,770 --> 00:17:20,300 I have two independent ODEs. 337 00:17:20,300 --> 00:17:23,594 I could solve these, subject to some set of initial conditions, 338 00:17:23,594 --> 00:17:25,010 and the solution might be the same 339 00:17:25,010 --> 00:17:27,760 as the solution to the DAE. 340 00:17:27,760 --> 00:17:29,690 Let's do stirred tank example 2 now. 341 00:17:33,780 --> 00:17:34,710 So here we go. 342 00:17:34,710 --> 00:17:36,390 They look very similar, but now c2 343 00:17:36,390 --> 00:17:38,700 appears in the algebraic equation instead of c1. 344 00:17:38,700 --> 00:17:42,360 So let's take a derivative of the algebraic equation, 345 00:17:42,360 --> 00:17:44,960 and I get a differential equation for c2. 346 00:17:44,960 --> 00:17:47,250 But I already had a differential equation for c2. 347 00:17:47,250 --> 00:17:49,620 I want a differential equation for c1. 348 00:17:49,620 --> 00:17:50,660 So what do I do? 349 00:17:50,660 --> 00:17:53,380 While I know dc2 dt from up here. 350 00:17:53,380 --> 00:17:56,740 So let's substitute that in and solve for c1. 351 00:17:56,740 --> 00:17:58,260 So a substitute equation 1 in here. 352 00:17:58,260 --> 00:17:59,590 I solve for c1. 353 00:17:59,590 --> 00:18:01,320 Now, let's take another derivative. 354 00:18:01,320 --> 00:18:02,220 Call this equation 3. 355 00:18:02,220 --> 00:18:04,770 Derivative of equation 3 now gives me 356 00:18:04,770 --> 00:18:07,080 a differential equation for c1. 357 00:18:07,080 --> 00:18:09,770 It's also in terms of dc2 dt. 358 00:18:09,770 --> 00:18:10,560 If I want-- 359 00:18:10,560 --> 00:18:12,570 I don't have to-- but if I want to, 360 00:18:12,570 --> 00:18:14,370 I can substitute for that again to just get 361 00:18:14,370 --> 00:18:17,490 dc1 dt in terms of c1 and c2. 362 00:18:17,490 --> 00:18:21,390 But it took two derivatives to generate 363 00:18:21,390 --> 00:18:26,430 a system of independent ODEs from the DAEs. 364 00:18:26,430 --> 00:18:30,249 And this is called index 2. 365 00:18:30,249 --> 00:18:32,040 So it's a different character from index 1. 366 00:18:32,040 --> 00:18:33,070 We saw what happens. 367 00:18:33,070 --> 00:18:35,940 We tend to lose an order of accuracy at index 2. 368 00:18:35,940 --> 00:18:39,600 Somehow index 2 problems are more sensitive than index 1 369 00:18:39,600 --> 00:18:41,552 problems. 370 00:18:41,552 --> 00:18:42,510 Here's another example. 371 00:18:45,370 --> 00:18:49,610 So DAE sample 3 is going to proceed the same way. 372 00:18:49,610 --> 00:18:53,549 So first we need a differential equation. 373 00:18:53,549 --> 00:18:54,590 Which one's missing here? 374 00:18:54,590 --> 00:18:56,381 We don't have differential equation for c1, 375 00:18:56,381 --> 00:18:58,280 and that's what we're hunting for. 376 00:18:58,280 --> 00:19:00,590 So let's take a derivative of the third equation, 377 00:19:00,590 --> 00:19:01,640 the algebraic equation. 378 00:19:01,640 --> 00:19:03,710 We get c3 dot is gamma dot. 379 00:19:03,710 --> 00:19:05,630 And we already have an equation for c3 dot. 380 00:19:05,630 --> 00:19:07,000 That was 2. 381 00:19:07,000 --> 00:19:08,780 So substitute 2 in. 382 00:19:08,780 --> 00:19:12,039 So now we got c2 is gamma dot, but we're really 383 00:19:12,039 --> 00:19:13,580 looking for something in terms of c1, 384 00:19:13,580 --> 00:19:15,920 so let's take the derivative again. 385 00:19:15,920 --> 00:19:19,140 We've got c2 [? dot ?] is gamma dot dot. 386 00:19:19,140 --> 00:19:21,930 Doesn't help us much, but c2 dot we know is equal to c1. 387 00:19:21,930 --> 00:19:25,110 So we substitute one more time and take a third derivative. 388 00:19:25,110 --> 00:19:29,080 Now, we have a differential equation for c1. 389 00:19:29,080 --> 00:19:31,170 And we have differential equations for c2 and c3. 390 00:19:33,910 --> 00:19:36,340 We can take this is the differential equation for c3. 391 00:19:36,340 --> 00:19:38,060 This is the differential equation for c2. 392 00:19:38,060 --> 00:19:41,260 This is the differential equation for c1. 393 00:19:41,260 --> 00:19:44,500 Some subset of these can be chosen, 394 00:19:44,500 --> 00:19:46,960 and we can replace this algebraic equation 395 00:19:46,960 --> 00:19:49,280 with a differential equation. 396 00:19:49,280 --> 00:19:53,154 This is an index 3 DAE. 397 00:19:53,154 --> 00:19:55,070 It's not always a good idea to replace the DAE 398 00:19:55,070 --> 00:19:58,370 system with these derivatives. 399 00:19:58,370 --> 00:19:59,330 I can write down-- 400 00:19:59,330 --> 00:20:03,246 I can have a particular function that's equal to 0 401 00:20:03,246 --> 00:20:05,120 and higher-order derivatives of that function 402 00:20:05,120 --> 00:20:06,620 are also equal to 0. 403 00:20:06,620 --> 00:20:08,900 But the solution to those differential equations 404 00:20:08,900 --> 00:20:11,477 are not necessarily equal to that function. 405 00:20:11,477 --> 00:20:14,060 There are lots of functions that might have this same property 406 00:20:14,060 --> 00:20:16,435 that its derivative-- certain number of derivatives of it 407 00:20:16,435 --> 00:20:17,750 are equal to 0. 408 00:20:17,750 --> 00:20:19,610 So there have to be extra initial conditions 409 00:20:19,610 --> 00:20:21,620 or constraints on the solution that 410 00:20:21,620 --> 00:20:24,350 confine me to the manifold of solutions 411 00:20:24,350 --> 00:20:25,740 that reflects the DAEs. 412 00:20:25,740 --> 00:20:28,679 We'll talk about consistent initialization later on. 413 00:20:28,679 --> 00:20:30,220 And that's going to fix this problem. 414 00:20:30,220 --> 00:20:32,000 So it's not always a good idea to do this, 415 00:20:32,000 --> 00:20:33,440 but it's one way of understanding, 416 00:20:33,440 --> 00:20:38,570 given this model, what sort of sensitivity is it can exhibit. 417 00:20:38,570 --> 00:20:40,503 We try to calculate its differential index. 418 00:20:40,503 --> 00:20:41,003 Yes? 419 00:20:41,003 --> 00:20:47,282 AUDIENCE: So [INAUDIBLE] 420 00:20:47,282 --> 00:20:49,216 PROFESSOR: Yes. 421 00:20:49,216 --> 00:20:49,715 Yes. 422 00:20:49,715 --> 00:20:50,772 AUDIENCE: [INAUDIBLE] 423 00:20:50,772 --> 00:20:52,730 PROFESSOR: Well, we have to be careful. 424 00:20:52,730 --> 00:20:55,880 We need to choose a set that are independent of the others. 425 00:20:55,880 --> 00:20:56,750 OK. 426 00:20:56,750 --> 00:20:59,135 So like these three are clearly going to be independent, 427 00:20:59,135 --> 00:21:01,850 and it's going to be OK. 428 00:21:01,850 --> 00:21:04,687 It's going to be a problem if I choose this one and these two. 429 00:21:04,687 --> 00:21:06,770 They're not going to be independent of each other. 430 00:21:06,770 --> 00:21:08,895 So we have to select independent ones from the set. 431 00:21:08,895 --> 00:21:11,525 It may not be obvious what independent means in general. 432 00:21:11,525 --> 00:21:13,211 AUDIENCE: [INAUDIBLE] 433 00:21:13,211 --> 00:21:15,631 PROFESSOR: 1, 2, 6 could also work, yes. 434 00:21:15,631 --> 00:21:16,130 Yes. 435 00:21:21,500 --> 00:21:25,906 So in general-- not in general. 436 00:21:25,906 --> 00:21:27,530 Let's talk about the differential index 437 00:21:27,530 --> 00:21:29,810 of a semi-explicit DAE system. 438 00:21:29,810 --> 00:21:35,090 Remember semi-explicit meant that the differential variables 439 00:21:35,090 --> 00:21:37,940 can be written as x dot or dx dt is 440 00:21:37,940 --> 00:21:40,310 equal to some function of x and y. 441 00:21:40,310 --> 00:21:43,010 Y are the algebraic states, and they 442 00:21:43,010 --> 00:21:46,770 satisfy a separate equation, g of x and y and t is equal to 0. 443 00:21:46,770 --> 00:21:49,170 Semi-explicit form. 444 00:21:49,170 --> 00:21:52,720 So with a semi-explicit DAE, the differential index 445 00:21:52,720 --> 00:21:55,880 is defined as the minimum number of differentiations 446 00:21:55,880 --> 00:22:00,270 required to convert the DAE to a system of independent ODEs. 447 00:22:00,270 --> 00:22:01,610 What does that mean? 448 00:22:01,610 --> 00:22:03,890 Means let's take the algebraic equations 449 00:22:03,890 --> 00:22:06,170 and let's take a time derivative of them. 450 00:22:06,170 --> 00:22:07,970 That will give us a new function. 451 00:22:07,970 --> 00:22:10,370 Call it g1. 452 00:22:10,370 --> 00:22:14,090 It's going to be a function of the differential 453 00:22:14,090 --> 00:22:18,020 state, the algebraic state, and the time derivative 454 00:22:18,020 --> 00:22:20,116 of the algebraic state. 455 00:22:20,116 --> 00:22:21,740 It could also be a function of the time 456 00:22:21,740 --> 00:22:23,600 derivative of the differential state, 457 00:22:23,600 --> 00:22:27,310 but we know that in terms of f, which is a function of x and y. 458 00:22:27,310 --> 00:22:30,169 So let's not try putting x dot in there for convenience. 459 00:22:30,169 --> 00:22:31,460 Let's just write out like this. 460 00:22:31,460 --> 00:22:33,466 In principle, we can do this. 461 00:22:33,466 --> 00:22:34,840 Let's take two derivatives of it. 462 00:22:34,840 --> 00:22:35,839 It will be the same way. 463 00:22:35,839 --> 00:22:37,790 It'll give us some function g2, which 464 00:22:37,790 --> 00:22:39,950 is a function of the differential 465 00:22:39,950 --> 00:22:42,410 state, the algebraic state, and the time derivative 466 00:22:42,410 --> 00:22:44,360 of the algebraic state. 467 00:22:44,360 --> 00:22:46,660 Let's take as many of these as we need. 468 00:22:46,660 --> 00:22:48,500 It will give us a system of equations 469 00:22:48,500 --> 00:22:51,860 that we can eventually solve for the time derivatives of all 470 00:22:51,860 --> 00:22:55,939 the algebraic states and convert this DAE to a system of ODEs. 471 00:22:55,939 --> 00:22:57,980 And the question is how many of these derivatives 472 00:22:57,980 --> 00:23:00,140 do I have to take? 473 00:23:00,140 --> 00:23:01,700 Index 1, I need to take one. 474 00:23:01,700 --> 00:23:04,730 Index 2, I need both sets of equations 475 00:23:04,730 --> 00:23:08,780 in order to get something that I can actually solve for dy dt. 476 00:23:08,780 --> 00:23:11,450 Index 3, I get to take three derivatives and higher. 477 00:23:15,890 --> 00:23:18,110 OK, so here's an example. 478 00:23:18,110 --> 00:23:20,090 Let's see if we can work through this. 479 00:23:20,090 --> 00:23:23,930 So we have a DAE system. 480 00:23:23,930 --> 00:23:28,620 Can you calculate the differential index 481 00:23:28,620 --> 00:23:29,310 of this system? 482 00:23:32,296 --> 00:23:32,796 Go ahead. 483 00:23:59,220 --> 00:24:01,546 Not 0. 484 00:24:01,546 --> 00:24:03,290 Not 0. 485 00:24:03,290 --> 00:24:06,430 ODE IVP is-- are DAEs of index 0. 486 00:24:06,430 --> 00:24:08,190 They require no derivatives to generate 487 00:24:08,190 --> 00:24:10,710 a system of ordinary differential equations. 488 00:24:14,651 --> 00:24:15,400 What do you think? 489 00:24:15,400 --> 00:24:17,160 AUDIENCE: [INAUDIBLE] 490 00:24:17,160 --> 00:24:19,010 PROFESSOR: Index 2 you say. 491 00:24:19,010 --> 00:24:20,270 How do you come to index 2? 492 00:24:20,270 --> 00:24:23,050 Why is it index 2? 493 00:24:23,050 --> 00:24:23,720 How'd you do it? 494 00:24:23,720 --> 00:24:25,354 AUDIENCE: [INAUDIBLE] 495 00:24:25,354 --> 00:24:27,520 PROFESSOR: OK, differentiate the algebraic equation. 496 00:24:27,520 --> 00:24:29,215 AUDIENCE: Stick it in [INAUDIBLE] 4. 497 00:24:29,215 --> 00:24:30,640 You get c3 for 0. 498 00:24:30,640 --> 00:24:33,015 Then you have to differentiate [INAUDIBLE] 499 00:24:33,015 --> 00:24:36,770 PROFESSOR: Good so differentiate this equation. 500 00:24:36,770 --> 00:24:39,960 Substitute for c1 dot and c2 dot. 501 00:24:39,960 --> 00:24:44,100 C1 minus c1 will be 0 c2 minus c2 will be 0. 502 00:24:44,100 --> 00:24:46,150 So you have 2c3 equals 0. 503 00:24:46,150 --> 00:24:47,690 That's still an algebraic equation. 504 00:24:47,690 --> 00:24:51,000 We need one more derivative to get a differential equation 505 00:24:51,000 --> 00:24:54,510 for the only algebraic state we have, c3. 506 00:24:54,510 --> 00:24:59,790 So one more derivative will tell us dc3 dt is 0. 507 00:24:59,790 --> 00:25:02,700 So two derivatives to get differential equations. 508 00:25:02,700 --> 00:25:04,221 It's an 2 to DAE. 509 00:25:04,221 --> 00:25:04,720 Sam? 510 00:25:04,720 --> 00:25:07,090 SAM: What if you can get c3 equals 0 511 00:25:07,090 --> 00:25:09,223 and just eliminate it then wouldn't you 512 00:25:09,223 --> 00:25:12,278 write an equivalent some simple equations, but that would be 513 00:25:12,278 --> 00:25:12,778 [INAUDIBLE]? 514 00:25:12,778 --> 00:25:13,867 PROFESSOR: You could. 515 00:25:13,867 --> 00:25:15,450 So we could write an equivalent system 516 00:25:15,450 --> 00:25:18,500 of equations that says instead of this equation, 517 00:25:18,500 --> 00:25:20,490 report c3 equals 0. 518 00:25:20,490 --> 00:25:21,780 That's model formulation. 519 00:25:21,780 --> 00:25:23,430 Here we've given-- we're given a model, 520 00:25:23,430 --> 00:25:26,160 and we're asked what index it is. 521 00:25:26,160 --> 00:25:28,340 We could formulate a different model, 522 00:25:28,340 --> 00:25:30,090 and the model will have a different index. 523 00:25:30,090 --> 00:25:33,600 If I were to replace this equation with c3 equals 0, 524 00:25:33,600 --> 00:25:36,750 what would the index of the DAE system be instead? 525 00:25:36,750 --> 00:25:39,270 Index 1 instead. 526 00:25:39,270 --> 00:25:41,760 I already told you index 2 is harder to solve than index 1, 527 00:25:41,760 --> 00:25:43,680 so if you can formulate an index 1 DAE, 528 00:25:43,680 --> 00:25:45,240 you should probably do it. 529 00:25:45,240 --> 00:25:47,767 But it may not be obvious whether you have or haven't. 530 00:25:47,767 --> 00:25:49,350 So it's an issue of model formulation. 531 00:25:49,350 --> 00:25:51,760 It's a great question. 532 00:25:51,760 --> 00:25:55,650 Does that distinction-- is that clear? 533 00:25:55,650 --> 00:25:57,370 Or is it a little-- 534 00:25:57,370 --> 00:26:01,720 am I not making it clear how this works? 535 00:26:01,720 --> 00:26:02,950 Unresponsive. 536 00:26:02,950 --> 00:26:04,730 OK. 537 00:26:04,730 --> 00:26:06,630 It's OK. 538 00:26:06,630 --> 00:26:09,360 Let's-- there's a generic index 1 example that we can talk 539 00:26:09,360 --> 00:26:10,170 about actually. 540 00:26:10,170 --> 00:26:15,360 So let's take a look at a semi-explicit DAE system. 541 00:26:15,360 --> 00:26:23,240 So we have x dot is f of xy and t and 0 is g of xy and t. 542 00:26:23,240 --> 00:26:28,770 Let's take the time derivative of the algebraic equation. 543 00:26:28,770 --> 00:26:32,010 So the total derivative of g with respect to t 544 00:26:32,010 --> 00:26:33,990 is the Jacobian of g with respect 545 00:26:33,990 --> 00:26:38,500 to x times dx dt plus the Jacobian of g 546 00:26:38,500 --> 00:26:43,260 with respect to y times dy dt plus the derivative of g-- 547 00:26:43,260 --> 00:26:46,510 partial derivative of g with respect to t. 548 00:26:46,510 --> 00:26:47,890 And now let's solve. 549 00:26:47,890 --> 00:26:52,030 Let's push the dy dt term to the other side of the equation, 550 00:26:52,030 --> 00:26:53,830 and let's substitute for dx dt. 551 00:26:53,830 --> 00:26:57,920 We know dx dt is f, so you get this equation here. 552 00:26:57,920 --> 00:27:01,520 If the Jacobian is full rank-- 553 00:27:01,520 --> 00:27:05,840 if dg dy is full rank, then I can invert this matrix, 554 00:27:05,840 --> 00:27:08,420 and I automatically get my system of ODEs 555 00:27:08,420 --> 00:27:09,905 for the algebraic states. 556 00:27:16,150 --> 00:27:18,970 What's problematic about this equation 557 00:27:18,970 --> 00:27:21,289 when dg dy is not full rank? 558 00:27:21,289 --> 00:27:22,330 That should be a partial. 559 00:27:22,330 --> 00:27:23,570 That's sloppy. 560 00:27:23,570 --> 00:27:25,140 I was doing this at 1:00 last night, 561 00:27:25,140 --> 00:27:27,040 so that's why that's there. 562 00:27:27,040 --> 00:27:31,600 What's the-- what's problematic about this root finding problem 563 00:27:31,600 --> 00:27:35,850 here if dg dy or the Jacobian of g with respect to y 564 00:27:35,850 --> 00:27:36,880 is not full rank? 565 00:27:39,760 --> 00:27:40,470 You recall? 566 00:27:45,401 --> 00:27:45,900 If-- 567 00:27:45,900 --> 00:27:48,597 AUDIENCE: [INAUDIBLE] 568 00:27:48,597 --> 00:27:50,930 PROFESSOR: Yes, you won't be able to compute the inverse 569 00:27:50,930 --> 00:27:52,070 to get your dy dt's. 570 00:27:52,070 --> 00:27:52,962 That's true. 571 00:27:52,962 --> 00:27:54,920 Remember this Jacobian-- if it's not full rank, 572 00:27:54,920 --> 00:27:57,950 it's determinant is 0. 573 00:27:57,950 --> 00:28:00,200 There is no inverse of this thing. 574 00:28:00,200 --> 00:28:03,350 We talked once about the systems of nonlinear equations 575 00:28:03,350 --> 00:28:08,120 and locally unique solutions, the so-called inverse function 576 00:28:08,120 --> 00:28:09,130 theorem. 577 00:28:09,130 --> 00:28:13,010 We can only guarantee there is a locally unique solution when 578 00:28:13,010 --> 00:28:15,040 I can invert this thing. 579 00:28:15,040 --> 00:28:17,000 That if I got in really close to the solution, 580 00:28:17,000 --> 00:28:18,890 it looked like a linear system of equations. 581 00:28:18,890 --> 00:28:21,820 That linear system of equations had a unique solution 582 00:28:21,820 --> 00:28:22,570 associated within. 583 00:28:22,570 --> 00:28:23,570 And everything was good. 584 00:28:23,570 --> 00:28:25,640 We were happy. 585 00:28:25,640 --> 00:28:28,250 So index 1 DAEs have that property. 586 00:28:28,250 --> 00:28:30,700 If I wanted to solve this equation for y, 587 00:28:30,700 --> 00:28:34,739 there will be a locally unique solution for y. 588 00:28:34,739 --> 00:28:36,530 Higher index DAEs don't have that property. 589 00:28:36,530 --> 00:28:39,970 We know that's a sort of unhappy generic circumstance 590 00:28:39,970 --> 00:28:42,370 to be in if we have to solve this equation. 591 00:28:42,370 --> 00:28:44,710 It can be hard to find those roots using, 592 00:28:44,710 --> 00:28:47,470 say, Newton-Raphson because the Jacobian [? e ?] 593 00:28:47,470 --> 00:28:49,870 would need to become singular during the root finding 594 00:28:49,870 --> 00:28:52,570 procedure. 595 00:28:52,570 --> 00:28:54,580 This is connected intimately to what 596 00:28:54,580 --> 00:28:57,280 we did for systems of nonlinear equations. 597 00:28:57,280 --> 00:28:59,410 So for an index 1 DAE, you can show its index 598 00:28:59,410 --> 00:29:04,880 1 because this Jacobian is invertible. 599 00:29:04,880 --> 00:29:07,480 Let's do some more examples on differential index. 600 00:29:07,480 --> 00:29:09,230 This seems to be the most important thing. 601 00:29:09,230 --> 00:29:10,790 If I have a model, what's its index 602 00:29:10,790 --> 00:29:13,040 because its index is going to tell me how sensitive it 603 00:29:13,040 --> 00:29:16,840 is to small perturbations. 604 00:29:16,840 --> 00:29:17,991 So here's a model. 605 00:29:17,991 --> 00:29:19,740 It's the same as the model you saw before, 606 00:29:19,740 --> 00:29:21,240 but now it's two mixing tanks. 607 00:29:21,240 --> 00:29:26,770 So here comes an inlet flow q1 carrying concentration c1. 608 00:29:26,770 --> 00:29:30,540 Out comes the same flow q1 carrying concentration c2. 609 00:29:30,540 --> 00:29:33,540 The mixer has volume v1. 610 00:29:33,540 --> 00:29:36,630 And then I blend this with some more 611 00:29:36,630 --> 00:29:40,640 of whatever the solute is at a flow rate of q2 612 00:29:40,640 --> 00:29:43,140 and a concentration c3. 613 00:29:43,140 --> 00:29:44,910 And those both go into this tank, 614 00:29:44,910 --> 00:29:48,910 and they come out concentration c4, flow rate q1 plus q2, 615 00:29:48,910 --> 00:29:52,076 and this has volume v2. 616 00:29:52,076 --> 00:29:53,320 Does that look good? 617 00:29:53,320 --> 00:29:54,760 Well-posed model? 618 00:29:54,760 --> 00:29:58,525 OK, so here's your material balances on the mixers. 619 00:29:58,525 --> 00:30:00,400 And I'm going to give you different algebraic 620 00:30:00,400 --> 00:30:01,450 constraints. 621 00:30:01,450 --> 00:30:06,120 So this is a problem now where we say measure c1 and c3. 622 00:30:09,400 --> 00:30:14,175 And we want to try to predict c2 and c4. 623 00:30:14,175 --> 00:30:16,120 Can you figure out the differential index 624 00:30:16,120 --> 00:30:17,160 of this DAE system? 625 00:30:26,640 --> 00:30:28,010 Feel free to talk to each other. 626 00:30:28,010 --> 00:30:28,510 It's OK. 627 00:30:59,210 --> 00:31:01,820 OK, the differential index is 1. 628 00:31:01,820 --> 00:31:03,410 This one's the easy one right? 629 00:31:03,410 --> 00:31:06,185 Just take one derivative of the algebraic equations. 630 00:31:06,185 --> 00:31:11,000 Now, I've got dc2 dt, dc4 dt, dc1 dt, dc3 dt. 631 00:31:11,000 --> 00:31:16,560 Differential index of 1 is easy to solve as an ODE IVP. 632 00:31:16,560 --> 00:31:18,800 This is the natural problem, too. 633 00:31:18,800 --> 00:31:22,640 It's useful to think about the inputs c1, c3 and asking 634 00:31:22,640 --> 00:31:23,930 about what the outputs are. 635 00:31:23,930 --> 00:31:26,765 Physically, this problem seems easier in nature. 636 00:31:32,570 --> 00:31:35,280 OK let's change it now. 637 00:31:35,280 --> 00:31:36,649 Same problem. 638 00:31:36,649 --> 00:31:38,440 But let's change the algebraic constraints. 639 00:31:38,440 --> 00:31:42,040 Now, the algebraic constraint is I measure c3 and c4, 640 00:31:42,040 --> 00:31:43,815 And I want to predict c1 and c2. 641 00:31:46,970 --> 00:31:48,450 What's the differential index here? 642 00:32:53,600 --> 00:32:54,100 Yes. 643 00:32:54,100 --> 00:32:54,980 AUDIENCE: [INAUDIBLE] 644 00:32:54,980 --> 00:32:56,250 PROFESSOR: Sure. 645 00:32:56,250 --> 00:32:58,600 AUDIENCE: Why was this a differential index of 1 646 00:32:58,600 --> 00:33:01,890 if you have to take a derivative for c1 and this equation? 647 00:33:01,890 --> 00:33:03,800 PROFESSOR: Oh, good question. 648 00:33:03,800 --> 00:33:06,810 I'm going to push the slide back, and I apologize for this. 649 00:33:06,810 --> 00:33:07,979 And I'll go back for it. 650 00:33:07,979 --> 00:33:10,270 So the question was why was this differential index one 651 00:33:10,270 --> 00:33:13,100 if I had to take a derivative of both this equation 652 00:33:13,100 --> 00:33:14,770 and that equation. 653 00:33:14,770 --> 00:33:17,260 So the way to think about this conceptually 654 00:33:17,260 --> 00:33:21,789 is I took a derivative of the algebraic equations, 655 00:33:21,789 --> 00:33:23,080 the set of algebraic equations. 656 00:33:23,080 --> 00:33:29,366 It was one time derivative of a vector valued function. 657 00:33:29,366 --> 00:33:30,990 I can see where the ambiguity is there. 658 00:33:30,990 --> 00:33:32,240 So it's important to be clear. 659 00:33:32,240 --> 00:33:35,230 I took one time derivative over the equations that 660 00:33:35,230 --> 00:33:37,400 prescribe the algebraic states. 661 00:33:41,260 --> 00:33:42,280 That's the distinction. 662 00:33:42,280 --> 00:33:43,780 So this is index 1 because I needed 663 00:33:43,780 --> 00:33:55,600 one time derivative of this type, dg dt of the entire set. 664 00:33:55,600 --> 00:33:56,100 Here we go. 665 00:33:56,100 --> 00:33:57,550 We're on this one. 666 00:33:57,550 --> 00:33:58,340 Prescribe c3. 667 00:33:58,340 --> 00:33:59,270 Prescribe c4. 668 00:33:59,270 --> 00:34:00,860 What is the differential index now? 669 00:34:06,920 --> 00:34:07,712 What do you think? 670 00:34:07,712 --> 00:34:08,420 We got an answer? 671 00:34:08,420 --> 00:34:10,670 [INTERPOSING VOICES] 672 00:34:10,670 --> 00:34:11,449 I hear 2. 673 00:34:11,449 --> 00:34:13,412 I hear 3. 674 00:34:13,412 --> 00:34:14,860 I saw 3 like this. 675 00:34:14,860 --> 00:34:16,719 Almost mistook it for a 2, but that's 3. 676 00:34:16,719 --> 00:34:17,590 Yes. 677 00:34:17,590 --> 00:34:18,337 It's 3. 678 00:34:18,337 --> 00:34:20,170 Let's see if we can work through why it's 3. 679 00:34:23,330 --> 00:34:26,300 So I take a derivative of the algebraic equations, 680 00:34:26,300 --> 00:34:27,260 one derivative. 681 00:34:27,260 --> 00:34:28,834 And after taking that derivative, 682 00:34:28,834 --> 00:34:30,500 I'll get a differential equation for one 683 00:34:30,500 --> 00:34:31,730 of the algebraic states. 684 00:34:31,730 --> 00:34:33,889 So c3 is taking care of. 685 00:34:33,889 --> 00:34:37,010 Now, I'm in the hunt for c1. 686 00:34:37,010 --> 00:34:38,780 So I have-- after that first derivative, 687 00:34:38,780 --> 00:34:42,320 I have an equation dc4 dt is gamma 2. 688 00:34:42,320 --> 00:34:45,530 And I know dc4 dt, so I drop that in, 689 00:34:45,530 --> 00:34:49,670 and I get an algebraic equation relating the derivative 690 00:34:49,670 --> 00:34:54,250 of gamma 2 to c2, c3, and c4. 691 00:34:54,250 --> 00:34:56,210 This new algebraic equation. 692 00:34:56,210 --> 00:34:59,120 This is like the g1 that I prescribed 693 00:34:59,120 --> 00:35:00,100 in the generic scene. 694 00:35:00,100 --> 00:35:02,725 I take a derivative of it again because I'm in the hunt for c1. 695 00:35:02,725 --> 00:35:04,760 So I take the derivative of this new equation, 696 00:35:04,760 --> 00:35:07,160 and I'll get a derivative of c2, a derivative of c3, 697 00:35:07,160 --> 00:35:08,150 and a derivative of c4. 698 00:35:08,150 --> 00:35:09,620 And I know all those. 699 00:35:09,620 --> 00:35:11,140 I know the derivative of c2. 700 00:35:11,140 --> 00:35:13,140 I know the derivative of c4. 701 00:35:13,140 --> 00:35:16,130 And I know the derivative of c3 from the previous level 702 00:35:16,130 --> 00:35:17,040 in the hierarchy. 703 00:35:17,040 --> 00:35:18,490 I figured that out already. 704 00:35:18,490 --> 00:35:19,810 So it's two derivatives. 705 00:35:19,810 --> 00:35:22,460 Still hasn't gotten me a time derivative of c1. 706 00:35:22,460 --> 00:35:23,960 But when I make those substitutions, 707 00:35:23,960 --> 00:35:27,390 I'll get an algebraic equation in terms of c1. 708 00:35:27,390 --> 00:35:29,780 And I can take one more derivative, 709 00:35:29,780 --> 00:35:31,740 and that will give me a dc1 dt. 710 00:35:31,740 --> 00:35:34,430 And I'll have an od for c1, c2, c3, and c4. 711 00:35:34,430 --> 00:35:36,900 So its differential index 3. 712 00:35:41,050 --> 00:35:42,340 It makes sense sort of. 713 00:35:42,340 --> 00:35:47,630 I'm taking a measurement of an output way down the line here. 714 00:35:47,630 --> 00:35:51,120 And I'm trying to predict what the input was 715 00:35:51,120 --> 00:35:51,966 in the first place. 716 00:35:51,966 --> 00:35:53,340 It's easy to imagine that there's 717 00:35:53,340 --> 00:35:56,220 a huge amount of sensitivity in that calculation. 718 00:36:00,570 --> 00:36:02,182 Here's another one. 719 00:36:02,182 --> 00:36:03,640 What's the differential index here? 720 00:36:10,310 --> 00:36:14,790 I'm now prescribing c1 and c2, and I 721 00:36:14,790 --> 00:36:16,060 want you to tell me c3 and c4. 722 00:36:16,060 --> 00:36:24,870 [INTERPOSING VOICES] 723 00:36:24,870 --> 00:36:25,670 Oh, good. 724 00:36:25,670 --> 00:36:26,870 Somebody noticed early on. 725 00:36:26,870 --> 00:36:29,300 So you say it's not possible. 726 00:36:29,300 --> 00:36:30,259 Right. 727 00:36:30,259 --> 00:36:31,175 Why isn't it possible? 728 00:36:31,175 --> 00:36:33,530 AUDIENCE: You take derivatives of both of those-- 729 00:36:33,530 --> 00:36:34,155 PROFESSOR: Yes. 730 00:36:34,155 --> 00:36:37,202 AUDIENCE: [INAUDIBLE] you can't isolate c3-- 731 00:36:37,202 --> 00:36:38,980 PROFESSOR: Right. 732 00:36:38,980 --> 00:36:42,040 We'll never be able to isolate a derivative of c3 here. 733 00:36:42,040 --> 00:36:44,230 There's actually there's something wrong physically 734 00:36:44,230 --> 00:36:45,070 with this problem. 735 00:36:47,830 --> 00:36:52,780 Yes, somehow I'm supposed to measure c1 and c2 736 00:36:52,780 --> 00:36:56,740 and use that to predict c3 and c4. 737 00:36:56,740 --> 00:36:59,592 I don't know c4, and c3 is an input. 738 00:36:59,592 --> 00:37:00,550 I don't know it either. 739 00:37:00,550 --> 00:37:05,260 How do I figure out an input when I don't know the output? 740 00:37:05,260 --> 00:37:06,430 It's impossible. 741 00:37:06,430 --> 00:37:10,554 So this model is flawed. 742 00:37:10,554 --> 00:37:12,220 We can formulate any model we want, pick 743 00:37:12,220 --> 00:37:15,920 any set of these variables to prescribe algebraically, 744 00:37:15,920 --> 00:37:20,050 but not all of them are going to admit solutions or make sense. 745 00:37:20,050 --> 00:37:21,820 There's no number of derivatives that's 746 00:37:21,820 --> 00:37:25,024 going to give us a closed system of ODEs. 747 00:37:25,024 --> 00:37:26,440 This again comes back to the point 748 00:37:26,440 --> 00:37:29,832 that really DAEs is all about model formulation. 749 00:37:29,832 --> 00:37:31,540 There are lots of good numerical methods. 750 00:37:31,540 --> 00:37:34,090 They work like the numerical methods for ODE IVPs. 751 00:37:34,090 --> 00:37:37,150 What's important is getting consistent models down, 752 00:37:37,150 --> 00:37:39,850 models that, for example, have solutions. 753 00:37:39,850 --> 00:37:41,650 You might feed this to a DAE solver 754 00:37:41,650 --> 00:37:45,460 and get nonsense because I can almost freely prescribe 755 00:37:45,460 --> 00:37:48,790 what c3 is, and I'll get some answer for c4. 756 00:37:48,790 --> 00:37:50,450 Any c3 can give me an answer here. 757 00:37:50,450 --> 00:37:53,740 So who knows whether this numerical solver is sensitive 758 00:37:53,740 --> 00:37:56,999 or not to this particular pathology in the problem 759 00:37:56,999 --> 00:37:57,790 we've written down. 760 00:38:02,080 --> 00:38:04,710 We got to speed up just a little bit, and I'm sorry for that. 761 00:38:04,710 --> 00:38:06,540 So we looked at index 1. 762 00:38:06,540 --> 00:38:07,590 That was stirred tank 1. 763 00:38:07,590 --> 00:38:09,450 Index 2, that was stirred tank 2. 764 00:38:09,450 --> 00:38:11,439 Index 3, that was stirred tank 3. 765 00:38:11,439 --> 00:38:13,230 And one thing we saw when we looked at them 766 00:38:13,230 --> 00:38:15,360 was this index 1 solution. 767 00:38:15,360 --> 00:38:17,970 It had pieces that were proportional 768 00:38:17,970 --> 00:38:22,470 to this forcing function gamma, this prescribed function 769 00:38:22,470 --> 00:38:26,710 in our system of equations and proportional to its derivative. 770 00:38:26,710 --> 00:38:29,510 So if gamma is jumping around, well, c1 will jump around. 771 00:38:29,510 --> 00:38:32,220 C2 is going to be smoothed out version of gamma 772 00:38:32,220 --> 00:38:34,630 because it depends on its integral. 773 00:38:34,630 --> 00:38:37,740 It's not very sensitive, kind of like an ODE IVP. 774 00:38:37,740 --> 00:38:40,260 Doesn't really show any more sensitivity than an ODE IVP 775 00:38:40,260 --> 00:38:41,250 does. 776 00:38:41,250 --> 00:38:43,350 The dynamics are different for index 2, though. 777 00:38:43,350 --> 00:38:46,050 When we look at the solution of the index 2 DAE, 778 00:38:46,050 --> 00:38:51,390 we found c2 goes like gamma, but c1 has to go like gamma dot. 779 00:38:51,390 --> 00:38:54,549 So that's pretty sensitive if you're making a measurement, 780 00:38:54,549 --> 00:38:56,090 and there's noise in the measurement. 781 00:38:56,090 --> 00:38:58,170 How do you even know what this derivative is? 782 00:38:58,170 --> 00:39:01,020 So c1 is wildly bouncing around. 783 00:39:01,020 --> 00:39:03,060 Our prediction of what c1 is is-- 784 00:39:03,060 --> 00:39:05,162 it's not going to be very good. 785 00:39:05,162 --> 00:39:09,330 Our solution of DAE example 3 when it told us c1 786 00:39:09,330 --> 00:39:11,730 was related to the second derivative 787 00:39:11,730 --> 00:39:15,332 of the forcing function, c2 to the first and c3 788 00:39:15,332 --> 00:39:17,040 was proportional to the forcing function. 789 00:39:17,040 --> 00:39:19,100 So this is hugely sensitive to changes 790 00:39:19,100 --> 00:39:20,100 in the forcing function. 791 00:39:20,100 --> 00:39:24,030 So the higher the index goes, the greater the sensitivity 792 00:39:24,030 --> 00:39:25,425 to perturbations in the system. 793 00:39:29,840 --> 00:39:31,540 Here's another simple example. 794 00:39:31,540 --> 00:39:33,073 You have all that data, actually, in your paper, 795 00:39:33,073 --> 00:39:35,614 so I won't ask you to do this one given our time constraints. 796 00:39:35,614 --> 00:39:39,040 But here's-- mechanical systems that have constraints are often 797 00:39:39,040 --> 00:39:41,260 indexed 3 DAEs it turns out. 798 00:39:41,260 --> 00:39:43,390 You can show this one is an index 3 DAE. 799 00:39:43,390 --> 00:39:47,170 This is the case of a pendulum swinging back and forth. 800 00:39:47,170 --> 00:39:50,740 So it's a-- change in position is its velocity. 801 00:39:50,740 --> 00:39:53,230 Its acceleration balances gravity, 802 00:39:53,230 --> 00:39:56,080 and there's some arm that holds the pendulum in place. 803 00:39:56,080 --> 00:39:58,420 We can imagine that that arm has some tension in it. 804 00:39:58,420 --> 00:40:01,510 It acts like a spring with a spring constant 805 00:40:01,510 --> 00:40:05,260 that changes in time in order to hold the pendulum at a fixed 806 00:40:05,260 --> 00:40:07,150 distance from its center. 807 00:40:07,150 --> 00:40:11,300 So two differential equations, one algebraic equation. 808 00:40:11,300 --> 00:40:13,660 This is a differential algebraic equation. 809 00:40:13,660 --> 00:40:16,270 You can imagine lots of mechanical systems 810 00:40:16,270 --> 00:40:18,130 work in exactly this way. 811 00:40:18,130 --> 00:40:20,490 They can only move along prescribed paths. 812 00:40:20,490 --> 00:40:23,679 They're constrained in how far they can stretch or go. 813 00:40:23,679 --> 00:40:24,220 I don't know. 814 00:40:24,220 --> 00:40:26,310 If you're trying to design a robot or something, 815 00:40:26,310 --> 00:40:29,180 DAEs are pretty important, and they're all of index 3 type 816 00:40:29,180 --> 00:40:30,370 it turns out. 817 00:40:30,370 --> 00:40:33,730 So differential variables here are x and v, 818 00:40:33,730 --> 00:40:36,700 and the algebraic variable is the spring constant, k, 819 00:40:36,700 --> 00:40:40,690 which has got to adjust dynamically in time in order 820 00:40:40,690 --> 00:40:41,460 to-- 821 00:40:41,460 --> 00:40:43,150 got to get away from the speaker-- 822 00:40:43,150 --> 00:40:44,710 adjust dynamically in time in order 823 00:40:44,710 --> 00:40:48,310 to provide just the right stiffness to keep this going 824 00:40:48,310 --> 00:40:51,170 around on a circular orbit. 825 00:40:51,170 --> 00:40:54,280 So when-- let's suppose I start with my pendulum down, 826 00:40:54,280 --> 00:40:56,170 and it's not moving, then this has 827 00:40:56,170 --> 00:40:59,530 to be just stiff enough to balance gravity. 828 00:40:59,530 --> 00:41:01,916 The pendulum is 90 degrees, and it's not moving. 829 00:41:01,916 --> 00:41:03,040 I don't need any stiffness. 830 00:41:03,040 --> 00:41:04,570 K and just be 0. 831 00:41:04,570 --> 00:41:06,640 There's no forces to balance. 832 00:41:06,640 --> 00:41:08,350 The pendulum is swinging around, and I 833 00:41:08,350 --> 00:41:10,570 have to be stiff enough to balance gravity 834 00:41:10,570 --> 00:41:14,500 where I am and also counteract any sort of centripetal 835 00:41:14,500 --> 00:41:16,180 [? acceleration. ?] So this k has 836 00:41:16,180 --> 00:41:18,970 got to be wildly fluctuating. 837 00:41:18,970 --> 00:41:21,280 If I was trying to control that k somehow 838 00:41:21,280 --> 00:41:25,000 to give this system these particular dynamics, 839 00:41:25,000 --> 00:41:29,270 you might imagine it's very difficult to do. 840 00:41:29,270 --> 00:41:31,340 So there's your differential variables. 841 00:41:31,340 --> 00:41:32,780 Here's your algebraic variable. 842 00:41:32,780 --> 00:41:34,437 You take a derivative of this equation, 843 00:41:34,437 --> 00:41:36,020 you'll get a constraint that tells you 844 00:41:36,020 --> 00:41:38,554 the velocity is orthogonal to the position. 845 00:41:38,554 --> 00:41:39,220 Of course it is. 846 00:41:39,220 --> 00:41:41,330 I'm on a circular trajectory. 847 00:41:41,330 --> 00:41:43,430 You take a derivative of this equation 848 00:41:43,430 --> 00:41:44,990 now, this algebraic equation. 849 00:41:44,990 --> 00:41:46,823 You'll get another constraint that gives you 850 00:41:46,823 --> 00:41:49,760 some relationship between k, x, and v 851 00:41:49,760 --> 00:41:51,860 but not a differential equation for k. 852 00:41:51,860 --> 00:41:54,456 Take another derivative of this algebraic equation, 853 00:41:54,456 --> 00:41:56,330 and you'll get a differential equation for k. 854 00:41:56,330 --> 00:41:58,660 So in principle, I could formulate a system of ODEs 855 00:41:58,660 --> 00:42:02,980 with dk dt, dx dt, dv dt will be equivalent to this. 856 00:42:02,980 --> 00:42:05,596 It took three derivatives to do this, though. 857 00:42:05,596 --> 00:42:07,220 So it's an incredibly sensitive system. 858 00:42:07,220 --> 00:42:08,460 It's index 3 in nature. 859 00:42:11,510 --> 00:42:13,030 So if I transform to this equivalent 860 00:42:13,030 --> 00:42:14,500 set of ODEs, the problem-- 861 00:42:14,500 --> 00:42:17,140 we discussed this earlier-- is that the solutions may 862 00:42:17,140 --> 00:42:21,197 drift away from the initial set of constraints. 863 00:42:21,197 --> 00:42:22,780 The solutions also need the right sort 864 00:42:22,780 --> 00:42:24,210 of initial conditions. 865 00:42:24,210 --> 00:42:31,400 Here I know that at time 0, c3 better be gamma 0. 866 00:42:31,400 --> 00:42:33,860 I know that in the actual solution to this problem, 867 00:42:33,860 --> 00:42:36,440 c3 dot was gamma dot. 868 00:42:36,440 --> 00:42:41,570 So at time 0, c2 better be gamma dot of 0. 869 00:42:41,570 --> 00:42:44,840 And at time 0, c1 better be gamma dot dot of 0. 870 00:42:44,840 --> 00:42:48,160 And if it's not, then there's going 871 00:42:48,160 --> 00:42:50,800 to be some step jump or some strange behavior 872 00:42:50,800 --> 00:42:53,210 in the solution to this ODE system. 873 00:42:53,210 --> 00:42:55,675 So I have to choose the right sort of initial conditions. 874 00:42:55,675 --> 00:42:57,850 If those initial conditions aren't chosen correctly, 875 00:42:57,850 --> 00:42:59,090 then I'm not going to get the right solution. 876 00:42:59,090 --> 00:43:01,990 I won't be constrained to the manifold of solutions that's 877 00:43:01,990 --> 00:43:03,655 given by the original DAE. 878 00:43:06,446 --> 00:43:08,070 So here's some things you need to know. 879 00:43:08,070 --> 00:43:12,680 Index 1, semi-explicit DAEs can be safely handled in MATLAB. 880 00:43:12,680 --> 00:43:16,670 So ode15s, ode23t, they have the ability 881 00:43:16,670 --> 00:43:18,870 to take an input of mass matrix. 882 00:43:18,870 --> 00:43:21,830 We discussed mass matrix last time. 883 00:43:21,830 --> 00:43:24,055 And a right hand side to the system of equations 884 00:43:24,055 --> 00:43:27,800 in f of x and t and solve it just 885 00:43:27,800 --> 00:43:30,380 like all the other solvers, all the other ODE IVP 886 00:43:30,380 --> 00:43:32,540 solvers in MATLAB. 887 00:43:32,540 --> 00:43:35,900 If it's not index 1, MATLAB can catch it. 888 00:43:35,900 --> 00:43:39,110 So it'll try to look at the Jacobian of f with respect 889 00:43:39,110 --> 00:43:40,970 to the variables and the mass matrix 890 00:43:40,970 --> 00:43:45,680 and to infer from that what the differential index is. 891 00:43:45,680 --> 00:43:48,759 And it will tell you often times-- 892 00:43:48,759 --> 00:43:50,550 I think depends which package you're using. 893 00:43:50,550 --> 00:43:52,910 But if you're using certain packages, it'll tell you. 894 00:43:52,910 --> 00:43:55,170 It's not index 1oe-- 895 00:43:55,170 --> 00:43:56,360 DAE. 896 00:43:56,360 --> 00:43:57,770 We can't handle this. 897 00:43:57,770 --> 00:44:01,220 The methods built into it aren't suitable. 898 00:44:01,220 --> 00:44:03,860 For generic DAEs, there are specific DAE solvers. 899 00:44:03,860 --> 00:44:06,020 So something like SUNDIALS or DAEPACK-- 900 00:44:06,020 --> 00:44:09,020 this is Professor Barton's software actually-- 901 00:44:09,020 --> 00:44:13,050 that are designed to handle DAEs up to some certain index 902 00:44:13,050 --> 00:44:13,580 instead. 903 00:44:13,580 --> 00:44:15,870 They're built to reliably solve these problems. 904 00:44:15,870 --> 00:44:18,500 So we have robots with constraints on them. 905 00:44:18,500 --> 00:44:21,779 We solve all sorts of problems for chemical process systems 906 00:44:21,779 --> 00:44:22,820 that are of higher index. 907 00:44:22,820 --> 00:44:25,589 And the way we do it is using specific numerical methods. 908 00:44:25,589 --> 00:44:27,380 They're based on the same sorts of schemes, 909 00:44:27,380 --> 00:44:30,020 backwards difference formulas, but careful analysis 910 00:44:30,020 --> 00:44:33,860 of the equations, which are in general unstructured, 911 00:44:33,860 --> 00:44:35,360 the software doesn't know beforehand 912 00:44:35,360 --> 00:44:38,630 what those equations look like, but a careful analysis of them 913 00:44:38,630 --> 00:44:41,180 in order to figure out how to minimize numerical errors 914 00:44:41,180 --> 00:44:45,620 and stay on the correct solution manifold. 915 00:44:45,620 --> 00:44:47,930 As input to these, though, we have 916 00:44:47,930 --> 00:44:49,440 to give initial conditions. 917 00:44:49,440 --> 00:44:52,675 So that's going to be the last thing that we talk about here. 918 00:44:52,675 --> 00:44:54,050 And those initial conditions have 919 00:44:54,050 --> 00:44:56,750 to be prescribed what we call consistently, 920 00:44:56,750 --> 00:45:00,320 or we can get numerical errors, or the software can just 921 00:45:00,320 --> 00:45:03,664 quit and throw an error, or it can run off 922 00:45:03,664 --> 00:45:06,080 on some other solution manifold that doesn't look anything 923 00:45:06,080 --> 00:45:08,060 like the solution we're after. 924 00:45:08,060 --> 00:45:11,360 The pendulum is an interesting example to think about. 925 00:45:11,360 --> 00:45:13,220 So you might say, well, what do I want 926 00:45:13,220 --> 00:45:14,712 to know about this pendulum. 927 00:45:14,712 --> 00:45:16,670 Maybe I want to know where it started initially 928 00:45:16,670 --> 00:45:20,030 and what its velocity initially was. 929 00:45:20,030 --> 00:45:23,240 Can I prescribe the initial position of the pendulum 930 00:45:23,240 --> 00:45:25,680 arbitrarily? 931 00:45:25,680 --> 00:45:27,358 What do you think yes or no? 932 00:45:27,358 --> 00:45:28,186 AUDIENCE: No. 933 00:45:28,186 --> 00:45:28,851 PROFESSOR: No. 934 00:45:28,851 --> 00:45:29,350 Why not? 935 00:45:29,350 --> 00:45:29,920 Why no? 936 00:45:29,920 --> 00:45:33,010 AUDIENCE: [INAUDIBLE] 937 00:45:33,010 --> 00:45:34,160 PROFESSOR: Good, yes. 938 00:45:34,160 --> 00:45:36,050 We know it's got a fixed arm length. 939 00:45:36,050 --> 00:45:37,069 Can't pick any position. 940 00:45:37,069 --> 00:45:38,110 That's not going to work. 941 00:45:38,110 --> 00:45:40,790 It's got to satisfy that algebraic constraint. 942 00:45:40,790 --> 00:45:45,110 Can I specify its velocity arbitrarily? 943 00:45:45,110 --> 00:45:46,130 What you think? 944 00:45:46,130 --> 00:45:48,426 Can I give it any initial velocity that mass 945 00:45:48,426 --> 00:45:49,550 at the end of the pendulum? 946 00:45:49,550 --> 00:45:50,350 AUDIENCE: No. 947 00:45:50,350 --> 00:45:51,001 PROFESSOR: No. 948 00:45:51,001 --> 00:45:51,500 Why not? 949 00:45:51,500 --> 00:45:52,448 AUDIENCE: [INAUDIBLE] 950 00:45:52,448 --> 00:45:53,870 PROFESSOR: What's that? 951 00:45:53,870 --> 00:45:56,714 AUDIENCE: [INAUDIBLE] 952 00:45:56,714 --> 00:45:58,540 PROFESSOR: Well, explain that to me, 953 00:45:58,540 --> 00:46:02,170 though, because I had equations for the velocity 954 00:46:02,170 --> 00:46:04,870 and for that acceleration and then 955 00:46:04,870 --> 00:46:08,390 an equation that told me that the length of this moment arm 956 00:46:08,390 --> 00:46:10,550 had a certain length associated with it. 957 00:46:10,550 --> 00:46:13,899 So why does that constrain the velocity? 958 00:46:13,899 --> 00:46:14,440 You're right. 959 00:46:14,440 --> 00:46:15,160 It does. 960 00:46:15,160 --> 00:46:16,326 I'm not saying you're wrong. 961 00:46:16,326 --> 00:46:17,370 You're right, but why? 962 00:46:17,370 --> 00:46:19,161 Those initial three equations don't tell me 963 00:46:19,161 --> 00:46:21,917 anything about the initial velocity was supposed to be. 964 00:46:21,917 --> 00:46:23,500 We know something about the trajectory 965 00:46:23,500 --> 00:46:25,083 this is supposed to sweep out, though. 966 00:46:25,083 --> 00:46:27,397 It's supposed to sweep out a circular trajectory. 967 00:46:27,397 --> 00:46:29,230 And the velocity in that circular trajectory 968 00:46:29,230 --> 00:46:32,470 is always going to be orthogonal to the moment time. 969 00:46:32,470 --> 00:46:35,770 And that equation actually-- it appeared. 970 00:46:35,770 --> 00:46:39,080 It appeared when we took a derivative of the constraint. 971 00:46:39,080 --> 00:46:41,950 So this was a hidden constraint. 972 00:46:41,950 --> 00:46:42,726 It's secret. 973 00:46:42,726 --> 00:46:45,100 You didn't know it was there until you took a derivative, 974 00:46:45,100 --> 00:46:47,980 and you tried to go and figure out what index this was. 975 00:46:47,980 --> 00:46:50,920 Then you discovered velocity and position or orthogonal 976 00:46:50,920 --> 00:46:52,070 to each other. 977 00:46:52,070 --> 00:46:55,150 My initial condition should respect that. 978 00:46:55,150 --> 00:46:56,050 How can they not? 979 00:46:56,050 --> 00:46:57,883 The initial conditions have to be a solution 980 00:46:57,883 --> 00:47:00,710 to the equations as well. 981 00:47:00,710 --> 00:47:05,110 So now can the initial stiffness be specified arbitrarily? 982 00:47:05,110 --> 00:47:06,590 No. 983 00:47:06,590 --> 00:47:08,650 There was also an algebraic equation 984 00:47:08,650 --> 00:47:10,600 that popped up at some point right here when 985 00:47:10,600 --> 00:47:12,433 I took the second derivative, which tells me 986 00:47:12,433 --> 00:47:13,990 how the stiffness has to be related 987 00:47:13,990 --> 00:47:15,520 to position and velocities. 988 00:47:15,520 --> 00:47:20,410 It's also hidden inside the structure of this DAE. 989 00:47:20,410 --> 00:47:22,866 So this-- again goes to model formulation, 990 00:47:22,866 --> 00:47:25,240 what are the initial conditions that I have to prescribe? 991 00:47:25,240 --> 00:47:26,860 If I prescribe the right ones, I'll 992 00:47:26,860 --> 00:47:30,670 get a solution that matches the dynamics of the problem I 993 00:47:30,670 --> 00:47:32,010 was actually interested in. 994 00:47:32,010 --> 00:47:33,970 If I prescribe them incorrectly, who 995 00:47:33,970 --> 00:47:35,830 knows what's going to result? 996 00:47:35,830 --> 00:47:38,352 We can't really predict. 997 00:47:38,352 --> 00:47:40,060 Depends on the solver we're working with. 998 00:47:42,820 --> 00:47:46,310 Here's a formal way to think of it. 999 00:47:46,310 --> 00:47:48,150 So usually in an initial value problem, 1000 00:47:48,150 --> 00:47:51,380 usually we want to know what's-- 1001 00:47:51,380 --> 00:47:53,510 for these first order difference-- 1002 00:47:53,510 --> 00:47:57,290 differential equations, what's the initial value of the state, 1003 00:47:57,290 --> 00:48:00,230 and what's the initial value of its derivative. 1004 00:48:00,230 --> 00:48:05,559 That completely specifies what's going on at the beginning. 1005 00:48:05,559 --> 00:48:07,100 One way to think about these problems 1006 00:48:07,100 --> 00:48:19,420 is I've got this equation that I want to solve. 1007 00:48:19,420 --> 00:48:22,381 I could think in some sense that x and x dot 1008 00:48:22,381 --> 00:48:23,630 are independent of each other. 1009 00:48:23,630 --> 00:48:30,840 I could always write this is something like an equation 1010 00:48:30,840 --> 00:48:33,150 for f of x and z. 1011 00:48:33,150 --> 00:48:35,340 And then I could add some extra information 1012 00:48:35,340 --> 00:48:42,430 that tells me actually x dot is the same as z. 1013 00:48:42,430 --> 00:48:46,040 So initially, I would like to know the values of x and z, 1014 00:48:46,040 --> 00:48:47,340 that z is x dot. 1015 00:48:47,340 --> 00:48:51,460 So I'd like to know initially the values of x dot and x. 1016 00:48:51,460 --> 00:48:54,250 And tell me where I start, and these 1017 00:48:54,250 --> 00:48:57,220 need to be consistent with the initial value problem 1018 00:48:57,220 --> 00:48:58,030 that I specified. 1019 00:48:58,030 --> 00:49:00,210 If I'm doing an index 0, an ODE IVP, at least 1020 00:49:00,210 --> 00:49:01,210 be consistent with this. 1021 00:49:01,210 --> 00:49:04,060 So if you tell me x naught, that I 1022 00:49:04,060 --> 00:49:06,280 should put x into this equation, then that 1023 00:49:06,280 --> 00:49:08,770 should tell me what x dot is. 1024 00:49:08,770 --> 00:49:14,140 If you tell me x dot, then I should solve my governing 1025 00:49:14,140 --> 00:49:17,710 equation for x naught. 1026 00:49:17,710 --> 00:49:22,150 You could tell me some equation that relates x naught and x 1027 00:49:22,150 --> 00:49:24,700 dot initially, and now I got to solve 1028 00:49:24,700 --> 00:49:27,760 this equation in conjunction with my governing equation 1029 00:49:27,760 --> 00:49:28,940 to determine these-- 1030 00:49:28,940 --> 00:49:31,210 this set of values. 1031 00:49:31,210 --> 00:49:34,600 I'd like to know both of these things initially. 1032 00:49:34,600 --> 00:49:38,770 The fully-implicit DAE, I really have two n unknowns here. 1033 00:49:38,770 --> 00:49:41,380 I don't know x, and I don't know x dot. 1034 00:49:41,380 --> 00:49:44,530 And I've only got n equations for those. 1035 00:49:44,530 --> 00:49:47,440 So apparently there's n degrees of freedom to specify here. 1036 00:49:47,440 --> 00:49:49,180 I need n more equations to say what 1037 00:49:49,180 --> 00:49:50,950 those initial conditions are. 1038 00:49:50,950 --> 00:49:53,110 And these hidden constraints I point out actually 1039 00:49:53,110 --> 00:49:55,690 reduce the number of degrees of freedom. 1040 00:49:55,690 --> 00:49:59,650 The fact that these are not index 0 problems 1041 00:49:59,650 --> 00:50:02,710 but have a finite index will introduce extra constraints 1042 00:50:02,710 --> 00:50:05,050 that reduce the number of degrees of freedom 1043 00:50:05,050 --> 00:50:08,734 or the number of other equations I can use to specify 1044 00:50:08,734 --> 00:50:09,775 those initial conditions. 1045 00:50:09,775 --> 00:50:11,020 Does that make sense? 1046 00:50:11,020 --> 00:50:14,230 Let me give you some examples, and then I'll let you go. 1047 00:50:14,230 --> 00:50:16,930 So if I have a problem that has separate differential 1048 00:50:16,930 --> 00:50:18,400 states and algebraic states, this 1049 00:50:18,400 --> 00:50:23,291 is like f of x dot x and y and t equal 0, the set of things 1050 00:50:23,291 --> 00:50:24,290 that I need to specify-- 1051 00:50:24,290 --> 00:50:27,787 I need the initial derivative of the differential state, 1052 00:50:27,787 --> 00:50:29,620 the initial value of the differential state, 1053 00:50:29,620 --> 00:50:32,130 and the initial value of the algebraic state. 1054 00:50:32,130 --> 00:50:34,930 I'd actually need to know the initial value of the derivative 1055 00:50:34,930 --> 00:50:36,820 of the algebraic state here because I 1056 00:50:36,820 --> 00:50:39,220 don't need that to satisfy this initial equation. 1057 00:50:41,950 --> 00:50:44,310 Let's look at a couple of examples. 1058 00:50:44,310 --> 00:50:47,400 So here's stirred tank example 1. 1059 00:50:47,400 --> 00:50:50,050 I converted to a system of ODEs by taking the derivative 1060 00:50:50,050 --> 00:50:51,680 of equation 2 here. 1061 00:50:51,680 --> 00:50:53,740 So here's my two ODEs. 1062 00:50:53,740 --> 00:50:55,682 I got to look back at my initial equations 1063 00:50:55,682 --> 00:50:57,890 and see how do they constrain the initial conditions. 1064 00:50:57,890 --> 00:51:01,080 So equation 2 tells me that c1 initially better 1065 00:51:01,080 --> 00:51:02,675 be gamma initially. 1066 00:51:02,675 --> 00:51:04,300 Otherwise, my initial condition doesn't 1067 00:51:04,300 --> 00:51:08,080 satisfy the governing equations of the model that I wrote down. 1068 00:51:08,080 --> 00:51:11,260 Equation 1 tells me that the derivative of c2 1069 00:51:11,260 --> 00:51:14,740 initially better be related to the derivative of c1 1070 00:51:14,740 --> 00:51:15,910 and the derivative of n-- 1071 00:51:15,910 --> 00:51:20,650 I'm sorry to the value of c1 and the value of c2 initially. 1072 00:51:20,650 --> 00:51:26,350 This third equation that I came up with, it doesn't constrain. 1073 00:51:26,350 --> 00:51:29,470 It's a derivative of one of the original algebraic variables, 1074 00:51:29,470 --> 00:51:31,877 so it doesn't constrain that original set of-- 1075 00:51:31,877 --> 00:51:34,210 I need to know the initial condition on the differential 1076 00:51:34,210 --> 00:51:36,824 variables, the derivative of the differential variables, 1077 00:51:36,824 --> 00:51:37,990 and the algebraic variables. 1078 00:51:37,990 --> 00:51:40,094 This only tells me something about the derivative 1079 00:51:40,094 --> 00:51:41,260 of this algebraic variables. 1080 00:51:41,260 --> 00:51:43,190 This equation doesn't give me any other constraints. 1081 00:51:43,190 --> 00:51:45,310 So I just got a free variable that I can specify. 1082 00:51:45,310 --> 00:51:46,160 One more equation. 1083 00:51:46,160 --> 00:51:47,690 It can be whatever I want. 1084 00:51:47,690 --> 00:51:50,500 So maybe I say c2 of 0c0, or maybe I 1085 00:51:50,500 --> 00:51:54,160 say there's some relationship between c1 of 0 and c2 0. 1086 00:51:54,160 --> 00:51:57,520 Or maybe I say that c2 of 0 and c2 prime of 0 1087 00:51:57,520 --> 00:51:59,170 are related in some way. 1088 00:51:59,170 --> 00:52:01,690 Whatever I specify, I have a system 1089 00:52:01,690 --> 00:52:05,837 of three equations for three unknowns, c1 of 0, c2 dot of 0, 1090 00:52:05,837 --> 00:52:07,270 and c2 of 0. 1091 00:52:07,270 --> 00:52:09,220 I need to have a unique solution to that. 1092 00:52:09,220 --> 00:52:11,890 Otherwise, there's going to be a problem solving for it in time, 1093 00:52:11,890 --> 00:52:13,010 the system of equations. 1094 00:52:13,010 --> 00:52:15,343 So I have to choose this one consistent with these other 1095 00:52:15,343 --> 00:52:18,220 ones, but beyond that I'm OK. 1096 00:52:22,549 --> 00:52:26,320 Here's stirred tank example 2. 1097 00:52:26,320 --> 00:52:29,230 So I take a derivative of equation 2, 1098 00:52:29,230 --> 00:52:32,053 and that gives me a differential equation for dc2 dt. 1099 00:52:32,053 --> 00:52:33,790 Actually, I already had one of those, 1100 00:52:33,790 --> 00:52:36,460 so I substitute for dc2 dt, and I 1101 00:52:36,460 --> 00:52:41,130 take another derivative to get a differential equation for c1. 1102 00:52:41,130 --> 00:52:43,070 And let's ask about the initial conditions. 1103 00:52:43,070 --> 00:52:47,730 So the initial algebraic equation, it constrained dc2-- 1104 00:52:47,730 --> 00:52:48,750 it constrained c2. 1105 00:52:48,750 --> 00:52:52,360 It said c2 of 0 has to be gamma 0. 1106 00:52:52,360 --> 00:52:57,130 The initial differential equation told me that c1 of 0 1107 00:52:57,130 --> 00:53:01,670 had to be equal to c2 of 0 plus v over q, c2 dot of 0. 1108 00:53:01,670 --> 00:53:05,260 Remember, I need three equations to specify c2 of 0, c1 of 0, 1109 00:53:05,260 --> 00:53:07,520 c2 dot of 0. 1110 00:53:07,520 --> 00:53:09,760 I actually can't specify c2 dot of 0 1111 00:53:09,760 --> 00:53:14,200 freely because in the process of deriving this differential 1112 00:53:14,200 --> 00:53:17,080 equation here so that I had dc1 dt and dc2 dt, 1113 00:53:17,080 --> 00:53:23,590 I introduced a constraint on the derivatives of c2. 1114 00:53:23,590 --> 00:53:28,150 And this tells me that c2 dot of 0 has to be gamma dot of 0. 1115 00:53:28,150 --> 00:53:31,380 There are no free variables to specify. 1116 00:53:31,380 --> 00:53:33,050 There's three equations. 1117 00:53:33,050 --> 00:53:35,230 Two of them come from the initial system of DAEs, 1118 00:53:35,230 --> 00:53:40,240 and the third one is related to this extra equation that 1119 00:53:40,240 --> 00:53:42,850 popped up as we tried to define the index. 1120 00:53:42,850 --> 00:53:46,370 It's an implicit constraint on the problem. 1121 00:53:46,370 --> 00:53:48,060 Does that make sense? 1122 00:53:48,060 --> 00:53:50,390 I can't specify these things freely. 1123 00:53:50,390 --> 00:53:55,690 They're specified by the signal that I measured gamma. 1124 00:53:55,690 --> 00:53:58,070 There are two examples here that you 1125 00:53:58,070 --> 00:54:01,250 should be able to work through. 1126 00:54:01,250 --> 00:54:03,695 One is an index 2 example, and one is an index 1 example. 1127 00:54:03,695 --> 00:54:05,570 We're not going to have time to go over them. 1128 00:54:05,570 --> 00:54:06,906 Sorry for that. 1129 00:54:06,906 --> 00:54:08,030 But they work the same way. 1130 00:54:08,030 --> 00:54:13,260 So convert to a system of ODEs, ask what initial conditions 1131 00:54:13,260 --> 00:54:15,690 will satisfy the initial constraints and were there 1132 00:54:15,690 --> 00:54:17,910 any other-- 1133 00:54:17,910 --> 00:54:19,770 were there any other constraining equations 1134 00:54:19,770 --> 00:54:20,784 that popped up? 1135 00:54:20,784 --> 00:54:22,950 So if you work through this example, you'll find out 1136 00:54:22,950 --> 00:54:27,370 there was a hidden strength that either the initial derivatives 1137 00:54:27,370 --> 00:54:29,535 of c1 and c2 sum together have to be 0 1138 00:54:29,535 --> 00:54:32,730 or equivalently the initial value of c3 1139 00:54:32,730 --> 00:54:33,840 had to be equal to 0. 1140 00:54:33,840 --> 00:54:35,400 You can take this constraint in either place. 1141 00:54:35,400 --> 00:54:36,690 Maybe you want to take both of them 1142 00:54:36,690 --> 00:54:38,400 to make sure that it's always satisfied, 1143 00:54:38,400 --> 00:54:40,969 that there is no numerical error that drives you off 1144 00:54:40,969 --> 00:54:41,760 of this constraint. 1145 00:54:41,760 --> 00:54:43,860 It's OK to have these things over constrained. 1146 00:54:43,860 --> 00:54:47,371 It's just not OK to have them under constrained. 1147 00:54:47,371 --> 00:54:48,870 And you'll find out that while there 1148 00:54:48,870 --> 00:54:52,140 were five things I could pick, the derivatives, and the values 1149 00:54:52,140 --> 00:54:54,810 of c1 and c2, and then the algebraic variable c3 1150 00:54:54,810 --> 00:54:55,330 initially. 1151 00:54:55,330 --> 00:54:57,510 So there are only four constraints-- for equations 1152 00:54:57,510 --> 00:54:58,890 here for the initial conditions. 1153 00:54:58,890 --> 00:54:59,848 I get to pick one more. 1154 00:54:59,848 --> 00:55:01,560 I can get it however I want. 1155 00:55:01,560 --> 00:55:04,410 Just has to be consistent with these. 1156 00:55:04,410 --> 00:55:05,740 There's one more example. 1157 00:55:05,740 --> 00:55:08,377 Again, should work through these on your own. 1158 00:55:08,377 --> 00:55:09,460 But you do the same thing. 1159 00:55:09,460 --> 00:55:12,650 Convert to a system of ODEs. 1160 00:55:12,650 --> 00:55:14,291 You'll find out that when you do this, 1161 00:55:14,291 --> 00:55:16,040 you'll get a differential equation for c3. 1162 00:55:16,040 --> 00:55:17,840 It didn't introduce any new constraints 1163 00:55:17,840 --> 00:55:21,094 on c1 dot, c2 dot, or c1, c2, and c3. 1164 00:55:21,094 --> 00:55:23,510 So you just need to satisfy-- your initial conditions have 1165 00:55:23,510 --> 00:55:25,310 to satisfy the initial equations, 1166 00:55:25,310 --> 00:55:28,130 and you can put in two more conditions, whatever 1167 00:55:28,130 --> 00:55:29,210 you want them to be. 1168 00:55:29,210 --> 00:55:31,130 Add five unknowns and only three equations 1169 00:55:31,130 --> 00:55:33,170 for the initial conditions. 1170 00:55:33,170 --> 00:55:34,500 And that's it. 1171 00:55:34,500 --> 00:55:36,810 So thank you very much.