WEBVTT 00:00:01.540 --> 00:00:03.910 The following content is provided under a Creative 00:00:03.910 --> 00:00:05.300 Commons license. 00:00:05.300 --> 00:00:07.510 Your support will help MIT OpenCourseWare 00:00:07.510 --> 00:00:11.600 continue to offer high quality educational resources for free. 00:00:11.600 --> 00:00:14.140 To make a donation or to view additional materials 00:00:14.140 --> 00:00:18.100 from hundreds of MIT courses, visit MIT OpenCourseWare 00:00:18.100 --> 00:00:19.310 at ocw.mit.edu. 00:00:24.488 --> 00:00:27.604 PROFESSOR: You get your quizzes back on Monday. 00:00:27.604 --> 00:00:29.520 Maybe you've had some time to think about them 00:00:29.520 --> 00:00:31.710 and reflect on your answers and the problems 00:00:31.710 --> 00:00:34.950 that were posed, maybe how they compared to the quiz that-- 00:00:34.950 --> 00:00:36.283 sample quiz that was put online. 00:00:36.283 --> 00:00:38.850 I think it was pretty fair exam overall. 00:00:38.850 --> 00:00:40.740 It looked very similar, actually, to the exam 00:00:40.740 --> 00:00:43.280 that we've given in the past, but we tried to-- 00:00:43.280 --> 00:00:46.860 tried not to put any sort of chemical engineering problems 00:00:46.860 --> 00:00:48.720 on top of the numerics, which tends 00:00:48.720 --> 00:00:51.160 to lead to confusion basically. 00:00:51.160 --> 00:00:52.830 Spent a lot of time reading and trying 00:00:52.830 --> 00:00:55.470 to understand what the underlying physical problem is 00:00:55.470 --> 00:01:00.450 instead of just testing the numerical methods. 00:01:00.450 --> 00:01:02.050 Are there any questions about the quiz 00:01:02.050 --> 00:01:02.940 now that you've gotten it back? 00:01:02.940 --> 00:01:04.920 You've seen your answers, things that were unclear that 00:01:04.920 --> 00:01:05.753 were asked on there? 00:01:05.753 --> 00:01:10.780 Now is a good time for reflection if there are. 00:01:10.780 --> 00:01:12.200 No. 00:01:12.200 --> 00:01:16.090 I heard this last homework assignment was really long, so 00:01:16.090 --> 00:01:16.800 sorry about that. 00:01:16.800 --> 00:01:17.800 I didn't write this one. 00:01:17.800 --> 00:01:19.285 [LAUGHS] 00:01:19.285 --> 00:01:21.730 I'm not going to Take credit for it. 00:01:21.730 --> 00:01:25.885 It's hard to judge sometimes how long these things take. 00:01:25.885 --> 00:01:30.280 A lot of it has to do with how facile you are with coding up 00:01:30.280 --> 00:01:31.192 some of the problems. 00:01:31.192 --> 00:01:32.650 So conceptually they can be simple, 00:01:32.650 --> 00:01:34.233 but the coding can be quite difficult, 00:01:34.233 --> 00:01:35.590 so it's hard to estimate. 00:01:35.590 --> 00:01:38.080 I think this week's, which should be posted 00:01:38.080 --> 00:01:41.020 and is on the topic of DAE, should be easier. 00:01:41.020 --> 00:01:43.380 It's only got two problems instead of three. 00:01:43.380 --> 00:01:45.130 You're likely going to be able to leverage 00:01:45.130 --> 00:01:46.780 some of the code from your solution 00:01:46.780 --> 00:01:49.160 last week for this week. 00:01:49.160 --> 00:01:52.210 So I don't think it should be quite so challenging. 00:01:52.210 --> 00:01:56.530 There is a last part to the first problem 00:01:56.530 --> 00:01:59.020 that seems to always give first-year graduate 00:01:59.020 --> 00:02:00.130 students trouble. 00:02:00.130 --> 00:02:02.470 Somehow you guys don't know how to write down 00:02:02.470 --> 00:02:05.050 energy balances on reactors. 00:02:05.050 --> 00:02:07.270 So you might ask pointed questions 00:02:07.270 --> 00:02:10.210 about this at your office hours, so you don't spend time trying 00:02:10.210 --> 00:02:14.200 to solve the wrong equations. 00:02:14.200 --> 00:02:14.780 This happens. 00:02:14.780 --> 00:02:19.060 I think many of you maybe were taught incorrectly or just 00:02:19.060 --> 00:02:22.690 forgot at some point how to write these things down. 00:02:22.690 --> 00:02:25.711 So the office hour is a good point to bring that up, 00:02:25.711 --> 00:02:27.460 and that way you don't spend a lot of time 00:02:27.460 --> 00:02:32.170 troubling yourself with the wrong system of equations. 00:02:32.170 --> 00:02:34.051 Any questions? 00:02:34.051 --> 00:02:34.550 Yes. 00:02:34.550 --> 00:02:39.542 AUDIENCE: [INAUDIBLE] 00:02:39.542 --> 00:02:40.500 PROFESSOR: What's that? 00:02:40.500 --> 00:02:44.980 I'd-- well, to be fair, almost all the lectures for the rest 00:02:44.980 --> 00:02:47.350 of the class are going to be taught by Professor Green. 00:02:47.350 --> 00:02:49.940 I'll step in and do some review for two sessions, 00:02:49.940 --> 00:02:53.710 review before the second quiz and review before the final. 00:02:53.710 --> 00:02:55.480 Since he's leading the lectures, he'll 00:02:55.480 --> 00:02:59.200 be responsible for largely generating the homework 00:02:59.200 --> 00:02:59.920 assignments. 00:02:59.920 --> 00:03:03.760 I'll try to help out with that as best I can, 00:03:03.760 --> 00:03:07.420 but this will be the last formal lecture 00:03:07.420 --> 00:03:10.210 you get from me this term. 00:03:10.210 --> 00:03:11.520 Other questions? 00:03:11.520 --> 00:03:13.720 Well, that's very, very sweet of you. 00:03:13.720 --> 00:03:17.040 Thank you so much. 00:03:17.040 --> 00:03:21.670 Least there wasn't any hushed yes! 00:03:21.670 --> 00:03:23.765 Today's lecture is our last one on DAE, 00:03:23.765 --> 00:03:26.512 so we're only going to do two. 00:03:26.512 --> 00:03:28.720 You're going to see today that differential algebraic 00:03:28.720 --> 00:03:31.960 equations are pretty complicated actually. 00:03:31.960 --> 00:03:34.870 And when they get sufficiently complicated, 00:03:34.870 --> 00:03:39.850 you really need to reach out to existing codes that 00:03:39.850 --> 00:03:42.790 are designed to solve particular classes of differential 00:03:42.790 --> 00:03:43.840 algebraic equations. 00:03:43.840 --> 00:03:46.270 So for you, it's more important to be 00:03:46.270 --> 00:03:49.600 able to identify in the models you formulate when 00:03:49.600 --> 00:03:52.420 these complications arise and what 00:03:52.420 --> 00:03:54.880 are the essential ingredients of the model that 00:03:54.880 --> 00:03:56.680 can be put into one of these solvers 00:03:56.680 --> 00:03:59.124 that you get a result that isn't nonsense. 00:03:59.124 --> 00:04:00.790 And that's what we're going to do today. 00:04:00.790 --> 00:04:02.510 There's going to be lots of examples to work through. 00:04:02.510 --> 00:04:04.180 So make sure you're sitting next to somebody 00:04:04.180 --> 00:04:05.800 you like because I'm going to ask you to try 00:04:05.800 --> 00:04:07.570 to think about these things and discuss things 00:04:07.570 --> 00:04:08.380 as we go through. 00:04:10.857 --> 00:04:12.940 Let me review where these complications come from. 00:04:12.940 --> 00:04:15.070 So last time we talked briefly about semi-explicit 00:04:15.070 --> 00:04:18.256 and fully-implicit differential algebraic equations. 00:04:18.256 --> 00:04:20.589 I told you in principle, you could simulate these things 00:04:20.589 --> 00:04:22.620 with backwards difference formulas 00:04:22.620 --> 00:04:25.327 in solving nonlinear equations at each time step 00:04:25.327 --> 00:04:26.410 and just marching forward. 00:04:26.410 --> 00:04:30.320 Did this in your homework for ODE IVP 00:04:30.320 --> 00:04:32.950 So you can do the same thing for differential algebraic 00:04:32.950 --> 00:04:33.670 equations. 00:04:33.670 --> 00:04:37.210 But there was a catch to that, and I illustrated that 00:04:37.210 --> 00:04:38.020 at the very end. 00:04:38.020 --> 00:04:40.980 So maybe your brain was fried at the end, and you missed this. 00:04:40.980 --> 00:04:42.740 It's good to recap. 00:04:42.740 --> 00:04:45.070 So I showed you an example of a stirred tank 00:04:45.070 --> 00:04:47.260 where you had some-- 00:04:47.260 --> 00:04:50.410 you had transport of some solute into the stirred tank, 00:04:50.410 --> 00:04:52.970 and then you're pulling it out at a different concentration. 00:04:52.970 --> 00:04:55.970 And we're trying to track the dynamics of the system, 00:04:55.970 --> 00:04:58.114 the concentration in and the concentration out, 00:04:58.114 --> 00:04:59.530 and we imagined a problem in which 00:04:59.530 --> 00:05:01.630 we measured and the concentration in 00:05:01.630 --> 00:05:04.390 and we tried to predict the concentration out. 00:05:04.390 --> 00:05:08.590 So we had a system of differential and algebraic 00:05:08.590 --> 00:05:10.540 equations to solve. 00:05:10.540 --> 00:05:13.600 And I showed you that if I used the backward Euler method, 00:05:13.600 --> 00:05:16.710 the lowest order backwards difference formula that I can-- 00:05:16.710 --> 00:05:18.730 of the canonical class of these things 00:05:18.730 --> 00:05:22.330 that I can craft, where I approximate the derivative 00:05:22.330 --> 00:05:27.250 of an unknown with a relative error or an error proportional 00:05:27.250 --> 00:05:31.502 to the time step delta t, that carrying out one time 00:05:31.502 --> 00:05:33.460 step with this backward differentiation formula 00:05:33.460 --> 00:05:35.990 would determine c1 in principle exactly-- 00:05:35.990 --> 00:05:38.765 if I knew gamma exactly, I would know c1 exactly. 00:05:38.765 --> 00:05:42.610 And it would determine c2 to within order delta t-squared. 00:05:42.610 --> 00:05:46.160 It's the local truncation error here was order delta t-squared. 00:05:46.160 --> 00:05:48.190 So you just substitute this formula in, 00:05:48.190 --> 00:05:50.860 and you ask does this error term change at all. 00:05:50.860 --> 00:05:52.900 At some point, I wind up multiplying by delta t, 00:05:52.900 --> 00:05:56.240 and so I go from order delta t to order delta t-squared It's 00:05:56.240 --> 00:05:58.240 the numerical error that gets carried around 00:05:58.240 --> 00:06:01.790 in this calculation. 00:06:01.790 --> 00:06:04.180 I just switch the model a little bit. 00:06:04.180 --> 00:06:07.140 It seems like a irrelevant change, 00:06:07.140 --> 00:06:09.210 but it turns out to be incredibly significant. 00:06:09.210 --> 00:06:11.520 So now I imagine a different problem 00:06:11.520 --> 00:06:15.490 in which I'm measuring c2 the outlet, 00:06:15.490 --> 00:06:19.740 and I want to predict c1 the inlet. 00:06:19.740 --> 00:06:22.290 I still have a system of DAEs that I have to solve, 00:06:22.290 --> 00:06:24.690 and if I applied the backwards Euler method, 00:06:24.690 --> 00:06:26.820 well, c2 is automatically determined 00:06:26.820 --> 00:06:28.350 by this algebraic equation. 00:06:28.350 --> 00:06:30.060 So I know that exactly. 00:06:30.060 --> 00:06:33.420 I got to go up to this first equation and solve it for c1, 00:06:33.420 --> 00:06:36.630 and when I do, I have to have an approximation 00:06:36.630 --> 00:06:37.500 for the derivative. 00:06:37.500 --> 00:06:40.576 And the derivative is proportional to delta t-- 00:06:40.576 --> 00:06:43.200 carries an error around with it that's proportional to delta t. 00:06:43.200 --> 00:06:47.879 So I know c1 to within order delta t, not order delta 00:06:47.879 --> 00:06:49.420 t-squared like in the previous model. 00:06:49.420 --> 00:06:51.300 So there's something fundamentally different 00:06:51.300 --> 00:06:53.130 about these two circumstances. 00:06:53.130 --> 00:06:57.550 And all I did was go from c1 to c2 in this algebraic equation. 00:06:57.550 --> 00:06:59.506 So that's peculiar, and that should 00:06:59.506 --> 00:07:01.380 make you really suspicious about your ability 00:07:01.380 --> 00:07:04.260 to solve these problems accurately. 00:07:04.260 --> 00:07:06.010 And here was the third example I gave you. 00:07:08.740 --> 00:07:12.010 So I said imagine this system of DAEs instead. 00:07:12.010 --> 00:07:13.750 So here's the differential equation. 00:07:13.750 --> 00:07:15.250 Here's an algebraic equation. 00:07:15.250 --> 00:07:17.340 Apply the backward Euler method. 00:07:17.340 --> 00:07:19.630 Well, c3 is determined automatically 00:07:19.630 --> 00:07:21.220 by this algebraic equation. 00:07:21.220 --> 00:07:22.680 I know it exactly. 00:07:22.680 --> 00:07:25.370 C2 is related to the derivative of c3, 00:07:25.370 --> 00:07:28.270 so I need to approximate it with my backward Euler derivative. 00:07:28.270 --> 00:07:31.222 And that picks up an error, order delta t. 00:07:31.222 --> 00:07:35.110 C1 is equal to the derivative of c2, 00:07:35.110 --> 00:07:37.270 so I need an approximation for c2. 00:07:37.270 --> 00:07:40.480 The derivative of c2 that's the backwards Euler approximation. 00:07:40.480 --> 00:07:44.070 That has an error that's proportional to order delta t. 00:07:44.070 --> 00:07:49.270 But c2 itself also has an error proportional to order delta t. 00:07:49.270 --> 00:07:53.110 And order delta t divided by delta t is order 1. 00:07:53.110 --> 00:07:54.250 So I get c3. 00:07:54.250 --> 00:07:55.660 I get c2. 00:07:55.660 --> 00:07:57.760 I even know what c1 is. 00:07:57.760 --> 00:08:00.280 I solve this problem with my backwards difference formula, 00:08:00.280 --> 00:08:06.230 and I have no resolution of c1, no clue what that value is. 00:08:06.230 --> 00:08:08.410 So this is some complicated control scheme 00:08:08.410 --> 00:08:11.260 that I've set up, and I need to know the value of c1 00:08:11.260 --> 00:08:13.210 to figure out how to operate this process. 00:08:13.210 --> 00:08:14.590 I'm lost. 00:08:14.590 --> 00:08:16.780 This is never going to work. 00:08:16.780 --> 00:08:19.330 So these three problems have fundamental differences 00:08:19.330 --> 00:08:20.560 between them. 00:08:20.560 --> 00:08:23.170 And I'm going to show you how to predict 00:08:23.170 --> 00:08:27.746 when these differences are going to occur, how to name them. 00:08:27.746 --> 00:08:30.180 So the stirred tank example 1, it 00:08:30.180 --> 00:08:33.110 carried a local truncation error, one time step order 00:08:33.110 --> 00:08:34.387 delta-t squared. 00:08:34.387 --> 00:08:36.720 Stirred tank example 2 carried a local truncation error. 00:08:36.720 --> 00:08:38.610 It's order delta-t. 00:08:38.610 --> 00:08:41.010 This third DAE example had a local truncation error 00:08:41.010 --> 00:08:41.730 that's order 1. 00:08:41.730 --> 00:08:45.030 No change in delta-t will improve my solution 00:08:45.030 --> 00:08:46.710 to that problem. 00:08:46.710 --> 00:08:50.570 It's independent of how I choose to do my time setting, which, 00:08:50.570 --> 00:08:54.757 of course, is ridiculous. 00:08:54.757 --> 00:08:55.840 So here's another problem. 00:08:59.490 --> 00:09:02.942 And I'd like you to try to do this. 00:09:02.942 --> 00:09:04.650 We'll see how far you can get through it. 00:09:04.650 --> 00:09:05.820 You don't have to get all the way through it, 00:09:05.820 --> 00:09:07.861 but see how far you can get through this problem. 00:09:07.861 --> 00:09:10.280 So can you apply the backward Euler method 00:09:10.280 --> 00:09:14.634 to this system of DAEs and try to predict how 00:09:14.634 --> 00:09:15.425 the air propagates. 00:09:15.425 --> 00:09:18.834 [INTERPOSING VOICES] 00:12:14.760 --> 00:12:17.410 OK, I don't want to break up the conversation too much, 00:12:17.410 --> 00:12:21.020 but tell me what you're finding as you try to do this. 00:12:21.020 --> 00:12:21.520 Yes. 00:12:21.520 --> 00:12:23.795 AUDIENCE: We found that c3 has to be 0. 00:12:23.795 --> 00:12:25.510 PROFESSOR: OK. 00:12:25.510 --> 00:12:27.610 How did you find that c3 has to be 0? 00:12:27.610 --> 00:12:29.466 AUDIENCE: Stick in this [INAUDIBLE] 00:12:29.466 --> 00:12:31.466 that's what the difference is of the [INAUDIBLE] 00:12:31.466 --> 00:12:32.250 PROFESSOR: Yes. 00:12:32.250 --> 00:12:34.400 AUDIENCE: [INAUDIBLE] 00:12:34.400 --> 00:12:35.130 PROFESSOR: Right. 00:12:35.130 --> 00:12:35.910 AUDIENCE: [INAUDIBLE] 00:12:35.910 --> 00:12:36.618 PROFESSOR: Right. 00:12:36.618 --> 00:12:38.536 AUDIENCE: [INAUDIBLE] 00:12:38.536 --> 00:12:40.655 PROFESSOR: OK, so this is very clever, 00:12:40.655 --> 00:12:45.500 and this has to do with how the model itself is formulated. 00:12:45.500 --> 00:12:49.040 So you did manipulations with the fundamental equations 00:12:49.040 --> 00:12:54.470 that I wrote down to determine that c3 had to be 0. 00:12:54.470 --> 00:12:55.465 That's good. 00:12:55.465 --> 00:12:56.090 We can do that. 00:12:56.090 --> 00:12:56.765 We can look at the equations. 00:12:56.765 --> 00:12:59.060 We can figure how to eliminate variables and find out-- 00:12:59.060 --> 00:13:01.604 c3 in this case always has to be 0. 00:13:01.604 --> 00:13:03.770 That may or may not be obvious, but one way to do it 00:13:03.770 --> 00:13:07.264 is take this constraint equation and take its derivative. 00:13:07.264 --> 00:13:08.930 This has to hold at every point in time, 00:13:08.930 --> 00:13:11.520 so its derivative must hold at every point in time. 00:13:11.520 --> 00:13:13.220 Substitute these other two equations 00:13:13.220 --> 00:13:18.080 in, c1 dot plus c2 dot eliminates c1 and c2, 00:13:18.080 --> 00:13:20.090 and all that will be left is c3 has to be 0. 00:13:20.090 --> 00:13:21.950 So if I manipulate these equations, 00:13:21.950 --> 00:13:24.369 I know c3 has to be 0. 00:13:24.369 --> 00:13:26.660 Let's suppose I don't manipulate the equations, though. 00:13:26.660 --> 00:13:30.080 Let's suppose I just put in the backward Euler approximation 00:13:30.080 --> 00:13:31.640 over here for the derivative, and I 00:13:31.640 --> 00:13:35.930 evaluate c1, c2 and c3 at the current time 00:13:35.930 --> 00:13:39.890 and ask what is c1, c2, and c3 at the current time in terms 00:13:39.890 --> 00:13:41.369 of what it was at the previous time 00:13:41.369 --> 00:13:43.160 And you get an answer that looks like this. 00:13:43.160 --> 00:13:44.930 So there's a simple system of equations 00:13:44.930 --> 00:13:45.960 that one has to solve. 00:13:45.960 --> 00:13:48.500 It's a three-by-three system, and it's easy to eliminate. 00:13:48.500 --> 00:13:49.910 I didn't really expect you-- 00:13:49.910 --> 00:13:51.230 actually expected this answer to come up 00:13:51.230 --> 00:13:52.771 because you guys are all very clever, 00:13:52.771 --> 00:13:54.770 so you try to make things easy beforehand. 00:13:54.770 --> 00:13:56.900 But it's-- there's a problem here. 00:13:56.900 --> 00:13:59.030 If I try to just do the backwards difference 00:13:59.030 --> 00:14:03.670 formula with the equations as they're written, 00:14:03.670 --> 00:14:05.420 I'll find out that when I solve this thing 00:14:05.420 --> 00:14:10.390 that c1 will be c1 minus c2 at the previous time step over 2 00:14:10.390 --> 00:14:12.990 plus an order delta t-squared error. 00:14:12.990 --> 00:14:15.950 C2 will take on the opposite sign plus an order delta 00:14:15.950 --> 00:14:18.590 t-squared error. 00:14:18.590 --> 00:14:25.280 C3 will be minus c1 plus c2 over 2 delta t 00:14:25.280 --> 00:14:28.580 plus an order delta t error. 00:14:28.580 --> 00:14:32.900 And actually if I were to apply successive approximations, 00:14:32.900 --> 00:14:34.820 different-- if I step-- 00:14:34.820 --> 00:14:37.310 take many steps in a row with this backward Euler method, 00:14:37.310 --> 00:14:39.601 there's nothing to guarantee that these constraints are 00:14:39.601 --> 00:14:42.260 satisfied exactly. 00:14:42.260 --> 00:14:44.330 I can have numerical error in my solution 00:14:44.330 --> 00:14:46.700 of this algebraic equation. 00:14:46.700 --> 00:14:49.340 So it isn't necessarily true that in my numerical solution 00:14:49.340 --> 00:14:52.400 c1 plus c2 is equal to 0. 00:14:52.400 --> 00:14:54.312 I can show you right here if I add c1 to c2. 00:14:54.312 --> 00:14:55.520 Well, these two are the same. 00:14:55.520 --> 00:14:58.370 They'll cancel, but am I guaranteed that these errors 00:14:58.370 --> 00:14:59.120 will cancel, too? 00:14:59.120 --> 00:15:00.786 Or will there be a small numerical error 00:15:00.786 --> 00:15:03.490 that propagates? 00:15:03.490 --> 00:15:06.170 So model formulation is key. 00:15:06.170 --> 00:15:09.530 This problem is like stirred tank problem 2. 00:15:09.530 --> 00:15:12.650 One of our solutions lost in order of accuracy 00:15:12.650 --> 00:15:14.279 in the local truncation error. 00:15:14.279 --> 00:15:15.320 It's not delta t-squared. 00:15:15.320 --> 00:15:18.160 It's order delta t. 00:15:18.160 --> 00:15:21.170 And what you want it to do, which is not necessarily 00:15:21.170 --> 00:15:23.840 a bad thing, was to try to use this constraint 00:15:23.840 --> 00:15:27.801 to somehow eliminate equations and simplify things. 00:15:27.801 --> 00:15:29.300 That can be helpful, but it can also 00:15:29.300 --> 00:15:31.460 be true that if I use this constraint to simplify 00:15:31.460 --> 00:15:34.580 equations and eliminate things, my numerical solutions 00:15:34.580 --> 00:15:36.860 may no longer satisfy the constraint at all. 00:15:36.860 --> 00:15:41.554 They may drift away from that algebraic equation. 00:15:41.554 --> 00:15:42.470 Does that makes sense? 00:15:42.470 --> 00:15:44.928 I'm not controlling the error with respect to that equation 00:15:44.928 --> 00:15:45.587 anymore. 00:15:50.950 --> 00:15:55.810 So how do you give a name to this kind of behavior? 00:15:55.810 --> 00:16:00.340 And the thing one talks about is the differential index 00:16:00.340 --> 00:16:02.017 of a DAE system. 00:16:02.017 --> 00:16:03.850 So let's look at the stirred tank example 1, 00:16:03.850 --> 00:16:06.550 and let's ask this question. 00:16:06.550 --> 00:16:09.400 How many time derivatives are needed 00:16:09.400 --> 00:16:12.910 to convert to a system of independent ODEs having 00:16:12.910 --> 00:16:14.770 differentials of all the unknowns. 00:16:14.770 --> 00:16:16.960 So there are two unknowns in this system. 00:16:16.960 --> 00:16:19.570 One is c1, and the other c2. 00:16:19.570 --> 00:16:20.980 Actually have one equation, which 00:16:20.980 --> 00:16:23.820 contains a differential of c2. 00:16:23.820 --> 00:16:25.960 How many derivatives of these equations with time 00:16:25.960 --> 00:16:31.660 do I need to take in order to get an equivalent ODE system? 00:16:31.660 --> 00:16:32.700 And it's just one. 00:16:32.700 --> 00:16:35.250 If I take the derivative of equation 2, 00:16:35.250 --> 00:16:37.140 I'll get a differential equation for c1. 00:16:37.140 --> 00:16:39.510 Now, I have a differential equation for c1 00:16:39.510 --> 00:16:42.540 and a differential equation for c2. 00:16:42.540 --> 00:16:44.280 DAEs of this type, where I only have 00:16:44.280 --> 00:16:47.790 to take the time derivative once of the algebraic constraints, 00:16:47.790 --> 00:16:51.664 are called index 1 DAEs. 00:16:51.664 --> 00:16:53.830 We saw that when we applied the backwards difference 00:16:53.830 --> 00:16:56.097 formula to stirred tank example 1, 00:16:56.097 --> 00:16:58.180 the local truncation error was the same as what we 00:16:58.180 --> 00:17:01.070 would get in an ODE IVP system. 00:17:01.070 --> 00:17:02.320 So index 1-- 00:17:02.320 --> 00:17:04.720 DAEs are easy to solve. 00:17:04.720 --> 00:17:08.619 They behave like ODE IVPs. 00:17:08.619 --> 00:17:10.626 You determine whether it's index 1 or not 00:17:10.626 --> 00:17:12.000 by asking how many derivatives do 00:17:12.000 --> 00:17:16.220 I have to take in order to get a system of independent ODEs. 00:17:16.220 --> 00:17:18.770 So I put this together with this, 00:17:18.770 --> 00:17:20.300 I have two independent ODEs. 00:17:20.300 --> 00:17:23.594 I could solve these, subject to some set of initial conditions, 00:17:23.594 --> 00:17:25.010 and the solution might be the same 00:17:25.010 --> 00:17:27.760 as the solution to the DAE. 00:17:27.760 --> 00:17:29.690 Let's do stirred tank example 2 now. 00:17:33.780 --> 00:17:34.710 So here we go. 00:17:34.710 --> 00:17:36.390 They look very similar, but now c2 00:17:36.390 --> 00:17:38.700 appears in the algebraic equation instead of c1. 00:17:38.700 --> 00:17:42.360 So let's take a derivative of the algebraic equation, 00:17:42.360 --> 00:17:44.960 and I get a differential equation for c2. 00:17:44.960 --> 00:17:47.250 But I already had a differential equation for c2. 00:17:47.250 --> 00:17:49.620 I want a differential equation for c1. 00:17:49.620 --> 00:17:50.660 So what do I do? 00:17:50.660 --> 00:17:53.380 While I know dc2 dt from up here. 00:17:53.380 --> 00:17:56.740 So let's substitute that in and solve for c1. 00:17:56.740 --> 00:17:58.260 So a substitute equation 1 in here. 00:17:58.260 --> 00:17:59.590 I solve for c1. 00:17:59.590 --> 00:18:01.320 Now, let's take another derivative. 00:18:01.320 --> 00:18:02.220 Call this equation 3. 00:18:02.220 --> 00:18:04.770 Derivative of equation 3 now gives me 00:18:04.770 --> 00:18:07.080 a differential equation for c1. 00:18:07.080 --> 00:18:09.770 It's also in terms of dc2 dt. 00:18:09.770 --> 00:18:10.560 If I want-- 00:18:10.560 --> 00:18:12.570 I don't have to-- but if I want to, 00:18:12.570 --> 00:18:14.370 I can substitute for that again to just get 00:18:14.370 --> 00:18:17.490 dc1 dt in terms of c1 and c2. 00:18:17.490 --> 00:18:21.390 But it took two derivatives to generate 00:18:21.390 --> 00:18:26.430 a system of independent ODEs from the DAEs. 00:18:26.430 --> 00:18:30.249 And this is called index 2. 00:18:30.249 --> 00:18:32.040 So it's a different character from index 1. 00:18:32.040 --> 00:18:33.070 We saw what happens. 00:18:33.070 --> 00:18:35.940 We tend to lose an order of accuracy at index 2. 00:18:35.940 --> 00:18:39.600 Somehow index 2 problems are more sensitive than index 1 00:18:39.600 --> 00:18:41.552 problems. 00:18:41.552 --> 00:18:42.510 Here's another example. 00:18:45.370 --> 00:18:49.610 So DAE sample 3 is going to proceed the same way. 00:18:49.610 --> 00:18:53.549 So first we need a differential equation. 00:18:53.549 --> 00:18:54.590 Which one's missing here? 00:18:54.590 --> 00:18:56.381 We don't have differential equation for c1, 00:18:56.381 --> 00:18:58.280 and that's what we're hunting for. 00:18:58.280 --> 00:19:00.590 So let's take a derivative of the third equation, 00:19:00.590 --> 00:19:01.640 the algebraic equation. 00:19:01.640 --> 00:19:03.710 We get c3 dot is gamma dot. 00:19:03.710 --> 00:19:05.630 And we already have an equation for c3 dot. 00:19:05.630 --> 00:19:07.000 That was 2. 00:19:07.000 --> 00:19:08.780 So substitute 2 in. 00:19:08.780 --> 00:19:12.039 So now we got c2 is gamma dot, but we're really 00:19:12.039 --> 00:19:13.580 looking for something in terms of c1, 00:19:13.580 --> 00:19:15.920 so let's take the derivative again. 00:19:15.920 --> 00:19:19.140 We've got c2 [? dot ?] is gamma dot dot. 00:19:19.140 --> 00:19:21.930 Doesn't help us much, but c2 dot we know is equal to c1. 00:19:21.930 --> 00:19:25.110 So we substitute one more time and take a third derivative. 00:19:25.110 --> 00:19:29.080 Now, we have a differential equation for c1. 00:19:29.080 --> 00:19:31.170 And we have differential equations for c2 and c3. 00:19:33.910 --> 00:19:36.340 We can take this is the differential equation for c3. 00:19:36.340 --> 00:19:38.060 This is the differential equation for c2. 00:19:38.060 --> 00:19:41.260 This is the differential equation for c1. 00:19:41.260 --> 00:19:44.500 Some subset of these can be chosen, 00:19:44.500 --> 00:19:46.960 and we can replace this algebraic equation 00:19:46.960 --> 00:19:49.280 with a differential equation. 00:19:49.280 --> 00:19:53.154 This is an index 3 DAE. 00:19:53.154 --> 00:19:55.070 It's not always a good idea to replace the DAE 00:19:55.070 --> 00:19:58.370 system with these derivatives. 00:19:58.370 --> 00:19:59.330 I can write down-- 00:19:59.330 --> 00:20:03.246 I can have a particular function that's equal to 0 00:20:03.246 --> 00:20:05.120 and higher-order derivatives of that function 00:20:05.120 --> 00:20:06.620 are also equal to 0. 00:20:06.620 --> 00:20:08.900 But the solution to those differential equations 00:20:08.900 --> 00:20:11.477 are not necessarily equal to that function. 00:20:11.477 --> 00:20:14.060 There are lots of functions that might have this same property 00:20:14.060 --> 00:20:16.435 that its derivative-- certain number of derivatives of it 00:20:16.435 --> 00:20:17.750 are equal to 0. 00:20:17.750 --> 00:20:19.610 So there have to be extra initial conditions 00:20:19.610 --> 00:20:21.620 or constraints on the solution that 00:20:21.620 --> 00:20:24.350 confine me to the manifold of solutions 00:20:24.350 --> 00:20:25.740 that reflects the DAEs. 00:20:25.740 --> 00:20:28.679 We'll talk about consistent initialization later on. 00:20:28.679 --> 00:20:30.220 And that's going to fix this problem. 00:20:30.220 --> 00:20:32.000 So it's not always a good idea to do this, 00:20:32.000 --> 00:20:33.440 but it's one way of understanding, 00:20:33.440 --> 00:20:38.570 given this model, what sort of sensitivity is it can exhibit. 00:20:38.570 --> 00:20:40.503 We try to calculate its differential index. 00:20:40.503 --> 00:20:41.003 Yes? 00:20:41.003 --> 00:20:47.282 AUDIENCE: So [INAUDIBLE] 00:20:47.282 --> 00:20:49.216 PROFESSOR: Yes. 00:20:49.216 --> 00:20:49.715 Yes. 00:20:49.715 --> 00:20:50.772 AUDIENCE: [INAUDIBLE] 00:20:50.772 --> 00:20:52.730 PROFESSOR: Well, we have to be careful. 00:20:52.730 --> 00:20:55.880 We need to choose a set that are independent of the others. 00:20:55.880 --> 00:20:56.750 OK. 00:20:56.750 --> 00:20:59.135 So like these three are clearly going to be independent, 00:20:59.135 --> 00:21:01.850 and it's going to be OK. 00:21:01.850 --> 00:21:04.687 It's going to be a problem if I choose this one and these two. 00:21:04.687 --> 00:21:06.770 They're not going to be independent of each other. 00:21:06.770 --> 00:21:08.895 So we have to select independent ones from the set. 00:21:08.895 --> 00:21:11.525 It may not be obvious what independent means in general. 00:21:11.525 --> 00:21:13.211 AUDIENCE: [INAUDIBLE] 00:21:13.211 --> 00:21:15.631 PROFESSOR: 1, 2, 6 could also work, yes. 00:21:15.631 --> 00:21:16.130 Yes. 00:21:21.500 --> 00:21:25.906 So in general-- not in general. 00:21:25.906 --> 00:21:27.530 Let's talk about the differential index 00:21:27.530 --> 00:21:29.810 of a semi-explicit DAE system. 00:21:29.810 --> 00:21:35.090 Remember semi-explicit meant that the differential variables 00:21:35.090 --> 00:21:37.940 can be written as x dot or dx dt is 00:21:37.940 --> 00:21:40.310 equal to some function of x and y. 00:21:40.310 --> 00:21:43.010 Y are the algebraic states, and they 00:21:43.010 --> 00:21:46.770 satisfy a separate equation, g of x and y and t is equal to 0. 00:21:46.770 --> 00:21:49.170 Semi-explicit form. 00:21:49.170 --> 00:21:52.720 So with a semi-explicit DAE, the differential index 00:21:52.720 --> 00:21:55.880 is defined as the minimum number of differentiations 00:21:55.880 --> 00:22:00.270 required to convert the DAE to a system of independent ODEs. 00:22:00.270 --> 00:22:01.610 What does that mean? 00:22:01.610 --> 00:22:03.890 Means let's take the algebraic equations 00:22:03.890 --> 00:22:06.170 and let's take a time derivative of them. 00:22:06.170 --> 00:22:07.970 That will give us a new function. 00:22:07.970 --> 00:22:10.370 Call it g1. 00:22:10.370 --> 00:22:14.090 It's going to be a function of the differential 00:22:14.090 --> 00:22:18.020 state, the algebraic state, and the time derivative 00:22:18.020 --> 00:22:20.116 of the algebraic state. 00:22:20.116 --> 00:22:21.740 It could also be a function of the time 00:22:21.740 --> 00:22:23.600 derivative of the differential state, 00:22:23.600 --> 00:22:27.310 but we know that in terms of f, which is a function of x and y. 00:22:27.310 --> 00:22:30.169 So let's not try putting x dot in there for convenience. 00:22:30.169 --> 00:22:31.460 Let's just write out like this. 00:22:31.460 --> 00:22:33.466 In principle, we can do this. 00:22:33.466 --> 00:22:34.840 Let's take two derivatives of it. 00:22:34.840 --> 00:22:35.839 It will be the same way. 00:22:35.839 --> 00:22:37.790 It'll give us some function g2, which 00:22:37.790 --> 00:22:39.950 is a function of the differential 00:22:39.950 --> 00:22:42.410 state, the algebraic state, and the time derivative 00:22:42.410 --> 00:22:44.360 of the algebraic state. 00:22:44.360 --> 00:22:46.660 Let's take as many of these as we need. 00:22:46.660 --> 00:22:48.500 It will give us a system of equations 00:22:48.500 --> 00:22:51.860 that we can eventually solve for the time derivatives of all 00:22:51.860 --> 00:22:55.939 the algebraic states and convert this DAE to a system of ODEs. 00:22:55.939 --> 00:22:57.980 And the question is how many of these derivatives 00:22:57.980 --> 00:23:00.140 do I have to take? 00:23:00.140 --> 00:23:01.700 Index 1, I need to take one. 00:23:01.700 --> 00:23:04.730 Index 2, I need both sets of equations 00:23:04.730 --> 00:23:08.780 in order to get something that I can actually solve for dy dt. 00:23:08.780 --> 00:23:11.450 Index 3, I get to take three derivatives and higher. 00:23:15.890 --> 00:23:18.110 OK, so here's an example. 00:23:18.110 --> 00:23:20.090 Let's see if we can work through this. 00:23:20.090 --> 00:23:23.930 So we have a DAE system. 00:23:23.930 --> 00:23:28.620 Can you calculate the differential index 00:23:28.620 --> 00:23:29.310 of this system? 00:23:32.296 --> 00:23:32.796 Go ahead. 00:23:59.220 --> 00:24:01.546 Not 0. 00:24:01.546 --> 00:24:03.290 Not 0. 00:24:03.290 --> 00:24:06.430 ODE IVP is-- are DAEs of index 0. 00:24:06.430 --> 00:24:08.190 They require no derivatives to generate 00:24:08.190 --> 00:24:10.710 a system of ordinary differential equations. 00:24:14.651 --> 00:24:15.400 What do you think? 00:24:15.400 --> 00:24:17.160 AUDIENCE: [INAUDIBLE] 00:24:17.160 --> 00:24:19.010 PROFESSOR: Index 2 you say. 00:24:19.010 --> 00:24:20.270 How do you come to index 2? 00:24:20.270 --> 00:24:23.050 Why is it index 2? 00:24:23.050 --> 00:24:23.720 How'd you do it? 00:24:23.720 --> 00:24:25.354 AUDIENCE: [INAUDIBLE] 00:24:25.354 --> 00:24:27.520 PROFESSOR: OK, differentiate the algebraic equation. 00:24:27.520 --> 00:24:29.215 AUDIENCE: Stick it in [INAUDIBLE] 4. 00:24:29.215 --> 00:24:30.640 You get c3 for 0. 00:24:30.640 --> 00:24:33.015 Then you have to differentiate [INAUDIBLE] 00:24:33.015 --> 00:24:36.770 PROFESSOR: Good so differentiate this equation. 00:24:36.770 --> 00:24:39.960 Substitute for c1 dot and c2 dot. 00:24:39.960 --> 00:24:44.100 C1 minus c1 will be 0 c2 minus c2 will be 0. 00:24:44.100 --> 00:24:46.150 So you have 2c3 equals 0. 00:24:46.150 --> 00:24:47.690 That's still an algebraic equation. 00:24:47.690 --> 00:24:51.000 We need one more derivative to get a differential equation 00:24:51.000 --> 00:24:54.510 for the only algebraic state we have, c3. 00:24:54.510 --> 00:24:59.790 So one more derivative will tell us dc3 dt is 0. 00:24:59.790 --> 00:25:02.700 So two derivatives to get differential equations. 00:25:02.700 --> 00:25:04.221 It's an 2 to DAE. 00:25:04.221 --> 00:25:04.720 Sam? 00:25:04.720 --> 00:25:07.090 SAM: What if you can get c3 equals 0 00:25:07.090 --> 00:25:09.223 and just eliminate it then wouldn't you 00:25:09.223 --> 00:25:12.278 write an equivalent some simple equations, but that would be 00:25:12.278 --> 00:25:12.778 [INAUDIBLE]? 00:25:12.778 --> 00:25:13.867 PROFESSOR: You could. 00:25:13.867 --> 00:25:15.450 So we could write an equivalent system 00:25:15.450 --> 00:25:18.500 of equations that says instead of this equation, 00:25:18.500 --> 00:25:20.490 report c3 equals 0. 00:25:20.490 --> 00:25:21.780 That's model formulation. 00:25:21.780 --> 00:25:23.430 Here we've given-- we're given a model, 00:25:23.430 --> 00:25:26.160 and we're asked what index it is. 00:25:26.160 --> 00:25:28.340 We could formulate a different model, 00:25:28.340 --> 00:25:30.090 and the model will have a different index. 00:25:30.090 --> 00:25:33.600 If I were to replace this equation with c3 equals 0, 00:25:33.600 --> 00:25:36.750 what would the index of the DAE system be instead? 00:25:36.750 --> 00:25:39.270 Index 1 instead. 00:25:39.270 --> 00:25:41.760 I already told you index 2 is harder to solve than index 1, 00:25:41.760 --> 00:25:43.680 so if you can formulate an index 1 DAE, 00:25:43.680 --> 00:25:45.240 you should probably do it. 00:25:45.240 --> 00:25:47.767 But it may not be obvious whether you have or haven't. 00:25:47.767 --> 00:25:49.350 So it's an issue of model formulation. 00:25:49.350 --> 00:25:51.760 It's a great question. 00:25:51.760 --> 00:25:55.650 Does that distinction-- is that clear? 00:25:55.650 --> 00:25:57.370 Or is it a little-- 00:25:57.370 --> 00:26:01.720 am I not making it clear how this works? 00:26:01.720 --> 00:26:02.950 Unresponsive. 00:26:02.950 --> 00:26:04.730 OK. 00:26:04.730 --> 00:26:06.630 It's OK. 00:26:06.630 --> 00:26:09.360 Let's-- there's a generic index 1 example that we can talk 00:26:09.360 --> 00:26:10.170 about actually. 00:26:10.170 --> 00:26:15.360 So let's take a look at a semi-explicit DAE system. 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. 00:26:23.240 --> 00:26:28.770 Let's take the time derivative of the algebraic equation. 00:26:28.770 --> 00:26:32.010 So the total derivative of g with respect to t 00:26:32.010 --> 00:26:33.990 is the Jacobian of g with respect 00:26:33.990 --> 00:26:38.500 to x times dx dt plus the Jacobian of g 00:26:38.500 --> 00:26:43.260 with respect to y times dy dt plus the derivative of g-- 00:26:43.260 --> 00:26:46.510 partial derivative of g with respect to t. 00:26:46.510 --> 00:26:47.890 And now let's solve. 00:26:47.890 --> 00:26:52.030 Let's push the dy dt term to the other side of the equation, 00:26:52.030 --> 00:26:53.830 and let's substitute for dx dt. 00:26:53.830 --> 00:26:57.920 We know dx dt is f, so you get this equation here. 00:26:57.920 --> 00:27:01.520 If the Jacobian is full rank-- 00:27:01.520 --> 00:27:05.840 if dg dy is full rank, then I can invert this matrix, 00:27:05.840 --> 00:27:08.420 and I automatically get my system of ODEs 00:27:08.420 --> 00:27:09.905 for the algebraic states. 00:27:16.150 --> 00:27:18.970 What's problematic about this equation 00:27:18.970 --> 00:27:21.289 when dg dy is not full rank? 00:27:21.289 --> 00:27:22.330 That should be a partial. 00:27:22.330 --> 00:27:23.570 That's sloppy. 00:27:23.570 --> 00:27:25.140 I was doing this at 1:00 last night, 00:27:25.140 --> 00:27:27.040 so that's why that's there. 00:27:27.040 --> 00:27:31.600 What's the-- what's problematic about this root finding problem 00:27:31.600 --> 00:27:35.850 here if dg dy or the Jacobian of g with respect to y 00:27:35.850 --> 00:27:36.880 is not full rank? 00:27:39.760 --> 00:27:40.470 You recall? 00:27:45.401 --> 00:27:45.900 If-- 00:27:45.900 --> 00:27:48.597 AUDIENCE: [INAUDIBLE] 00:27:48.597 --> 00:27:50.930 PROFESSOR: Yes, you won't be able to compute the inverse 00:27:50.930 --> 00:27:52.070 to get your dy dt's. 00:27:52.070 --> 00:27:52.962 That's true. 00:27:52.962 --> 00:27:54.920 Remember this Jacobian-- if it's not full rank, 00:27:54.920 --> 00:27:57.950 it's determinant is 0. 00:27:57.950 --> 00:28:00.200 There is no inverse of this thing. 00:28:00.200 --> 00:28:03.350 We talked once about the systems of nonlinear equations 00:28:03.350 --> 00:28:08.120 and locally unique solutions, the so-called inverse function 00:28:08.120 --> 00:28:09.130 theorem. 00:28:09.130 --> 00:28:13.010 We can only guarantee there is a locally unique solution when 00:28:13.010 --> 00:28:15.040 I can invert this thing. 00:28:15.040 --> 00:28:17.000 That if I got in really close to the solution, 00:28:17.000 --> 00:28:18.890 it looked like a linear system of equations. 00:28:18.890 --> 00:28:21.820 That linear system of equations had a unique solution 00:28:21.820 --> 00:28:22.570 associated within. 00:28:22.570 --> 00:28:23.570 And everything was good. 00:28:23.570 --> 00:28:25.640 We were happy. 00:28:25.640 --> 00:28:28.250 So index 1 DAEs have that property. 00:28:28.250 --> 00:28:30.700 If I wanted to solve this equation for y, 00:28:30.700 --> 00:28:34.739 there will be a locally unique solution for y. 00:28:34.739 --> 00:28:36.530 Higher index DAEs don't have that property. 00:28:36.530 --> 00:28:39.970 We know that's a sort of unhappy generic circumstance 00:28:39.970 --> 00:28:42.370 to be in if we have to solve this equation. 00:28:42.370 --> 00:28:44.710 It can be hard to find those roots using, 00:28:44.710 --> 00:28:47.470 say, Newton-Raphson because the Jacobian [? e ?] 00:28:47.470 --> 00:28:49.870 would need to become singular during the root finding 00:28:49.870 --> 00:28:52.570 procedure. 00:28:52.570 --> 00:28:54.580 This is connected intimately to what 00:28:54.580 --> 00:28:57.280 we did for systems of nonlinear equations. 00:28:57.280 --> 00:28:59.410 So for an index 1 DAE, you can show its index 00:28:59.410 --> 00:29:04.880 1 because this Jacobian is invertible. 00:29:04.880 --> 00:29:07.480 Let's do some more examples on differential index. 00:29:07.480 --> 00:29:09.230 This seems to be the most important thing. 00:29:09.230 --> 00:29:10.790 If I have a model, what's its index 00:29:10.790 --> 00:29:13.040 because its index is going to tell me how sensitive it 00:29:13.040 --> 00:29:16.840 is to small perturbations. 00:29:16.840 --> 00:29:17.991 So here's a model. 00:29:17.991 --> 00:29:19.740 It's the same as the model you saw before, 00:29:19.740 --> 00:29:21.240 but now it's two mixing tanks. 00:29:21.240 --> 00:29:26.770 So here comes an inlet flow q1 carrying concentration c1. 00:29:26.770 --> 00:29:30.540 Out comes the same flow q1 carrying concentration c2. 00:29:30.540 --> 00:29:33.540 The mixer has volume v1. 00:29:33.540 --> 00:29:36.630 And then I blend this with some more 00:29:36.630 --> 00:29:40.640 of whatever the solute is at a flow rate of q2 00:29:40.640 --> 00:29:43.140 and a concentration c3. 00:29:43.140 --> 00:29:44.910 And those both go into this tank, 00:29:44.910 --> 00:29:48.910 and they come out concentration c4, flow rate q1 plus q2, 00:29:48.910 --> 00:29:52.076 and this has volume v2. 00:29:52.076 --> 00:29:53.320 Does that look good? 00:29:53.320 --> 00:29:54.760 Well-posed model? 00:29:54.760 --> 00:29:58.525 OK, so here's your material balances on the mixers. 00:29:58.525 --> 00:30:00.400 And I'm going to give you different algebraic 00:30:00.400 --> 00:30:01.450 constraints. 00:30:01.450 --> 00:30:06.120 So this is a problem now where we say measure c1 and c3. 00:30:09.400 --> 00:30:14.175 And we want to try to predict c2 and c4. 00:30:14.175 --> 00:30:16.120 Can you figure out the differential index 00:30:16.120 --> 00:30:17.160 of this DAE system? 00:30:26.640 --> 00:30:28.010 Feel free to talk to each other. 00:30:28.010 --> 00:30:28.510 It's OK. 00:30:59.210 --> 00:31:01.820 OK, the differential index is 1. 00:31:01.820 --> 00:31:03.410 This one's the easy one right? 00:31:03.410 --> 00:31:06.185 Just take one derivative of the algebraic equations. 00:31:06.185 --> 00:31:11.000 Now, I've got dc2 dt, dc4 dt, dc1 dt, dc3 dt. 00:31:11.000 --> 00:31:16.560 Differential index of 1 is easy to solve as an ODE IVP. 00:31:16.560 --> 00:31:18.800 This is the natural problem, too. 00:31:18.800 --> 00:31:22.640 It's useful to think about the inputs c1, c3 and asking 00:31:22.640 --> 00:31:23.930 about what the outputs are. 00:31:23.930 --> 00:31:26.765 Physically, this problem seems easier in nature. 00:31:32.570 --> 00:31:35.280 OK let's change it now. 00:31:35.280 --> 00:31:36.649 Same problem. 00:31:36.649 --> 00:31:38.440 But let's change the algebraic constraints. 00:31:38.440 --> 00:31:42.040 Now, the algebraic constraint is I measure c3 and c4, 00:31:42.040 --> 00:31:43.815 And I want to predict c1 and c2. 00:31:46.970 --> 00:31:48.450 What's the differential index here? 00:32:53.600 --> 00:32:54.100 Yes. 00:32:54.100 --> 00:32:54.980 AUDIENCE: [INAUDIBLE] 00:32:54.980 --> 00:32:56.250 PROFESSOR: Sure. 00:32:56.250 --> 00:32:58.600 AUDIENCE: Why was this a differential index of 1 00:32:58.600 --> 00:33:01.890 if you have to take a derivative for c1 and this equation? 00:33:01.890 --> 00:33:03.800 PROFESSOR: Oh, good question. 00:33:03.800 --> 00:33:06.810 I'm going to push the slide back, and I apologize for this. 00:33:06.810 --> 00:33:07.979 And I'll go back for it. 00:33:07.979 --> 00:33:10.270 So the question was why was this differential index one 00:33:10.270 --> 00:33:13.100 if I had to take a derivative of both this equation 00:33:13.100 --> 00:33:14.770 and that equation. 00:33:14.770 --> 00:33:17.260 So the way to think about this conceptually 00:33:17.260 --> 00:33:21.789 is I took a derivative of the algebraic equations, 00:33:21.789 --> 00:33:23.080 the set of algebraic equations. 00:33:23.080 --> 00:33:29.366 It was one time derivative of a vector valued function. 00:33:29.366 --> 00:33:30.990 I can see where the ambiguity is there. 00:33:30.990 --> 00:33:32.240 So it's important to be clear. 00:33:32.240 --> 00:33:35.230 I took one time derivative over the equations that 00:33:35.230 --> 00:33:37.400 prescribe the algebraic states. 00:33:41.260 --> 00:33:42.280 That's the distinction. 00:33:42.280 --> 00:33:43.780 So this is index 1 because I needed 00:33:43.780 --> 00:33:55.600 one time derivative of this type, dg dt of the entire set. 00:33:55.600 --> 00:33:56.100 Here we go. 00:33:56.100 --> 00:33:57.550 We're on this one. 00:33:57.550 --> 00:33:58.340 Prescribe c3. 00:33:58.340 --> 00:33:59.270 Prescribe c4. 00:33:59.270 --> 00:34:00.860 What is the differential index now? 00:34:06.920 --> 00:34:07.712 What do you think? 00:34:07.712 --> 00:34:08.420 We got an answer? 00:34:08.420 --> 00:34:10.670 [INTERPOSING VOICES] 00:34:10.670 --> 00:34:11.449 I hear 2. 00:34:11.449 --> 00:34:13.412 I hear 3. 00:34:13.412 --> 00:34:14.860 I saw 3 like this. 00:34:14.860 --> 00:34:16.719 Almost mistook it for a 2, but that's 3. 00:34:16.719 --> 00:34:17.590 Yes. 00:34:17.590 --> 00:34:18.337 It's 3. 00:34:18.337 --> 00:34:20.170 Let's see if we can work through why it's 3. 00:34:23.330 --> 00:34:26.300 So I take a derivative of the algebraic equations, 00:34:26.300 --> 00:34:27.260 one derivative. 00:34:27.260 --> 00:34:28.834 And after taking that derivative, 00:34:28.834 --> 00:34:30.500 I'll get a differential equation for one 00:34:30.500 --> 00:34:31.730 of the algebraic states. 00:34:31.730 --> 00:34:33.889 So c3 is taking care of. 00:34:33.889 --> 00:34:37.010 Now, I'm in the hunt for c1. 00:34:37.010 --> 00:34:38.780 So I have-- after that first derivative, 00:34:38.780 --> 00:34:42.320 I have an equation dc4 dt is gamma 2. 00:34:42.320 --> 00:34:45.530 And I know dc4 dt, so I drop that in, 00:34:45.530 --> 00:34:49.670 and I get an algebraic equation relating the derivative 00:34:49.670 --> 00:34:54.250 of gamma 2 to c2, c3, and c4. 00:34:54.250 --> 00:34:56.210 This new algebraic equation. 00:34:56.210 --> 00:34:59.120 This is like the g1 that I prescribed 00:34:59.120 --> 00:35:00.100 in the generic scene. 00:35:00.100 --> 00:35:02.725 I take a derivative of it again because I'm in the hunt for c1. 00:35:02.725 --> 00:35:04.760 So I take the derivative of this new equation, 00:35:04.760 --> 00:35:07.160 and I'll get a derivative of c2, a derivative of c3, 00:35:07.160 --> 00:35:08.150 and a derivative of c4. 00:35:08.150 --> 00:35:09.620 And I know all those. 00:35:09.620 --> 00:35:11.140 I know the derivative of c2. 00:35:11.140 --> 00:35:13.140 I know the derivative of c4. 00:35:13.140 --> 00:35:16.130 And I know the derivative of c3 from the previous level 00:35:16.130 --> 00:35:17.040 in the hierarchy. 00:35:17.040 --> 00:35:18.490 I figured that out already. 00:35:18.490 --> 00:35:19.810 So it's two derivatives. 00:35:19.810 --> 00:35:22.460 Still hasn't gotten me a time derivative of c1. 00:35:22.460 --> 00:35:23.960 But when I make those substitutions, 00:35:23.960 --> 00:35:27.390 I'll get an algebraic equation in terms of c1. 00:35:27.390 --> 00:35:29.780 And I can take one more derivative, 00:35:29.780 --> 00:35:31.740 and that will give me a dc1 dt. 00:35:31.740 --> 00:35:34.430 And I'll have an od for c1, c2, c3, and c4. 00:35:34.430 --> 00:35:36.900 So its differential index 3. 00:35:41.050 --> 00:35:42.340 It makes sense sort of. 00:35:42.340 --> 00:35:47.630 I'm taking a measurement of an output way down the line here. 00:35:47.630 --> 00:35:51.120 And I'm trying to predict what the input was 00:35:51.120 --> 00:35:51.966 in the first place. 00:35:51.966 --> 00:35:53.340 It's easy to imagine that there's 00:35:53.340 --> 00:35:56.220 a huge amount of sensitivity in that calculation. 00:36:00.570 --> 00:36:02.182 Here's another one. 00:36:02.182 --> 00:36:03.640 What's the differential index here? 00:36:10.310 --> 00:36:14.790 I'm now prescribing c1 and c2, and I 00:36:14.790 --> 00:36:16.060 want you to tell me c3 and c4. 00:36:16.060 --> 00:36:24.870 [INTERPOSING VOICES] 00:36:24.870 --> 00:36:25.670 Oh, good. 00:36:25.670 --> 00:36:26.870 Somebody noticed early on. 00:36:26.870 --> 00:36:29.300 So you say it's not possible. 00:36:29.300 --> 00:36:30.259 Right. 00:36:30.259 --> 00:36:31.175 Why isn't it possible? 00:36:31.175 --> 00:36:33.530 AUDIENCE: You take derivatives of both of those-- 00:36:33.530 --> 00:36:34.155 PROFESSOR: Yes. 00:36:34.155 --> 00:36:37.202 AUDIENCE: [INAUDIBLE] you can't isolate c3-- 00:36:37.202 --> 00:36:38.980 PROFESSOR: Right. 00:36:38.980 --> 00:36:42.040 We'll never be able to isolate a derivative of c3 here. 00:36:42.040 --> 00:36:44.230 There's actually there's something wrong physically 00:36:44.230 --> 00:36:45.070 with this problem. 00:36:47.830 --> 00:36:52.780 Yes, somehow I'm supposed to measure c1 and c2 00:36:52.780 --> 00:36:56.740 and use that to predict c3 and c4. 00:36:56.740 --> 00:36:59.592 I don't know c4, and c3 is an input. 00:36:59.592 --> 00:37:00.550 I don't know it either. 00:37:00.550 --> 00:37:05.260 How do I figure out an input when I don't know the output? 00:37:05.260 --> 00:37:06.430 It's impossible. 00:37:06.430 --> 00:37:10.554 So this model is flawed. 00:37:10.554 --> 00:37:12.220 We can formulate any model we want, pick 00:37:12.220 --> 00:37:15.920 any set of these variables to prescribe algebraically, 00:37:15.920 --> 00:37:20.050 but not all of them are going to admit solutions or make sense. 00:37:20.050 --> 00:37:21.820 There's no number of derivatives that's 00:37:21.820 --> 00:37:25.024 going to give us a closed system of ODEs. 00:37:25.024 --> 00:37:26.440 This again comes back to the point 00:37:26.440 --> 00:37:29.832 that really DAEs is all about model formulation. 00:37:29.832 --> 00:37:31.540 There are lots of good numerical methods. 00:37:31.540 --> 00:37:34.090 They work like the numerical methods for ODE IVPs. 00:37:34.090 --> 00:37:37.150 What's important is getting consistent models down, 00:37:37.150 --> 00:37:39.850 models that, for example, have solutions. 00:37:39.850 --> 00:37:41.650 You might feed this to a DAE solver 00:37:41.650 --> 00:37:45.460 and get nonsense because I can almost freely prescribe 00:37:45.460 --> 00:37:48.790 what c3 is, and I'll get some answer for c4. 00:37:48.790 --> 00:37:50.450 Any c3 can give me an answer here. 00:37:50.450 --> 00:37:53.740 So who knows whether this numerical solver is sensitive 00:37:53.740 --> 00:37:56.999 or not to this particular pathology in the problem 00:37:56.999 --> 00:37:57.790 we've written down. 00:38:02.080 --> 00:38:04.710 We got to speed up just a little bit, and I'm sorry for that. 00:38:04.710 --> 00:38:06.540 So we looked at index 1. 00:38:06.540 --> 00:38:07.590 That was stirred tank 1. 00:38:07.590 --> 00:38:09.450 Index 2, that was stirred tank 2. 00:38:09.450 --> 00:38:11.439 Index 3, that was stirred tank 3. 00:38:11.439 --> 00:38:13.230 And one thing we saw when we looked at them 00:38:13.230 --> 00:38:15.360 was this index 1 solution. 00:38:15.360 --> 00:38:17.970 It had pieces that were proportional 00:38:17.970 --> 00:38:22.470 to this forcing function gamma, this prescribed function 00:38:22.470 --> 00:38:26.710 in our system of equations and proportional to its derivative. 00:38:26.710 --> 00:38:29.510 So if gamma is jumping around, well, c1 will jump around. 00:38:29.510 --> 00:38:32.220 C2 is going to be smoothed out version of gamma 00:38:32.220 --> 00:38:34.630 because it depends on its integral. 00:38:34.630 --> 00:38:37.740 It's not very sensitive, kind of like an ODE IVP. 00:38:37.740 --> 00:38:40.260 Doesn't really show any more sensitivity than an ODE IVP 00:38:40.260 --> 00:38:41.250 does. 00:38:41.250 --> 00:38:43.350 The dynamics are different for index 2, though. 00:38:43.350 --> 00:38:46.050 When we look at the solution of the index 2 DAE, 00:38:46.050 --> 00:38:51.390 we found c2 goes like gamma, but c1 has to go like gamma dot. 00:38:51.390 --> 00:38:54.549 So that's pretty sensitive if you're making a measurement, 00:38:54.549 --> 00:38:56.090 and there's noise in the measurement. 00:38:56.090 --> 00:38:58.170 How do you even know what this derivative is? 00:38:58.170 --> 00:39:01.020 So c1 is wildly bouncing around. 00:39:01.020 --> 00:39:03.060 Our prediction of what c1 is is-- 00:39:03.060 --> 00:39:05.162 it's not going to be very good. 00:39:05.162 --> 00:39:09.330 Our solution of DAE example 3 when it told us c1 00:39:09.330 --> 00:39:11.730 was related to the second derivative 00:39:11.730 --> 00:39:15.332 of the forcing function, c2 to the first and c3 00:39:15.332 --> 00:39:17.040 was proportional to the forcing function. 00:39:17.040 --> 00:39:19.100 So this is hugely sensitive to changes 00:39:19.100 --> 00:39:20.100 in the forcing function. 00:39:20.100 --> 00:39:24.030 So the higher the index goes, the greater the sensitivity 00:39:24.030 --> 00:39:25.425 to perturbations in the system. 00:39:29.840 --> 00:39:31.540 Here's another simple example. 00:39:31.540 --> 00:39:33.073 You have all that data, actually, in your paper, 00:39:33.073 --> 00:39:35.614 so I won't ask you to do this one given our time constraints. 00:39:35.614 --> 00:39:39.040 But here's-- mechanical systems that have constraints are often 00:39:39.040 --> 00:39:41.260 indexed 3 DAEs it turns out. 00:39:41.260 --> 00:39:43.390 You can show this one is an index 3 DAE. 00:39:43.390 --> 00:39:47.170 This is the case of a pendulum swinging back and forth. 00:39:47.170 --> 00:39:50.740 So it's a-- change in position is its velocity. 00:39:50.740 --> 00:39:53.230 Its acceleration balances gravity, 00:39:53.230 --> 00:39:56.080 and there's some arm that holds the pendulum in place. 00:39:56.080 --> 00:39:58.420 We can imagine that that arm has some tension in it. 00:39:58.420 --> 00:40:01.510 It acts like a spring with a spring constant 00:40:01.510 --> 00:40:05.260 that changes in time in order to hold the pendulum at a fixed 00:40:05.260 --> 00:40:07.150 distance from its center. 00:40:07.150 --> 00:40:11.300 So two differential equations, one algebraic equation. 00:40:11.300 --> 00:40:13.660 This is a differential algebraic equation. 00:40:13.660 --> 00:40:16.270 You can imagine lots of mechanical systems 00:40:16.270 --> 00:40:18.130 work in exactly this way. 00:40:18.130 --> 00:40:20.490 They can only move along prescribed paths. 00:40:20.490 --> 00:40:23.679 They're constrained in how far they can stretch or go. 00:40:23.679 --> 00:40:24.220 I don't know. 00:40:24.220 --> 00:40:26.310 If you're trying to design a robot or something, 00:40:26.310 --> 00:40:29.180 DAEs are pretty important, and they're all of index 3 type 00:40:29.180 --> 00:40:30.370 it turns out. 00:40:30.370 --> 00:40:33.730 So differential variables here are x and v, 00:40:33.730 --> 00:40:36.700 and the algebraic variable is the spring constant, k, 00:40:36.700 --> 00:40:40.690 which has got to adjust dynamically in time in order 00:40:40.690 --> 00:40:41.460 to-- 00:40:41.460 --> 00:40:43.150 got to get away from the speaker-- 00:40:43.150 --> 00:40:44.710 adjust dynamically in time in order 00:40:44.710 --> 00:40:48.310 to provide just the right stiffness to keep this going 00:40:48.310 --> 00:40:51.170 around on a circular orbit. 00:40:51.170 --> 00:40:54.280 So when-- let's suppose I start with my pendulum down, 00:40:54.280 --> 00:40:56.170 and it's not moving, then this has 00:40:56.170 --> 00:40:59.530 to be just stiff enough to balance gravity. 00:40:59.530 --> 00:41:01.916 The pendulum is 90 degrees, and it's not moving. 00:41:01.916 --> 00:41:03.040 I don't need any stiffness. 00:41:03.040 --> 00:41:04.570 K and just be 0. 00:41:04.570 --> 00:41:06.640 There's no forces to balance. 00:41:06.640 --> 00:41:08.350 The pendulum is swinging around, and I 00:41:08.350 --> 00:41:10.570 have to be stiff enough to balance gravity 00:41:10.570 --> 00:41:14.500 where I am and also counteract any sort of centripetal 00:41:14.500 --> 00:41:16.180 [? acceleration. ?] So this k has 00:41:16.180 --> 00:41:18.970 got to be wildly fluctuating. 00:41:18.970 --> 00:41:21.280 If I was trying to control that k somehow 00:41:21.280 --> 00:41:25.000 to give this system these particular dynamics, 00:41:25.000 --> 00:41:29.270 you might imagine it's very difficult to do. 00:41:29.270 --> 00:41:31.340 So there's your differential variables. 00:41:31.340 --> 00:41:32.780 Here's your algebraic variable. 00:41:32.780 --> 00:41:34.437 You take a derivative of this equation, 00:41:34.437 --> 00:41:36.020 you'll get a constraint that tells you 00:41:36.020 --> 00:41:38.554 the velocity is orthogonal to the position. 00:41:38.554 --> 00:41:39.220 Of course it is. 00:41:39.220 --> 00:41:41.330 I'm on a circular trajectory. 00:41:41.330 --> 00:41:43.430 You take a derivative of this equation 00:41:43.430 --> 00:41:44.990 now, this algebraic equation. 00:41:44.990 --> 00:41:46.823 You'll get another constraint that gives you 00:41:46.823 --> 00:41:49.760 some relationship between k, x, and v 00:41:49.760 --> 00:41:51.860 but not a differential equation for k. 00:41:51.860 --> 00:41:54.456 Take another derivative of this algebraic equation, 00:41:54.456 --> 00:41:56.330 and you'll get a differential equation for k. 00:41:56.330 --> 00:41:58.660 So in principle, I could formulate a system of ODEs 00:41:58.660 --> 00:42:02.980 with dk dt, dx dt, dv dt will be equivalent to this. 00:42:02.980 --> 00:42:05.596 It took three derivatives to do this, though. 00:42:05.596 --> 00:42:07.220 So it's an incredibly sensitive system. 00:42:07.220 --> 00:42:08.460 It's index 3 in nature. 00:42:11.510 --> 00:42:13.030 So if I transform to this equivalent 00:42:13.030 --> 00:42:14.500 set of ODEs, the problem-- 00:42:14.500 --> 00:42:17.140 we discussed this earlier-- is that the solutions may 00:42:17.140 --> 00:42:21.197 drift away from the initial set of constraints. 00:42:21.197 --> 00:42:22.780 The solutions also need the right sort 00:42:22.780 --> 00:42:24.210 of initial conditions. 00:42:24.210 --> 00:42:31.400 Here I know that at time 0, c3 better be gamma 0. 00:42:31.400 --> 00:42:33.860 I know that in the actual solution to this problem, 00:42:33.860 --> 00:42:36.440 c3 dot was gamma dot. 00:42:36.440 --> 00:42:41.570 So at time 0, c2 better be gamma dot of 0. 00:42:41.570 --> 00:42:44.840 And at time 0, c1 better be gamma dot dot of 0. 00:42:44.840 --> 00:42:48.160 And if it's not, then there's going 00:42:48.160 --> 00:42:50.800 to be some step jump or some strange behavior 00:42:50.800 --> 00:42:53.210 in the solution to this ODE system. 00:42:53.210 --> 00:42:55.675 So I have to choose the right sort of initial conditions. 00:42:55.675 --> 00:42:57.850 If those initial conditions aren't chosen correctly, 00:42:57.850 --> 00:42:59.090 then I'm not going to get the right solution. 00:42:59.090 --> 00:43:01.990 I won't be constrained to the manifold of solutions that's 00:43:01.990 --> 00:43:03.655 given by the original DAE. 00:43:06.446 --> 00:43:08.070 So here's some things you need to know. 00:43:08.070 --> 00:43:12.680 Index 1, semi-explicit DAEs can be safely handled in MATLAB. 00:43:12.680 --> 00:43:16.670 So ode15s, ode23t, they have the ability 00:43:16.670 --> 00:43:18.870 to take an input of mass matrix. 00:43:18.870 --> 00:43:21.830 We discussed mass matrix last time. 00:43:21.830 --> 00:43:24.055 And a right hand side to the system of equations 00:43:24.055 --> 00:43:27.800 in f of x and t and solve it just 00:43:27.800 --> 00:43:30.380 like all the other solvers, all the other ODE IVP 00:43:30.380 --> 00:43:32.540 solvers in MATLAB. 00:43:32.540 --> 00:43:35.900 If it's not index 1, MATLAB can catch it. 00:43:35.900 --> 00:43:39.110 So it'll try to look at the Jacobian of f with respect 00:43:39.110 --> 00:43:40.970 to the variables and the mass matrix 00:43:40.970 --> 00:43:45.680 and to infer from that what the differential index is. 00:43:45.680 --> 00:43:48.759 And it will tell you often times-- 00:43:48.759 --> 00:43:50.550 I think depends which package you're using. 00:43:50.550 --> 00:43:52.910 But if you're using certain packages, it'll tell you. 00:43:52.910 --> 00:43:55.170 It's not index 1oe-- 00:43:55.170 --> 00:43:56.360 DAE. 00:43:56.360 --> 00:43:57.770 We can't handle this. 00:43:57.770 --> 00:44:01.220 The methods built into it aren't suitable. 00:44:01.220 --> 00:44:03.860 For generic DAEs, there are specific DAE solvers. 00:44:03.860 --> 00:44:06.020 So something like SUNDIALS or DAEPACK-- 00:44:06.020 --> 00:44:09.020 this is Professor Barton's software actually-- 00:44:09.020 --> 00:44:13.050 that are designed to handle DAEs up to some certain index 00:44:13.050 --> 00:44:13.580 instead. 00:44:13.580 --> 00:44:15.870 They're built to reliably solve these problems. 00:44:15.870 --> 00:44:18.500 So we have robots with constraints on them. 00:44:18.500 --> 00:44:21.779 We solve all sorts of problems for chemical process systems 00:44:21.779 --> 00:44:22.820 that are of higher index. 00:44:22.820 --> 00:44:25.589 And the way we do it is using specific numerical methods. 00:44:25.589 --> 00:44:27.380 They're based on the same sorts of schemes, 00:44:27.380 --> 00:44:30.020 backwards difference formulas, but careful analysis 00:44:30.020 --> 00:44:33.860 of the equations, which are in general unstructured, 00:44:33.860 --> 00:44:35.360 the software doesn't know beforehand 00:44:35.360 --> 00:44:38.630 what those equations look like, but a careful analysis of them 00:44:38.630 --> 00:44:41.180 in order to figure out how to minimize numerical errors 00:44:41.180 --> 00:44:45.620 and stay on the correct solution manifold. 00:44:45.620 --> 00:44:47.930 As input to these, though, we have 00:44:47.930 --> 00:44:49.440 to give initial conditions. 00:44:49.440 --> 00:44:52.675 So that's going to be the last thing that we talk about here. 00:44:52.675 --> 00:44:54.050 And those initial conditions have 00:44:54.050 --> 00:44:56.750 to be prescribed what we call consistently, 00:44:56.750 --> 00:45:00.320 or we can get numerical errors, or the software can just 00:45:00.320 --> 00:45:03.664 quit and throw an error, or it can run off 00:45:03.664 --> 00:45:06.080 on some other solution manifold that doesn't look anything 00:45:06.080 --> 00:45:08.060 like the solution we're after. 00:45:08.060 --> 00:45:11.360 The pendulum is an interesting example to think about. 00:45:11.360 --> 00:45:13.220 So you might say, well, what do I want 00:45:13.220 --> 00:45:14.712 to know about this pendulum. 00:45:14.712 --> 00:45:16.670 Maybe I want to know where it started initially 00:45:16.670 --> 00:45:20.030 and what its velocity initially was. 00:45:20.030 --> 00:45:23.240 Can I prescribe the initial position of the pendulum 00:45:23.240 --> 00:45:25.680 arbitrarily? 00:45:25.680 --> 00:45:27.358 What do you think yes or no? 00:45:27.358 --> 00:45:28.186 AUDIENCE: No. 00:45:28.186 --> 00:45:28.851 PROFESSOR: No. 00:45:28.851 --> 00:45:29.350 Why not? 00:45:29.350 --> 00:45:29.920 Why no? 00:45:29.920 --> 00:45:33.010 AUDIENCE: [INAUDIBLE] 00:45:33.010 --> 00:45:34.160 PROFESSOR: Good, yes. 00:45:34.160 --> 00:45:36.050 We know it's got a fixed arm length. 00:45:36.050 --> 00:45:37.069 Can't pick any position. 00:45:37.069 --> 00:45:38.110 That's not going to work. 00:45:38.110 --> 00:45:40.790 It's got to satisfy that algebraic constraint. 00:45:40.790 --> 00:45:45.110 Can I specify its velocity arbitrarily? 00:45:45.110 --> 00:45:46.130 What you think? 00:45:46.130 --> 00:45:48.426 Can I give it any initial velocity that mass 00:45:48.426 --> 00:45:49.550 at the end of the pendulum? 00:45:49.550 --> 00:45:50.350 AUDIENCE: No. 00:45:50.350 --> 00:45:51.001 PROFESSOR: No. 00:45:51.001 --> 00:45:51.500 Why not? 00:45:51.500 --> 00:45:52.448 AUDIENCE: [INAUDIBLE] 00:45:52.448 --> 00:45:53.870 PROFESSOR: What's that? 00:45:53.870 --> 00:45:56.714 AUDIENCE: [INAUDIBLE] 00:45:56.714 --> 00:45:58.540 PROFESSOR: Well, explain that to me, 00:45:58.540 --> 00:46:02.170 though, because I had equations for the velocity 00:46:02.170 --> 00:46:04.870 and for that acceleration and then 00:46:04.870 --> 00:46:08.390 an equation that told me that the length of this moment arm 00:46:08.390 --> 00:46:10.550 had a certain length associated with it. 00:46:10.550 --> 00:46:13.899 So why does that constrain the velocity? 00:46:13.899 --> 00:46:14.440 You're right. 00:46:14.440 --> 00:46:15.160 It does. 00:46:15.160 --> 00:46:16.326 I'm not saying you're wrong. 00:46:16.326 --> 00:46:17.370 You're right, but why? 00:46:17.370 --> 00:46:19.161 Those initial three equations don't tell me 00:46:19.161 --> 00:46:21.917 anything about the initial velocity was supposed to be. 00:46:21.917 --> 00:46:23.500 We know something about the trajectory 00:46:23.500 --> 00:46:25.083 this is supposed to sweep out, though. 00:46:25.083 --> 00:46:27.397 It's supposed to sweep out a circular trajectory. 00:46:27.397 --> 00:46:29.230 And the velocity in that circular trajectory 00:46:29.230 --> 00:46:32.470 is always going to be orthogonal to the moment time. 00:46:32.470 --> 00:46:35.770 And that equation actually-- it appeared. 00:46:35.770 --> 00:46:39.080 It appeared when we took a derivative of the constraint. 00:46:39.080 --> 00:46:41.950 So this was a hidden constraint. 00:46:41.950 --> 00:46:42.726 It's secret. 00:46:42.726 --> 00:46:45.100 You didn't know it was there until you took a derivative, 00:46:45.100 --> 00:46:47.980 and you tried to go and figure out what index this was. 00:46:47.980 --> 00:46:50.920 Then you discovered velocity and position or orthogonal 00:46:50.920 --> 00:46:52.070 to each other. 00:46:52.070 --> 00:46:55.150 My initial condition should respect that. 00:46:55.150 --> 00:46:56.050 How can they not? 00:46:56.050 --> 00:46:57.883 The initial conditions have to be a solution 00:46:57.883 --> 00:47:00.710 to the equations as well. 00:47:00.710 --> 00:47:05.110 So now can the initial stiffness be specified arbitrarily? 00:47:05.110 --> 00:47:06.590 No. 00:47:06.590 --> 00:47:08.650 There was also an algebraic equation 00:47:08.650 --> 00:47:10.600 that popped up at some point right here when 00:47:10.600 --> 00:47:12.433 I took the second derivative, which tells me 00:47:12.433 --> 00:47:13.990 how the stiffness has to be related 00:47:13.990 --> 00:47:15.520 to position and velocities. 00:47:15.520 --> 00:47:20.410 It's also hidden inside the structure of this DAE. 00:47:20.410 --> 00:47:22.866 So this-- again goes to model formulation, 00:47:22.866 --> 00:47:25.240 what are the initial conditions that I have to prescribe? 00:47:25.240 --> 00:47:26.860 If I prescribe the right ones, I'll 00:47:26.860 --> 00:47:30.670 get a solution that matches the dynamics of the problem I 00:47:30.670 --> 00:47:32.010 was actually interested in. 00:47:32.010 --> 00:47:33.970 If I prescribe them incorrectly, who 00:47:33.970 --> 00:47:35.830 knows what's going to result? 00:47:35.830 --> 00:47:38.352 We can't really predict. 00:47:38.352 --> 00:47:40.060 Depends on the solver we're working with. 00:47:42.820 --> 00:47:46.310 Here's a formal way to think of it. 00:47:46.310 --> 00:47:48.150 So usually in an initial value problem, 00:47:48.150 --> 00:47:51.380 usually we want to know what's-- 00:47:51.380 --> 00:47:53.510 for these first order difference-- 00:47:53.510 --> 00:47:57.290 differential equations, what's the initial value of the state, 00:47:57.290 --> 00:48:00.230 and what's the initial value of its derivative. 00:48:00.230 --> 00:48:05.559 That completely specifies what's going on at the beginning. 00:48:05.559 --> 00:48:07.100 One way to think about these problems 00:48:07.100 --> 00:48:19.420 is I've got this equation that I want to solve. 00:48:19.420 --> 00:48:22.381 I could think in some sense that x and x dot 00:48:22.381 --> 00:48:23.630 are independent of each other. 00:48:23.630 --> 00:48:30.840 I could always write this is something like an equation 00:48:30.840 --> 00:48:33.150 for f of x and z. 00:48:33.150 --> 00:48:35.340 And then I could add some extra information 00:48:35.340 --> 00:48:42.430 that tells me actually x dot is the same as z. 00:48:42.430 --> 00:48:46.040 So initially, I would like to know the values of x and z, 00:48:46.040 --> 00:48:47.340 that z is x dot. 00:48:47.340 --> 00:48:51.460 So I'd like to know initially the values of x dot and x. 00:48:51.460 --> 00:48:54.250 And tell me where I start, and these 00:48:54.250 --> 00:48:57.220 need to be consistent with the initial value problem 00:48:57.220 --> 00:48:58.030 that I specified. 00:48:58.030 --> 00:49:00.210 If I'm doing an index 0, an ODE IVP, at least 00:49:00.210 --> 00:49:01.210 be consistent with this. 00:49:01.210 --> 00:49:04.060 So if you tell me x naught, that I 00:49:04.060 --> 00:49:06.280 should put x into this equation, then that 00:49:06.280 --> 00:49:08.770 should tell me what x dot is. 00:49:08.770 --> 00:49:14.140 If you tell me x dot, then I should solve my governing 00:49:14.140 --> 00:49:17.710 equation for x naught. 00:49:17.710 --> 00:49:22.150 You could tell me some equation that relates x naught and x 00:49:22.150 --> 00:49:24.700 dot initially, and now I got to solve 00:49:24.700 --> 00:49:27.760 this equation in conjunction with my governing equation 00:49:27.760 --> 00:49:28.940 to determine these-- 00:49:28.940 --> 00:49:31.210 this set of values. 00:49:31.210 --> 00:49:34.600 I'd like to know both of these things initially. 00:49:34.600 --> 00:49:38.770 The fully-implicit DAE, I really have two n unknowns here. 00:49:38.770 --> 00:49:41.380 I don't know x, and I don't know x dot. 00:49:41.380 --> 00:49:44.530 And I've only got n equations for those. 00:49:44.530 --> 00:49:47.440 So apparently there's n degrees of freedom to specify here. 00:49:47.440 --> 00:49:49.180 I need n more equations to say what 00:49:49.180 --> 00:49:50.950 those initial conditions are. 00:49:50.950 --> 00:49:53.110 And these hidden constraints I point out actually 00:49:53.110 --> 00:49:55.690 reduce the number of degrees of freedom. 00:49:55.690 --> 00:49:59.650 The fact that these are not index 0 problems 00:49:59.650 --> 00:50:02.710 but have a finite index will introduce extra constraints 00:50:02.710 --> 00:50:05.050 that reduce the number of degrees of freedom 00:50:05.050 --> 00:50:08.734 or the number of other equations I can use to specify 00:50:08.734 --> 00:50:09.775 those initial conditions. 00:50:09.775 --> 00:50:11.020 Does that make sense? 00:50:11.020 --> 00:50:14.230 Let me give you some examples, and then I'll let you go. 00:50:14.230 --> 00:50:16.930 So if I have a problem that has separate differential 00:50:16.930 --> 00:50:18.400 states and algebraic states, this 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 00:50:23.291 --> 00:50:24.290 that I need to specify-- 00:50:24.290 --> 00:50:27.787 I need the initial derivative of the differential state, 00:50:27.787 --> 00:50:29.620 the initial value of the differential state, 00:50:29.620 --> 00:50:32.130 and the initial value of the algebraic state. 00:50:32.130 --> 00:50:34.930 I'd actually need to know the initial value of the derivative 00:50:34.930 --> 00:50:36.820 of the algebraic state here because I 00:50:36.820 --> 00:50:39.220 don't need that to satisfy this initial equation. 00:50:41.950 --> 00:50:44.310 Let's look at a couple of examples. 00:50:44.310 --> 00:50:47.400 So here's stirred tank example 1. 00:50:47.400 --> 00:50:50.050 I converted to a system of ODEs by taking the derivative 00:50:50.050 --> 00:50:51.680 of equation 2 here. 00:50:51.680 --> 00:50:53.740 So here's my two ODEs. 00:50:53.740 --> 00:50:55.682 I got to look back at my initial equations 00:50:55.682 --> 00:50:57.890 and see how do they constrain the initial conditions. 00:50:57.890 --> 00:51:01.080 So equation 2 tells me that c1 initially better 00:51:01.080 --> 00:51:02.675 be gamma initially. 00:51:02.675 --> 00:51:04.300 Otherwise, my initial condition doesn't 00:51:04.300 --> 00:51:08.080 satisfy the governing equations of the model that I wrote down. 00:51:08.080 --> 00:51:11.260 Equation 1 tells me that the derivative of c2 00:51:11.260 --> 00:51:14.740 initially better be related to the derivative of c1 00:51:14.740 --> 00:51:15.910 and the derivative of n-- 00:51:15.910 --> 00:51:20.650 I'm sorry to the value of c1 and the value of c2 initially. 00:51:20.650 --> 00:51:26.350 This third equation that I came up with, it doesn't constrain. 00:51:26.350 --> 00:51:29.470 It's a derivative of one of the original algebraic variables, 00:51:29.470 --> 00:51:31.877 so it doesn't constrain that original set of-- 00:51:31.877 --> 00:51:34.210 I need to know the initial condition on the differential 00:51:34.210 --> 00:51:36.824 variables, the derivative of the differential variables, 00:51:36.824 --> 00:51:37.990 and the algebraic variables. 00:51:37.990 --> 00:51:40.094 This only tells me something about the derivative 00:51:40.094 --> 00:51:41.260 of this algebraic variables. 00:51:41.260 --> 00:51:43.190 This equation doesn't give me any other constraints. 00:51:43.190 --> 00:51:45.310 So I just got a free variable that I can specify. 00:51:45.310 --> 00:51:46.160 One more equation. 00:51:46.160 --> 00:51:47.690 It can be whatever I want. 00:51:47.690 --> 00:51:50.500 So maybe I say c2 of 0c0, or maybe I 00:51:50.500 --> 00:51:54.160 say there's some relationship between c1 of 0 and c2 0. 00:51:54.160 --> 00:51:57.520 Or maybe I say that c2 of 0 and c2 prime of 0 00:51:57.520 --> 00:51:59.170 are related in some way. 00:51:59.170 --> 00:52:01.690 Whatever I specify, I have a system 00:52:01.690 --> 00:52:05.837 of three equations for three unknowns, c1 of 0, c2 dot of 0, 00:52:05.837 --> 00:52:07.270 and c2 of 0. 00:52:07.270 --> 00:52:09.220 I need to have a unique solution to that. 00:52:09.220 --> 00:52:11.890 Otherwise, there's going to be a problem solving for it in time, 00:52:11.890 --> 00:52:13.010 the system of equations. 00:52:13.010 --> 00:52:15.343 So I have to choose this one consistent with these other 00:52:15.343 --> 00:52:18.220 ones, but beyond that I'm OK. 00:52:22.549 --> 00:52:26.320 Here's stirred tank example 2. 00:52:26.320 --> 00:52:29.230 So I take a derivative of equation 2, 00:52:29.230 --> 00:52:32.053 and that gives me a differential equation for dc2 dt. 00:52:32.053 --> 00:52:33.790 Actually, I already had one of those, 00:52:33.790 --> 00:52:36.460 so I substitute for dc2 dt, and I 00:52:36.460 --> 00:52:41.130 take another derivative to get a differential equation for c1. 00:52:41.130 --> 00:52:43.070 And let's ask about the initial conditions. 00:52:43.070 --> 00:52:47.730 So the initial algebraic equation, it constrained dc2-- 00:52:47.730 --> 00:52:48.750 it constrained c2. 00:52:48.750 --> 00:52:52.360 It said c2 of 0 has to be gamma 0. 00:52:52.360 --> 00:52:57.130 The initial differential equation told me that c1 of 0 00:52:57.130 --> 00:53:01.670 had to be equal to c2 of 0 plus v over q, c2 dot of 0. 00:53:01.670 --> 00:53:05.260 Remember, I need three equations to specify c2 of 0, c1 of 0, 00:53:05.260 --> 00:53:07.520 c2 dot of 0. 00:53:07.520 --> 00:53:09.760 I actually can't specify c2 dot of 0 00:53:09.760 --> 00:53:14.200 freely because in the process of deriving this differential 00:53:14.200 --> 00:53:17.080 equation here so that I had dc1 dt and dc2 dt, 00:53:17.080 --> 00:53:23.590 I introduced a constraint on the derivatives of c2. 00:53:23.590 --> 00:53:28.150 And this tells me that c2 dot of 0 has to be gamma dot of 0. 00:53:28.150 --> 00:53:31.380 There are no free variables to specify. 00:53:31.380 --> 00:53:33.050 There's three equations. 00:53:33.050 --> 00:53:35.230 Two of them come from the initial system of DAEs, 00:53:35.230 --> 00:53:40.240 and the third one is related to this extra equation that 00:53:40.240 --> 00:53:42.850 popped up as we tried to define the index. 00:53:42.850 --> 00:53:46.370 It's an implicit constraint on the problem. 00:53:46.370 --> 00:53:48.060 Does that make sense? 00:53:48.060 --> 00:53:50.390 I can't specify these things freely. 00:53:50.390 --> 00:53:55.690 They're specified by the signal that I measured gamma. 00:53:55.690 --> 00:53:58.070 There are two examples here that you 00:53:58.070 --> 00:54:01.250 should be able to work through. 00:54:01.250 --> 00:54:03.695 One is an index 2 example, and one is an index 1 example. 00:54:03.695 --> 00:54:05.570 We're not going to have time to go over them. 00:54:05.570 --> 00:54:06.906 Sorry for that. 00:54:06.906 --> 00:54:08.030 But they work the same way. 00:54:08.030 --> 00:54:13.260 So convert to a system of ODEs, ask what initial conditions 00:54:13.260 --> 00:54:15.690 will satisfy the initial constraints and were there 00:54:15.690 --> 00:54:17.910 any other-- 00:54:17.910 --> 00:54:19.770 were there any other constraining equations 00:54:19.770 --> 00:54:20.784 that popped up? 00:54:20.784 --> 00:54:22.950 So if you work through this example, you'll find out 00:54:22.950 --> 00:54:27.370 there was a hidden strength that either the initial derivatives 00:54:27.370 --> 00:54:29.535 of c1 and c2 sum together have to be 0 00:54:29.535 --> 00:54:32.730 or equivalently the initial value of c3 00:54:32.730 --> 00:54:33.840 had to be equal to 0. 00:54:33.840 --> 00:54:35.400 You can take this constraint in either place. 00:54:35.400 --> 00:54:36.690 Maybe you want to take both of them 00:54:36.690 --> 00:54:38.400 to make sure that it's always satisfied, 00:54:38.400 --> 00:54:40.969 that there is no numerical error that drives you off 00:54:40.969 --> 00:54:41.760 of this constraint. 00:54:41.760 --> 00:54:43.860 It's OK to have these things over constrained. 00:54:43.860 --> 00:54:47.371 It's just not OK to have them under constrained. 00:54:47.371 --> 00:54:48.870 And you'll find out that while there 00:54:48.870 --> 00:54:52.140 were five things I could pick, the derivatives, and the values 00:54:52.140 --> 00:54:54.810 of c1 and c2, and then the algebraic variable c3 00:54:54.810 --> 00:54:55.330 initially. 00:54:55.330 --> 00:54:57.510 So there are only four constraints-- for equations 00:54:57.510 --> 00:54:58.890 here for the initial conditions. 00:54:58.890 --> 00:54:59.848 I get to pick one more. 00:54:59.848 --> 00:55:01.560 I can get it however I want. 00:55:01.560 --> 00:55:04.410 Just has to be consistent with these. 00:55:04.410 --> 00:55:05.740 There's one more example. 00:55:05.740 --> 00:55:08.377 Again, should work through these on your own. 00:55:08.377 --> 00:55:09.460 But you do the same thing. 00:55:09.460 --> 00:55:12.650 Convert to a system of ODEs. 00:55:12.650 --> 00:55:14.291 You'll find out that when you do this, 00:55:14.291 --> 00:55:16.040 you'll get a differential equation for c3. 00:55:16.040 --> 00:55:17.840 It didn't introduce any new constraints 00:55:17.840 --> 00:55:21.094 on c1 dot, c2 dot, or c1, c2, and c3. 00:55:21.094 --> 00:55:23.510 So you just need to satisfy-- your initial conditions have 00:55:23.510 --> 00:55:25.310 to satisfy the initial equations, 00:55:25.310 --> 00:55:28.130 and you can put in two more conditions, whatever 00:55:28.130 --> 00:55:29.210 you want them to be. 00:55:29.210 --> 00:55:31.130 Add five unknowns and only three equations 00:55:31.130 --> 00:55:33.170 for the initial conditions. 00:55:33.170 --> 00:55:34.500 And that's it. 00:55:34.500 --> 00:55:36.810 So thank you very much.