WEBVTT 00:00:01.580 --> 00:00:03.920 The following content is provided under a Creative 00:00:03.920 --> 00:00:05.340 Commons license. 00:00:05.340 --> 00:00:07.550 Your support will help MIT OpenCourseWare 00:00:07.550 --> 00:00:11.640 continue to offer high quality educational resources for free. 00:00:11.640 --> 00:00:14.180 To make a donation or to view additional materials 00:00:14.180 --> 00:00:18.110 from hundreds of MIT courses, visit MIT OpenCourseWare 00:00:18.110 --> 00:00:19.340 at ocw.mit.edu. 00:00:23.562 --> 00:00:25.270 WILLIAM GREEN, JR.: Shall we get started? 00:00:25.270 --> 00:00:27.340 Yes. 00:00:27.340 --> 00:00:29.710 Apologies for the notes being posted so late. 00:00:29.710 --> 00:00:33.020 I had trouble with the new two factor authentication system, 00:00:33.020 --> 00:00:35.980 which I need to get into Stellar and put things up there. 00:00:35.980 --> 00:00:39.610 My phone wasn't working so I had to find my iPad someplace 00:00:39.610 --> 00:00:42.542 and use that to get credentials to log in. 00:00:42.542 --> 00:00:44.500 So if you don't have them don't worry about it. 00:00:44.500 --> 00:00:47.540 The full notes are posted online, no missing pieces. 00:00:47.540 --> 00:00:49.480 So you can utilize those afterwards 00:00:49.480 --> 00:00:51.060 if something is unclear. 00:00:51.060 --> 00:00:51.560 All right? 00:00:56.250 --> 00:00:58.070 So your quizzes are graded. 00:00:58.070 --> 00:01:00.390 The TAs have them. 00:01:00.390 --> 00:01:02.610 I think overall, we were really pleased 00:01:02.610 --> 00:01:04.349 with the outcome of this quiz. 00:01:04.349 --> 00:01:06.360 It is very conceptual in nature and we 00:01:06.360 --> 00:01:08.790 thought people really understood a lot of the material 00:01:08.790 --> 00:01:09.540 quite well. 00:01:09.540 --> 00:01:12.440 And you can see that from the distribution of scores. 00:01:12.440 --> 00:01:15.650 So the mean on the quiz was about 77, 00:01:15.650 --> 00:01:17.250 a standard deviation of 12. 00:01:17.250 --> 00:01:19.930 I think everyone did really quite well. 00:01:19.930 --> 00:01:23.580 And on a number of the problems, perfect solutions were posted. 00:01:23.580 --> 00:01:26.670 So I'll take a second here. 00:01:26.670 --> 00:01:28.800 If there are any questions about the material that 00:01:28.800 --> 00:01:31.464 was on the quiz, or the format of it 00:01:31.464 --> 00:01:32.880 that you'd like to have addressed, 00:01:32.880 --> 00:01:35.351 I'm happy to answer them. 00:01:35.351 --> 00:01:35.850 No? 00:01:38.580 --> 00:01:41.440 Does it seem like it took too much time 00:01:41.440 --> 00:01:44.080 to complete all three questions in the two hours allotted? 00:01:44.080 --> 00:01:49.510 Or was it about the right length of material for the time slot? 00:01:49.510 --> 00:01:51.030 So we vetted it with the two TAs. 00:01:51.030 --> 00:01:53.720 And it took each of them less than an hour to complete it. 00:01:53.720 --> 00:01:55.969 So that was a good sign that we thought 00:01:55.969 --> 00:01:57.760 you guys should be able to get through most 00:01:57.760 --> 00:01:59.140 of the material in two hours. 00:02:05.560 --> 00:02:08.012 So let's recap. 00:02:08.012 --> 00:02:10.470 We're going to move on today to a topic called differential 00:02:10.470 --> 00:02:11.650 algebraic equations. 00:02:11.650 --> 00:02:14.310 But before doing that, there are a couple of concepts 00:02:14.310 --> 00:02:16.260 I want to review from numerical integration, 00:02:16.260 --> 00:02:19.669 or actually concepts that weren't covered when we talked 00:02:19.669 --> 00:02:21.210 about numerical integration on Friday 00:02:21.210 --> 00:02:22.480 that I think are important. 00:02:22.480 --> 00:02:23.855 And then I want to review briefly 00:02:23.855 --> 00:02:26.190 implicit methods for ODE-IVPs because those 00:02:26.190 --> 00:02:28.710 are going to be important for differential algebraic 00:02:28.710 --> 00:02:29.250 equations. 00:02:29.250 --> 00:02:32.550 So for numerical integration we talked 00:02:32.550 --> 00:02:35.250 about quadrature of various integrals, 00:02:35.250 --> 00:02:36.960 how to develop quadrature formulas. 00:02:36.960 --> 00:02:38.130 But there's actually one type of integral 00:02:38.130 --> 00:02:39.150 we didn't talk about there, which 00:02:39.150 --> 00:02:41.910 is sort of improper integrals, where the limits of integration 00:02:41.910 --> 00:02:43.440 are unbounded. 00:02:43.440 --> 00:02:46.670 And these come up all the time in engineering problems. 00:02:46.670 --> 00:02:49.560 We'd like to know what the, say, integrated 00:02:49.560 --> 00:02:53.576 response of some function is over all possible times. 00:02:53.576 --> 00:02:55.450 And we may have to evaluate that numerically. 00:02:55.450 --> 00:02:57.860 So how do you do those sorts of integrals? 00:02:57.860 --> 00:03:00.090 And it turns out the strategy that's most often used 00:03:00.090 --> 00:03:03.420 is to divide this domain up into two pieces. 00:03:03.420 --> 00:03:05.370 One piece, which is proper integral 00:03:05.370 --> 00:03:07.680 over a finite domain, and another piece, 00:03:07.680 --> 00:03:09.930 which is an improper integral over an infinite domain. 00:03:09.930 --> 00:03:11.513 This first integral, you can handle it 00:03:11.513 --> 00:03:15.120 with an ODE-IVP method or some higher order polynomial 00:03:15.120 --> 00:03:16.560 interpolation. 00:03:16.560 --> 00:03:19.494 But the second one we have to do differently. 00:03:19.494 --> 00:03:21.160 And there are couple of ways to do this. 00:03:21.160 --> 00:03:23.370 One is to transform variables. 00:03:23.370 --> 00:03:27.070 So you map this domain onto a finite domain. 00:03:27.070 --> 00:03:28.920 And the other option is to substitute 00:03:28.920 --> 00:03:30.420 an asymptotic approximation. 00:03:30.420 --> 00:03:34.740 So we say, when, say, tau is large, f of tau 00:03:34.740 --> 00:03:37.410 has a characteristic asymptotic approximation. 00:03:37.410 --> 00:03:40.200 Maybe it's exponentially decaying. 00:03:40.200 --> 00:03:43.080 And we take the integral of that asymptotic approximation. 00:03:43.080 --> 00:03:44.910 And the more accurate that approximation 00:03:44.910 --> 00:03:48.150 is the better our approximation for this whole 00:03:48.150 --> 00:03:49.800 integral will be. 00:03:49.800 --> 00:03:52.370 There's another type of interval that we have to do, 00:03:52.370 --> 00:03:54.960 where the same idea can be applied. 00:03:54.960 --> 00:03:58.485 That's the integral over integrable singularities. 00:03:58.485 --> 00:03:59.360 So here's an example. 00:04:02.810 --> 00:04:06.300 So say we want to integrate the function cosine of tau 00:04:06.300 --> 00:04:10.185 over square root of tau from 0 to some finite positive time 00:04:10.185 --> 00:04:12.800 point. 00:04:12.800 --> 00:04:15.260 As tau goes to 0, cosine goes to 1 and square root 00:04:15.260 --> 00:04:17.430 of tau diverges. 00:04:17.430 --> 00:04:20.390 But this integral actually has a finite value. 00:04:20.390 --> 00:04:21.980 This 1 over square root of tau, that's 00:04:21.980 --> 00:04:25.650 an integrable singularity. 00:04:25.650 --> 00:04:26.670 So how do you do this? 00:04:26.670 --> 00:04:29.730 Well, you want to split the domain again into two parts. 00:04:29.730 --> 00:04:32.840 One part that contains the integrable singularity 00:04:32.840 --> 00:04:35.270 and the other part, which excludes it. 00:04:35.270 --> 00:04:38.210 And in the part that contains the singularity 00:04:38.210 --> 00:04:41.540 we might do an asymptotic expansion of the function, 00:04:41.540 --> 00:04:44.510 and integrate each of those asymptotically accurate terms 00:04:44.510 --> 00:04:45.740 separately. 00:04:45.740 --> 00:04:48.530 So the integral of 1 over square root of tau 00:04:48.530 --> 00:04:51.680 is going to give us 2 tau 0 to the 1/2. 00:04:51.680 --> 00:04:54.230 The integral of this ratio here is 00:04:54.230 --> 00:04:57.840 going to give us minus 1/2 t0 to the 5/2. 00:04:57.840 --> 00:05:00.600 And then we have an integral over a finite domain, 00:05:00.600 --> 00:05:02.580 which doesn't include the singularity. 00:05:02.580 --> 00:05:05.390 And the accuracy of this method can be dictated by two things 00:05:05.390 --> 00:05:06.140 now. 00:05:06.140 --> 00:05:08.370 One is how accurately can we do this integral. 00:05:08.370 --> 00:05:13.190 And the other is how accurate is our asymptotic expansion here. 00:05:13.190 --> 00:05:16.130 So if we want higher accuracy we need more terms. 00:05:16.130 --> 00:05:17.810 We might need to be selective about how 00:05:17.810 --> 00:05:20.000 we choose the end point for this domain 00:05:20.000 --> 00:05:21.489 in order to minimize the error. 00:05:21.489 --> 00:05:23.030 There are lots of methods to do this. 00:05:23.030 --> 00:05:23.490 Sam. 00:05:23.490 --> 00:05:25.031 AUDIENCE: Will you talk a little more 00:05:25.031 --> 00:05:28.360 about what typical integral [INAUDIBLE] is? 00:05:28.360 --> 00:05:31.210 WILLIAM GREEN, JR.: Yeah. 00:05:31.210 --> 00:05:33.850 So you know these things. 00:05:33.850 --> 00:05:40.260 If I try to integrate 1 over x, as x goes to 0 00:05:40.260 --> 00:05:42.310 this thing will always diverge. 00:05:42.310 --> 00:05:46.180 You know the antiderivative of 1 over x is going to be the log. 00:05:46.180 --> 00:05:47.890 And the log diverges. 00:05:47.890 --> 00:05:53.950 So power is weaker than 1 over x, like 1 over x to the 1/2. 00:05:53.950 --> 00:05:55.290 Don't diverge like the log. 00:05:55.290 --> 00:05:56.830 They actually go to a finite value. 00:05:56.830 --> 00:05:59.770 The log is like the most weakly singular function there is. 00:05:59.770 --> 00:06:02.800 And integrals over weaker powers of 1 over x actually 00:06:02.800 --> 00:06:04.994 don't produce a singularity anymore. 00:06:04.994 --> 00:06:06.910 They're what we call integrable singularities. 00:06:06.910 --> 00:06:11.200 The function itself diverges but its integral does not. 00:06:11.200 --> 00:06:15.100 So you usually have this asymptotic power law type 00:06:15.100 --> 00:06:16.270 behavior. 00:06:16.270 --> 00:06:18.140 OK. 00:06:18.140 --> 00:06:20.100 So they come up all the time. 00:06:20.100 --> 00:06:23.910 We see these things in places like, say, 00:06:23.910 --> 00:06:27.460 squeezing a thin film of fluid between two plates. 00:06:27.460 --> 00:06:31.567 How long does it take for these two plates to come together? 00:06:31.567 --> 00:06:33.900 You might say, well, as the plates get closer and closer 00:06:33.900 --> 00:06:36.737 together the pressure starts to diverge in the gap 00:06:36.737 --> 00:06:38.070 and they'll never come together. 00:06:38.070 --> 00:06:40.950 But maybe depending on the geometry of the plates 00:06:40.950 --> 00:06:42.750 there can be a sort of finite time 00:06:42.750 --> 00:06:44.190 at which they'll come together. 00:06:44.190 --> 00:06:45.273 Depends on their geometry. 00:06:48.557 --> 00:06:49.390 So that's one topic. 00:06:49.390 --> 00:06:51.473 I think that's useful because these things come up 00:06:51.473 --> 00:06:52.030 all the time. 00:06:52.030 --> 00:06:55.330 And you can't actually handle those cases with quadrature 00:06:55.330 --> 00:06:56.800 by itself. 00:06:56.800 --> 00:06:57.910 The function diverges. 00:06:57.910 --> 00:06:59.680 There's no polynomial interpolation 00:06:59.680 --> 00:07:01.150 that's suitable to match it. 00:07:01.150 --> 00:07:02.777 So the quadrature can't handle it. 00:07:02.777 --> 00:07:04.360 And if you're trying to integrate over 00:07:04.360 --> 00:07:07.210 an infinite domain, well, good luck. 00:07:07.210 --> 00:07:09.580 You're never going to be able to fit enough points in 00:07:09.580 --> 00:07:14.440 to know that you've accurately integrated out to infinity. 00:07:14.440 --> 00:07:17.260 The other thing I want to recap here, 00:07:17.260 --> 00:07:19.750 and this is a problem for you to try. 00:07:19.750 --> 00:07:23.890 So here we have a simple first order ordinary differential 00:07:23.890 --> 00:07:24.390 equation. 00:07:24.390 --> 00:07:27.080 We want to use the implicit Euler method to solve it. 00:07:27.080 --> 00:07:29.840 So can you give a closed form solution 00:07:29.840 --> 00:07:33.260 for one step of the implicit Euler method, 00:07:33.260 --> 00:07:36.620 or for n steps of the implicit Euler method? 00:07:36.620 --> 00:07:38.300 What result does that produce? 00:07:38.300 --> 00:07:40.520 Can you compute this? 00:07:40.520 --> 00:07:48.556 Do you guys talk about backward Euler methods for for ODE-IVPs? 00:07:48.556 --> 00:07:51.180 Backwards difference methods in general and the backwards Euler 00:07:51.180 --> 00:07:54.600 method specifically? 00:07:54.600 --> 00:07:55.940 No. 00:07:55.940 --> 00:07:57.150 OK. 00:07:57.150 --> 00:08:01.770 Well, clearly I should have been here. 00:08:01.770 --> 00:08:03.150 So what's the strategy? 00:08:03.150 --> 00:08:05.100 We have to approximate the derivative. 00:08:05.100 --> 00:08:06.960 Yes? 00:08:06.960 --> 00:08:08.880 So the approximation for the derivative we use 00:08:08.880 --> 00:08:11.160 is a finite difference approximation. 00:08:11.160 --> 00:08:14.940 So we write dx dt as x at a point k 00:08:14.940 --> 00:08:18.870 plus 1 minus x at the point k divided by delta t. 00:08:18.870 --> 00:08:20.700 Backwards difference. 00:08:20.700 --> 00:08:24.930 We evaluate the right-hand side of the ODE-IVP at k 00:08:24.930 --> 00:08:27.940 plus 1, xk plus 1. 00:08:27.940 --> 00:08:30.360 And then we solve for xk plus 1. 00:08:30.360 --> 00:08:32.159 We just substitute in here. 00:08:32.159 --> 00:08:33.690 Backwards difference approximation 00:08:33.690 --> 00:08:35.320 for the differential. 00:08:35.320 --> 00:08:37.979 And xk plus 1 for the right-hand side. 00:08:37.979 --> 00:08:41.010 And we find that xk plus 1 will be 1 over 1 00:08:41.010 --> 00:08:46.950 minus delta t times this lambda times xk. 00:08:46.950 --> 00:08:51.680 And if we iterate this from k equals 0 up to some finite k, 00:08:51.680 --> 00:08:54.240 it's like multiplying by powers of 1 over 1 00:08:54.240 --> 00:08:55.980 minus delta t times lambda. 00:08:55.980 --> 00:08:57.044 Yes? 00:08:57.044 --> 00:08:57.960 Does this makes sense? 00:09:01.240 --> 00:09:05.910 So our backwards Euler solution to this fundamental problem 00:09:05.910 --> 00:09:08.386 is going to look like this. 00:09:08.386 --> 00:09:10.260 And a lot of times we ask about the stability 00:09:10.260 --> 00:09:12.240 of these sorts of solutions. 00:09:12.240 --> 00:09:16.850 So if for example, the solution to this equation 00:09:16.850 --> 00:09:20.000 was supposed to decay to zero, we 00:09:20.000 --> 00:09:21.770 would expect our numerical method 00:09:21.770 --> 00:09:23.890 to also yield results which decay to zero. 00:09:23.890 --> 00:09:27.110 We would call that stable. 00:09:27.110 --> 00:09:29.660 So we all know under what conditions 00:09:29.660 --> 00:09:34.662 does xk, for example, go to 0 as I change the product delta 00:09:34.662 --> 00:09:36.350 t times lambda. 00:09:36.350 --> 00:09:40.850 Lambda could be a real number or it can be a complex number too. 00:09:40.850 --> 00:09:42.360 Doesn't really matter. 00:09:42.360 --> 00:09:46.830 So what needs to be true for xk to the k to 0 as k gets bigger? 00:09:46.830 --> 00:09:49.340 Well what needs to be true is this quantity in parentheses 00:09:49.340 --> 00:09:50.960 needs to be smaller than 1 because I'm 00:09:50.960 --> 00:09:52.400 taking lots of products of things 00:09:52.400 --> 00:09:53.483 that are smaller than one. 00:09:53.483 --> 00:09:57.770 That will continue to shrink xk from x 0 until it goes to zero. 00:09:57.770 --> 00:10:01.790 If this is smaller than 1, then this quantity down here 00:10:01.790 --> 00:10:04.610 needs to be bigger than 1 in absolute value. 00:10:04.610 --> 00:10:07.010 And then I'll get solutions that slowly decay to zero. 00:10:07.010 --> 00:10:10.490 So if I ask, for what values of delta t times lambda 00:10:10.490 --> 00:10:13.770 is the absolute value of this thing bigger than 1, 00:10:13.770 --> 00:10:17.245 and I allow lambda to be a complex number, for example? 00:10:17.245 --> 00:10:20.110 I'll find out that's true when 1 minus delta t 00:10:20.110 --> 00:10:22.490 times the real part of lambda squared 00:10:22.490 --> 00:10:26.630 plus delta t times the imaginary part of lambda squared 00:10:26.630 --> 00:10:27.530 is bigger than 1. 00:10:30.100 --> 00:10:33.440 And if I plot the real part of delta t times 00:10:33.440 --> 00:10:37.630 lambda versus the imaginary part of delta times lambda, 00:10:37.630 --> 00:10:41.160 this inequality is represented by this pink zone 00:10:41.160 --> 00:10:42.470 on the outside of the circle. 00:10:42.470 --> 00:10:44.510 So any values of delta t times lambda 00:10:44.510 --> 00:10:49.460 that live in this pink area will give stable integration. 00:10:49.460 --> 00:10:54.500 And the solution xk, well, the k to zero as k goes to infinity. 00:10:57.020 --> 00:10:59.530 So it makes sense? 00:10:59.530 --> 00:11:02.010 This is a little funny though. 00:11:02.010 --> 00:11:03.780 What does the solution to this equation 00:11:03.780 --> 00:11:08.012 do when lambda is bigger than 1? 00:11:08.012 --> 00:11:09.470 When it's bigger than 0, excuse me. 00:11:09.470 --> 00:11:10.886 When lambda is bigger than 0, what 00:11:10.886 --> 00:11:14.860 does the solution to this ODE-IVP do? 00:11:14.860 --> 00:11:15.560 Does it decay? 00:11:15.560 --> 00:11:16.640 AUDIENCE: [INAUDIBLE] 00:11:16.640 --> 00:11:17.973 WILLIAM GREEN, JR.: It blows up. 00:11:17.973 --> 00:11:19.100 It grows exponentially. 00:11:19.100 --> 00:11:21.920 But the numerical method here, when 00:11:21.920 --> 00:11:25.070 I have lambda is bigger than zero, 00:11:25.070 --> 00:11:28.659 can be stable and produce decaying solutions. 00:11:28.659 --> 00:11:29.825 So this is sort of peculiar. 00:11:32.890 --> 00:11:35.290 So we're often concerned when we solve ODE-IVPs 00:11:35.290 --> 00:11:37.640 and we try to use these sorts of backwards difference 00:11:37.640 --> 00:11:40.880 methods, these implicit methods, to solve them. 00:11:40.880 --> 00:11:44.390 We're often concerned with getting stability. 00:11:44.390 --> 00:11:46.299 We want to get solutions that decay 00:11:46.299 --> 00:11:47.590 when they're supposed to decay. 00:11:47.590 --> 00:11:49.340 Sometimes we get solutions that decay even 00:11:49.340 --> 00:11:52.690 when the solution is supposed to blow up. 00:11:52.690 --> 00:11:55.240 So that's sort of funny. 00:11:55.240 --> 00:11:57.835 So the exact solution should look like this. 00:11:57.835 --> 00:11:59.710 But there are going to be circumstances where 00:11:59.710 --> 00:12:00.940 lambda is bigger than zero. 00:12:00.940 --> 00:12:02.530 In fact, lambda is bigger than-- 00:12:02.530 --> 00:12:06.340 lambda delta t is bigger than 2 along the real axis, 00:12:06.340 --> 00:12:11.060 where the solution will actually blow up rather than decay away. 00:12:11.060 --> 00:12:12.970 So stability and accuracy don't always 00:12:12.970 --> 00:12:15.760 correlate with each other. 00:12:15.760 --> 00:12:18.760 We really have to understand the underlying dynamics 00:12:18.760 --> 00:12:21.040 of the problem we're trying to solve 00:12:21.040 --> 00:12:24.110 before we choose a numerical method to apply to it. 00:12:24.110 --> 00:12:27.280 And this is going to be true of differential algebraic 00:12:27.280 --> 00:12:29.740 equations as well. 00:12:29.740 --> 00:12:33.310 Are there any questions about this? 00:12:33.310 --> 00:12:37.030 There's a discussion of these things in your textbook, 00:12:37.030 --> 00:12:39.160 these sorts of stability regions associated 00:12:39.160 --> 00:12:41.290 with different backwards difference methods. 00:12:41.290 --> 00:12:46.450 The implicit Euler method is one backwards difference formula. 00:12:46.450 --> 00:12:48.760 It's the simplest one. 00:12:48.760 --> 00:12:52.160 It's the least accurate one that you can derive 00:12:52.160 --> 00:12:53.420 but it's an easy one to use. 00:12:56.074 --> 00:12:58.490 But now we want to move on to a different sort of problem. 00:12:58.490 --> 00:13:01.737 And it's one that we come across in engineering all the time. 00:13:01.737 --> 00:13:03.320 These problems are called differential 00:13:03.320 --> 00:13:04.670 algebraic equations. 00:13:04.670 --> 00:13:08.270 And they're problems of the sort a function of some state 00:13:08.270 --> 00:13:12.500 variables, and the derivatives of those state variables 00:13:12.500 --> 00:13:15.147 and time is equal to zero. 00:13:15.147 --> 00:13:16.730 And there are some initial conditions. 00:13:16.730 --> 00:13:19.825 The state variables at that time 0 is equal to some x0. 00:13:24.840 --> 00:13:26.490 So this is a vector valued function 00:13:26.490 --> 00:13:29.320 of a vector valued state and it's time derivatives. 00:13:29.320 --> 00:13:31.650 And we want to understand the dynamics of this system. 00:13:31.650 --> 00:13:33.660 How does x vary with time in a way 00:13:33.660 --> 00:13:37.500 that's consistent with this governing equation? 00:13:37.500 --> 00:13:39.840 Usually we call these well-posed problems 00:13:39.840 --> 00:13:42.690 when the dimension of the state variable and the function 00:13:42.690 --> 00:13:43.590 are the same. 00:13:43.590 --> 00:13:46.730 They both live in the vector space rn. 00:13:46.730 --> 00:13:48.375 So we have n different states. 00:13:48.375 --> 00:13:50.550 I have n different functions that I have to satisfy. 00:13:50.550 --> 00:13:55.090 And I want to understand how x varies with time. 00:13:55.090 --> 00:13:57.320 Here's an example. 00:13:57.320 --> 00:13:58.200 A stirred tank. 00:13:58.200 --> 00:14:00.460 I call this stirred tank, example one. 00:14:00.460 --> 00:14:06.684 So into the stirred tank comes some concentration of solute. 00:14:06.684 --> 00:14:09.100 Out of the stirred tank comes some different concentration 00:14:09.100 --> 00:14:09.891 of the same solute. 00:14:09.891 --> 00:14:14.120 And we want to model the dynamics as a function of time. 00:14:14.120 --> 00:14:16.900 So we know that a material balance 00:14:16.900 --> 00:14:19.390 of the solute on the stirred tank, if it's well mixed, 00:14:19.390 --> 00:14:26.320 will tell us that d c2 dt is related to the volumetric flow 00:14:26.320 --> 00:14:28.960 rate into the tank divided by the volume of the tank 00:14:28.960 --> 00:14:32.140 multiplied by the difference between c1 and c2, 00:14:32.140 --> 00:14:37.112 the amount carried in and the amount carried out of the tank. 00:14:37.112 --> 00:14:38.890 And let's put a constraint in here too. 00:14:38.890 --> 00:14:42.430 Let's suppose we're trying to control this tank in some way. 00:14:42.430 --> 00:14:45.820 So we say that c1 of time has got 00:14:45.820 --> 00:14:49.190 to be equal to some control signal, gamma of t. 00:14:52.390 --> 00:14:55.810 And then we want to understand the dynamics of this system. 00:14:55.810 --> 00:14:59.710 At time 0, c2 is some c0, the initial concentration 00:14:59.710 --> 00:15:00.570 of the tank. 00:15:00.570 --> 00:15:03.490 And c1 is whatever the initial value of this control signal 00:15:03.490 --> 00:15:04.960 is. 00:15:04.960 --> 00:15:08.350 And the solution is c1 equals gamma of t. 00:15:08.350 --> 00:15:11.470 And c2 is, well, the initial concentration's 00:15:11.470 --> 00:15:13.090 got to decay out exponentially. 00:15:13.090 --> 00:15:15.220 And then I've got to integrate over how 00:15:15.220 --> 00:15:17.770 the input is changing in time. 00:15:17.770 --> 00:15:20.060 Those also produce exponential decay. 00:15:20.060 --> 00:15:22.870 So if the input is a little delta function spike 00:15:22.870 --> 00:15:25.330 then I'll get an exponential decay leaving the tank. 00:15:25.330 --> 00:15:27.310 If the input is changing dynamically in time, 00:15:27.310 --> 00:15:30.550 I need this convolution of the input with this integral. 00:15:30.550 --> 00:15:33.130 So you just to solve this system, 00:15:33.130 --> 00:15:34.900 this ordinary differential equation 00:15:34.900 --> 00:15:37.310 by substituting c1 in there. 00:15:37.310 --> 00:15:40.280 But this system of equations here is not just 00:15:40.280 --> 00:15:42.040 a differential equation. 00:15:42.040 --> 00:15:45.847 It's a differential equation and an algebraic equation. 00:15:45.847 --> 00:15:47.680 There's something peculiar about that that's 00:15:47.680 --> 00:15:50.140 different than just solving systems 00:15:50.140 --> 00:15:52.440 of ordinary differential equations. 00:15:52.440 --> 00:15:57.400 Oftentimes in models, we write down all 00:15:57.400 --> 00:15:59.720 the governing equations associated with the model. 00:15:59.720 --> 00:16:01.636 Some of them may be differential, some of them 00:16:01.636 --> 00:16:02.800 may be algebraic. 00:16:02.800 --> 00:16:04.960 There are times when we can look at those equations 00:16:04.960 --> 00:16:07.480 and see clever substitutions that we can make. 00:16:07.480 --> 00:16:10.840 But if you have a hundred equations, or a thousand 00:16:10.840 --> 00:16:12.900 equations, or a million equations, 00:16:12.900 --> 00:16:14.820 there's no way to do that reliably. 00:16:14.820 --> 00:16:18.310 We just kind of wind up with a system of ordinary differential 00:16:18.310 --> 00:16:21.730 equations and algebraic equations mixed together. 00:16:21.730 --> 00:16:24.777 And we have to solve these reliably. 00:16:24.777 --> 00:16:25.735 Here's another example. 00:16:28.650 --> 00:16:29.720 Oh, excuse me. 00:16:29.720 --> 00:16:31.800 So let me write this in this form. 00:16:31.800 --> 00:16:34.110 So we said there should be a vector valued 00:16:34.110 --> 00:16:38.460 function of the states and the derivatives of the states. 00:16:38.460 --> 00:16:41.490 The states are the concentration c1 and c2. 00:16:41.490 --> 00:16:45.090 And the derivative here is dc 2 dt. 00:16:45.090 --> 00:16:46.720 So here's my vector valued function. 00:16:46.720 --> 00:16:48.900 The first element is this, and it's equal to 0. 00:16:48.900 --> 00:16:51.660 The second element is this, and it's equal to 0. 00:16:51.660 --> 00:16:54.900 Here's my initial condition over the states. 00:16:54.900 --> 00:16:57.930 So I can just transform that list of equations 00:16:57.930 --> 00:17:00.450 to the sort of canonical form for a differential 00:17:00.450 --> 00:17:01.420 algebraic equation. 00:17:04.690 --> 00:17:07.619 Here's a separate example. 00:17:07.619 --> 00:17:12.790 Suppose instead, I'm trying to either control 00:17:12.790 --> 00:17:16.450 or make a measurement of the outlet concentrations c2. 00:17:16.450 --> 00:17:18.220 And we call that gamma of t. 00:17:22.060 --> 00:17:26.329 And I want to know now what is c2 and what is c1? 00:17:26.329 --> 00:17:29.340 So I've got to solve this system of equations for c2 and c1. 00:17:29.340 --> 00:17:31.740 The solution for c2 is easy. 00:17:31.740 --> 00:17:32.360 I know that. 00:17:32.360 --> 00:17:35.060 c2 is gamma, by definition. 00:17:35.060 --> 00:17:37.430 I got to substitute this into the first equation 00:17:37.430 --> 00:17:38.420 and solve for c1. 00:17:38.420 --> 00:17:43.199 And I see, c1 is gamma plus v over q gamma dot. 00:17:43.199 --> 00:17:45.240 And I have to have initial conditions to go along 00:17:45.240 --> 00:17:47.992 with this. 00:17:47.992 --> 00:17:49.200 There's something funny here. 00:17:49.200 --> 00:17:52.800 We can see the initial condition for c2 has to be gamma 0. 00:17:52.800 --> 00:17:57.710 The initial condition for c1, what does that have to be? 00:17:57.710 --> 00:18:01.430 Well, it needs to be consistent with the solution here. 00:18:01.430 --> 00:18:03.800 There were no free variables when 00:18:03.800 --> 00:18:05.484 I solved this equation for c1. 00:18:05.484 --> 00:18:07.400 It isn't like solving a differential equation, 00:18:07.400 --> 00:18:10.220 having a free coefficient to specify. 00:18:10.220 --> 00:18:13.130 So it better be the case that this c0 here 00:18:13.130 --> 00:18:18.680 is the same as gamma 0 plus v over q gamma dot 0. 00:18:18.680 --> 00:18:21.410 Somehow I have to know this input signal to prescribe 00:18:21.410 --> 00:18:23.030 the initial conditions for c1. 00:18:23.030 --> 00:18:28.805 So that's peculiar and different from just ODE-IVPs. 00:18:28.805 --> 00:18:29.680 You see this picture? 00:18:29.680 --> 00:18:32.940 You see how this goes together? 00:18:32.940 --> 00:18:34.630 The solution is funny too. 00:18:34.630 --> 00:18:36.880 Suppose we are trying to use this system of equations 00:18:36.880 --> 00:18:39.040 to do the following. 00:18:39.040 --> 00:18:43.720 We are trying to measure c2 and use the measurement of c2 00:18:43.720 --> 00:18:46.030 to predict c1. 00:18:46.030 --> 00:18:48.580 So gamma is the measurement, and we're 00:18:48.580 --> 00:18:52.060 trying to solve this system of differential algebraic 00:18:52.060 --> 00:18:55.120 equations to predict what c1 is. 00:18:55.120 --> 00:18:58.300 All measurements incur numerical error. 00:18:58.300 --> 00:19:01.560 So even though our signal for gamma that we measure 00:19:01.560 --> 00:19:06.670 may be continuous, it's bouncing around wildly. 00:19:06.670 --> 00:19:08.290 It's not a constant signal. 00:19:08.290 --> 00:19:10.900 It's going to move around a great deal. 00:19:10.900 --> 00:19:13.960 And c1, our predicted value for c1, 00:19:13.960 --> 00:19:15.670 depends on the derivatives of gamma, 00:19:15.670 --> 00:19:18.490 which means c1 is going to be incredibly sensitive to how 00:19:18.490 --> 00:19:21.440 that measurement is made. 00:19:21.440 --> 00:19:22.395 So it's peculiar. 00:19:22.395 --> 00:19:24.520 It means that there's going to be a lot of problems 00:19:24.520 --> 00:19:28.840 with stably solving these equations because the solutions 00:19:28.840 --> 00:19:34.390 can admit cases where they're not integrals of the input 00:19:34.390 --> 00:19:37.990 but they're actually related to derivatives of the input. 00:19:37.990 --> 00:19:39.960 Look at the solution back here again. 00:19:39.960 --> 00:19:43.140 Oop, sorry. 00:19:43.140 --> 00:19:45.960 The solution here was related to an integral of the input. 00:19:45.960 --> 00:19:48.030 Suppose I was making a measurement of c1 00:19:48.030 --> 00:19:50.940 in trying to predict c2 instead? 00:19:50.940 --> 00:19:53.480 My prediction for c2 would be related to the integral 00:19:53.480 --> 00:19:54.480 of that measured signal. 00:19:54.480 --> 00:19:56.850 And we know integrals smooth signals out. 00:19:56.850 --> 00:19:59.340 So it's not going to be so sensitive to the measurement. 00:19:59.340 --> 00:20:01.740 So one way of doing this seems really sensible. 00:20:01.740 --> 00:20:06.540 And the other way of doing it seems a little bizarre. 00:20:06.540 --> 00:20:09.317 Well, we can construct these equations however we want. 00:20:09.317 --> 00:20:11.400 So there's something about how we formulate models 00:20:11.400 --> 00:20:13.140 that are going to make them either easily 00:20:13.140 --> 00:20:18.850 solvable, as DAEs, or quite difficult to solve. 00:20:18.850 --> 00:20:23.190 So these pop up all the time in engineering problems. 00:20:23.190 --> 00:20:25.530 They pop up in dynamic simulations, 00:20:25.530 --> 00:20:27.900 where we have some sort of underlying conservation 00:20:27.900 --> 00:20:31.950 principle, or constraints, or equilibria. 00:20:31.950 --> 00:20:33.450 So if we have to conserve something, 00:20:33.450 --> 00:20:38.550 like total energy, total mass, total momentum, total particle 00:20:38.550 --> 00:20:41.760 number, total atomic species, total charge, 00:20:41.760 --> 00:20:44.040 there will always be an algebraic constraint that 00:20:44.040 --> 00:20:48.210 tells us the total mass has to stay this, 00:20:48.210 --> 00:20:51.500 or the total number of atoms of a certain type has to say this. 00:20:51.500 --> 00:20:52.230 Yeah, Ian. 00:20:52.230 --> 00:20:56.950 AUDIENCE: Just to be clear, were you saying, on the stirred tank 00:20:56.950 --> 00:20:59.722 example, that those were the same problem 00:20:59.722 --> 00:21:02.119 with different solution approaches? 00:21:02.119 --> 00:21:03.660 WILLIAM GREEN, JR.: So good question. 00:21:03.660 --> 00:21:04.285 Let me go back. 00:21:07.260 --> 00:21:10.960 So they are fundamentally different problems. 00:21:10.960 --> 00:21:17.450 So in one case I specified the value of c2 to be this gamma. 00:21:17.450 --> 00:21:21.450 I tried to make a measurement of c2 and then predict c1. 00:21:21.450 --> 00:21:24.180 In the other case, I specified c1. 00:21:24.180 --> 00:21:28.100 I tried to measure c1 and predict c2. 00:21:28.100 --> 00:21:29.850 So one case I measured the input and tried 00:21:29.850 --> 00:21:31.869 to predict the output to the tank. 00:21:31.869 --> 00:21:33.660 And in the other case I measured the output 00:21:33.660 --> 00:21:36.600 and tried to use it to predict the input. 00:21:36.600 --> 00:21:37.273 Yeah? 00:21:37.273 --> 00:21:39.125 AUDIENCE: So physically the same system 00:21:39.125 --> 00:21:40.790 with a choice by the operator. 00:21:40.790 --> 00:21:41.790 WILLIAM GREEN, JR.: Yes. 00:21:41.790 --> 00:21:42.569 That's right. 00:21:42.569 --> 00:21:43.110 That's right. 00:21:43.110 --> 00:21:45.270 So I formulated the model differently. 00:21:45.270 --> 00:21:49.080 The same underlying system but I formulated 00:21:49.080 --> 00:21:50.700 the model differently. 00:21:50.700 --> 00:21:52.950 And it led to completely different behavior 00:21:52.950 --> 00:21:55.770 in the underlying solution to the differential algebraic 00:21:55.770 --> 00:21:57.520 equations. 00:21:57.520 --> 00:21:59.160 So you've already solved problems 00:21:59.160 --> 00:22:00.360 that involve conservation. 00:22:00.360 --> 00:22:04.770 You've done things like ODE-IVPs on batch reactors, 00:22:04.770 --> 00:22:07.320 where the total amount of material in the reactor 00:22:07.320 --> 00:22:09.015 had to remain constant. 00:22:09.015 --> 00:22:10.890 Instead, you probably solved for the dynamics 00:22:10.890 --> 00:22:14.130 of each of the components undergoing the reaction. 00:22:14.130 --> 00:22:18.540 You might have checked to see if at each point in time 00:22:18.540 --> 00:22:21.200 the total mass in the reactor stayed the same. 00:22:21.200 --> 00:22:23.754 Because you solved all these differential equations, 00:22:23.754 --> 00:22:25.420 they were interconnected with each other 00:22:25.420 --> 00:22:27.450 but they all incurred some numerical error. 00:22:27.450 --> 00:22:29.370 It's possible that numerical error 00:22:29.370 --> 00:22:33.150 accrues in a way that actually has you losing mass 00:22:33.150 --> 00:22:35.430 from your system. 00:22:35.430 --> 00:22:37.230 So there may be benefits to actually trying 00:22:37.230 --> 00:22:39.360 to solve these sorts of systems of equations 00:22:39.360 --> 00:22:42.150 with, say, one less differential equation 00:22:42.150 --> 00:22:45.390 for one of the components and instead an algebraic constraint 00:22:45.390 --> 00:22:48.330 governing the total mass or total number of moles 00:22:48.330 --> 00:22:49.050 in the system. 00:22:53.010 --> 00:22:55.320 So these pop up all the time. 00:22:55.320 --> 00:22:58.080 You see them in models of reaction networks, 00:22:58.080 --> 00:23:00.540 where we try to use a pseudo steady state approximation. 00:23:00.540 --> 00:23:03.960 We say some reactions go so fast that they equilibrate. 00:23:03.960 --> 00:23:06.120 So they're not governed by differential equations 00:23:06.120 --> 00:23:11.790 but by algebraic equations for equilibria. 00:23:11.790 --> 00:23:14.520 They can pop up in models of control, where we neglect 00:23:14.520 --> 00:23:16.680 the controller dynamics. 00:23:16.680 --> 00:23:19.440 So you say I'm going to try to control this process 00:23:19.440 --> 00:23:23.010 and I get to instantaneously turn on or off this valve 00:23:23.010 --> 00:23:24.930 in response to a measurement. 00:23:24.930 --> 00:23:28.472 Actually, controllers have inherent dynamics in them. 00:23:28.472 --> 00:23:30.180 But if I don't put those dynamics in then 00:23:30.180 --> 00:23:32.760 I may get an algebraic equation instead of a differential 00:23:32.760 --> 00:23:34.755 equation for how the control process occurs. 00:23:39.690 --> 00:23:44.250 There are different types of DAEs that we talk about. 00:23:44.250 --> 00:23:48.640 So there's one type that's called semi-explicit DAEs. 00:23:48.640 --> 00:23:50.980 And in these types of differential algebraic 00:23:50.980 --> 00:23:53.710 equations we can write them in the form 00:23:53.710 --> 00:24:00.310 M dx dt is equal to f of x and t, some initial condition. 00:24:00.310 --> 00:24:03.410 You say that this looks like an ODE-IVP. 00:24:03.410 --> 00:24:06.430 M is called the mass matrix. 00:24:06.430 --> 00:24:10.990 And it may or may not be the case that M is invertible. 00:24:10.990 --> 00:24:16.150 M may or may not be a full rank matrix. 00:24:16.150 --> 00:24:20.600 So let's look at that stirred tank example one. 00:24:20.600 --> 00:24:24.890 This was the underlying equation, f of x and dx dt 00:24:24.890 --> 00:24:25.790 and t. 00:24:25.790 --> 00:24:27.860 Can write it in semi-explicit form. 00:24:27.860 --> 00:24:35.590 If c has two components, c1 and c2, then dc 2 dt-- oops, 00:24:35.590 --> 00:24:36.640 this is a typo here. 00:24:36.640 --> 00:24:37.810 I really apologize. 00:24:37.810 --> 00:24:42.180 This should be a 0, and this element should be a 1. 00:24:42.180 --> 00:24:43.230 I'll correct that online. 00:24:43.230 --> 00:24:46.200 I apologize for that. 00:24:46.200 --> 00:24:48.709 So this should be this matrix multiplied 00:24:48.709 --> 00:24:50.250 with this differential should give me 00:24:50.250 --> 00:24:52.750 dc 2 dt, this term here. 00:24:52.750 --> 00:24:56.410 And a 0 for the second equation. 00:24:56.410 --> 00:25:02.450 And then I've got q over v multiplying c1 and minus c2. 00:25:02.450 --> 00:25:05.010 So that gives me the first line here. 00:25:05.010 --> 00:25:09.850 I've got minus c1 coming on the second equation here. 00:25:09.850 --> 00:25:12.820 And those balance with a 0 and a gamma. 00:25:12.820 --> 00:25:16.000 So the lower line reads like this equation here 00:25:16.000 --> 00:25:17.650 and the upper line-- 00:25:17.650 --> 00:25:19.660 fix this typo-- replace 0 with one, 00:25:19.660 --> 00:25:21.790 reads like the upper equation here. 00:25:21.790 --> 00:25:24.950 This is the semi-explicit form of the differential equation. 00:25:24.950 --> 00:25:27.970 You can see this mass matrix is singular and has 00:25:27.970 --> 00:25:29.200 a row that's all zeros. 00:25:29.200 --> 00:25:31.000 There's no way to invert this matrix 00:25:31.000 --> 00:25:34.540 and convert it into a system of ODE-IVPs, which you already 00:25:34.540 --> 00:25:37.360 know how to solve. 00:25:37.360 --> 00:25:41.080 If the mass matrix is full rank, which can happen, 00:25:41.080 --> 00:25:43.570 you can formulate your model in such a way 00:25:43.570 --> 00:25:45.980 that it takes on this semi-explicit form. 00:25:45.980 --> 00:25:48.530 But if the mass matrix is full rank and invertible, 00:25:48.530 --> 00:25:50.730 then you just invert the mass matrix. 00:25:50.730 --> 00:25:53.260 And now you can apply all your ODE-IVP techniques 00:25:53.260 --> 00:25:56.710 to solve this thing. 00:25:56.710 --> 00:25:58.960 The mass matrix doesn't need to be constant. 00:25:58.960 --> 00:26:01.720 Could depend on the concentration. 00:26:01.720 --> 00:26:06.397 It can depend on the states as well in this form. 00:26:06.397 --> 00:26:08.480 I just chose an example where that isn't the case. 00:26:08.480 --> 00:26:10.345 But in general, it can depend on the states. 00:26:10.345 --> 00:26:12.662 If it depends on the states and we invert it, 00:26:12.662 --> 00:26:14.620 then we need to be sure that the mass matrix is 00:26:14.620 --> 00:26:18.100 invertible for all values of the states, which may or may not 00:26:18.100 --> 00:26:19.395 be easy to ensure. 00:26:24.700 --> 00:26:28.010 So here's what I'd like you to try to do. 00:26:28.010 --> 00:26:30.850 So here's that stirred tank example two. 00:26:30.850 --> 00:26:34.240 The only difference between this and the previous one 00:26:34.240 --> 00:26:38.130 is that c2 now is set equal to gamma instead of c1. 00:26:38.130 --> 00:26:40.430 I want you to try to write it in semi-explicit form. 00:26:40.430 --> 00:26:42.790 A sort of test of understanding. 00:26:42.790 --> 00:26:46.700 Can you define the mass matrix? 00:26:46.700 --> 00:26:48.380 Can you write out the right-hand side 00:26:48.380 --> 00:26:50.175 of the equation in semi-explicit form? 00:26:50.175 --> 00:26:51.800 Take a second, work with your neighbor. 00:26:51.800 --> 00:26:53.508 See if you can figure out how to do that. 00:28:04.658 --> 00:28:06.020 AUDIENCE: Is this clear? 00:28:06.020 --> 00:28:07.296 OK to do this? 00:28:12.920 --> 00:28:15.420 WILLIAM GREEN, JR.: You can look back on the previous slide. 00:28:15.420 --> 00:28:18.690 Actually, the semi-exclusive form is going to be the same. 00:28:18.690 --> 00:28:20.910 This is a 0, that's a 1. 00:28:20.910 --> 00:28:23.850 The only difference is going to be c2 is now 00:28:23.850 --> 00:28:25.970 balanced with gamma, the input. 00:28:25.970 --> 00:28:27.880 So this minus 1 is got to shift over. 00:28:27.880 --> 00:28:29.460 So I switched the 0 for the minus 1, 00:28:29.460 --> 00:28:32.040 you got the second, the semi-explicit form. 00:28:38.270 --> 00:28:43.140 There's another way that these equations can pop up. 00:28:43.140 --> 00:28:44.880 So we have this mass matrix. 00:28:44.880 --> 00:28:47.790 And it might be the case that M is 00:28:47.790 --> 00:28:51.450 diagonal over many of the states and then 00:28:51.450 --> 00:28:52.789 0 for all the other states. 00:28:52.789 --> 00:28:55.080 So we might be able to write the equation in this form. 00:28:55.080 --> 00:28:56.204 This comes up all the time. 00:28:56.204 --> 00:28:58.260 So we might be able to write dx dt is 00:28:58.260 --> 00:29:01.110 a function of x, y, and t. 00:29:01.110 --> 00:29:06.950 And 0 is equal to another function of x, y, and t. 00:29:06.950 --> 00:29:09.300 And we have some initial conditions 00:29:09.300 --> 00:29:13.410 prescribed for the x's, of which we're 00:29:13.410 --> 00:29:16.560 taking the differential in one of these equations. 00:29:16.560 --> 00:29:19.890 And we've also got to satisfy those initial conditions 00:29:19.890 --> 00:29:23.400 for these variables y associated with the second nonlinear 00:29:23.400 --> 00:29:24.390 equation. 00:29:24.390 --> 00:29:27.750 The x's are called the differential states 00:29:27.750 --> 00:29:30.144 because they're involved in equations where they're 00:29:30.144 --> 00:29:32.310 difference with respect to time, or the differential 00:29:32.310 --> 00:29:34.260 with respect to time is taken. 00:29:34.260 --> 00:29:37.740 y are called the algebraic states 00:29:37.740 --> 00:29:41.280 because they only appear as themselves and not 00:29:41.280 --> 00:29:44.470 as their time differential in any of these equations. 00:29:44.470 --> 00:29:47.590 So we might want to solve for both the differential 00:29:47.590 --> 00:29:50.040 and the algebraic states simultaneously. 00:29:50.040 --> 00:29:52.710 We have to know these algebraic states to determine 00:29:52.710 --> 00:29:54.024 the differential states. 00:29:54.024 --> 00:29:56.190 We have to know the differential states to determine 00:29:56.190 --> 00:29:57.120 the algebraic states. 00:30:00.771 --> 00:30:01.740 AUDIENCE: Professor? 00:30:01.740 --> 00:30:02.739 WILLIAM GREEN, JR.: Yes. 00:30:02.739 --> 00:30:07.695 AUDIENCE: [INAUDIBLE] 00:30:07.695 --> 00:30:09.070 WILLIAM GREEN, JR.: That's right. 00:30:09.070 --> 00:30:13.370 So here the x that I have is not all of the states. 00:30:13.370 --> 00:30:14.690 It's not all of the unknowns. 00:30:14.690 --> 00:30:17.390 It's only the set of the unknowns of which I 00:30:17.390 --> 00:30:18.576 have to take a differential. 00:30:18.576 --> 00:30:19.700 That's right, that's right. 00:30:25.040 --> 00:30:29.030 So we've sort of tooled down the kinds of equations 00:30:29.030 --> 00:30:30.810 we might want to look at. 00:30:30.810 --> 00:30:33.860 We've gone from the generic differential algebraic equation 00:30:33.860 --> 00:30:38.390 to a semi-explicit form to the splitting, this differential 00:30:38.390 --> 00:30:39.320 algebraic splitting. 00:30:39.320 --> 00:30:43.400 Many problems can be formulated in this fashion. 00:30:43.400 --> 00:30:45.251 This one's the easiest one to think 00:30:45.251 --> 00:30:47.000 about so that's the one that we'll discuss 00:30:47.000 --> 00:30:48.530 for most of the lecture. 00:30:51.740 --> 00:30:53.720 So let's look at some examples. 00:30:53.720 --> 00:30:57.340 I think examples are interesting to think about. 00:30:57.340 --> 00:30:58.420 So here's one. 00:30:58.420 --> 00:31:02.140 Think about a mass spring system. 00:31:02.140 --> 00:31:04.360 There's an algebraic equation that 00:31:04.360 --> 00:31:07.270 governs a conserved quantity, the total energy 00:31:07.270 --> 00:31:09.370 of this mass spring system. 00:31:09.370 --> 00:31:13.060 So the kinetic energy, 1/2 m times the velocity of the mass 00:31:13.060 --> 00:31:14.200 squared. 00:31:14.200 --> 00:31:18.444 Plus the energy stored in the spring, 1/2 kx squared 00:31:18.444 --> 00:31:20.110 has got to be equal to the total energy. 00:31:20.110 --> 00:31:22.670 And the total energy can change in time. 00:31:22.670 --> 00:31:27.010 So we have a differential algebraic equation. 00:31:27.010 --> 00:31:29.860 We have a differential state because we're taking 00:31:29.860 --> 00:31:31.420 the time derivative of x. 00:31:31.420 --> 00:31:33.610 And we'd like to be able to determine 00:31:33.610 --> 00:31:36.040 x as a function of time. 00:31:36.040 --> 00:31:38.920 So f of x, x dot and t. 00:31:38.920 --> 00:31:42.250 That's this equation, 1/2 mv squared plus one half kx 00:31:42.250 --> 00:31:45.254 squared minus E is equal to 0. 00:31:45.254 --> 00:31:47.170 We'd like to solve this differential algebraic 00:31:47.170 --> 00:31:48.070 equation. 00:31:48.070 --> 00:31:49.300 It's well-posed. 00:31:49.300 --> 00:31:54.580 We have one equation for one unknown, one state. 00:31:54.580 --> 00:31:58.260 The equation has a solution, which 00:31:58.260 --> 00:32:00.520 is a times cosine omega t. 00:32:00.520 --> 00:32:02.020 It's not the only possible solution. 00:32:02.020 --> 00:32:04.290 It's a solution to this equation. 00:32:04.290 --> 00:32:06.660 We can substitute that solution into the equation 00:32:06.660 --> 00:32:10.560 and determine, for example what this oscillation frequency 00:32:10.560 --> 00:32:11.490 omega has to be. 00:32:11.490 --> 00:32:12.810 It's square root of k over m. 00:32:12.810 --> 00:32:14.250 Or what the amplitude has to be. 00:32:14.250 --> 00:32:17.130 It's the square root of E over k. 00:32:17.130 --> 00:32:20.496 So you've seen this from physics before. 00:32:20.496 --> 00:32:21.870 We've already solved differential 00:32:21.870 --> 00:32:24.030 algebraic equations. 00:32:24.030 --> 00:32:26.520 The thing you saw before though, was actually 00:32:26.520 --> 00:32:28.500 conservation of momentum. 00:32:28.500 --> 00:32:32.190 You converted this entirely to a second order 00:32:32.190 --> 00:32:34.380 differential equation in x and solved it. 00:32:34.380 --> 00:32:37.980 But in principle, we could have tried to find various solutions 00:32:37.980 --> 00:32:39.150 to this equation instead. 00:32:39.150 --> 00:32:41.802 We would have said, well, we've observed masses and springs 00:32:41.802 --> 00:32:43.260 and we know they oscillate in time. 00:32:43.260 --> 00:32:45.450 So let's propose some oscillatory solutions 00:32:45.450 --> 00:32:48.300 and see what values of the amplitude and omega we 00:32:48.300 --> 00:32:50.356 need to satisfy this equation. 00:32:50.356 --> 00:32:52.230 It would have been a perfectly acceptable way 00:32:52.230 --> 00:32:53.280 of solving this problem. 00:33:01.010 --> 00:33:03.940 The thing is, this mass spring problem, 00:33:03.940 --> 00:33:05.980 the differential algebraic equation 00:33:05.980 --> 00:33:09.880 contains nonlinearities with respect to the states. 00:33:09.880 --> 00:33:11.470 So it's nonlinear in the velocity 00:33:11.470 --> 00:33:12.940 and it's nonlinear in the position. 00:33:12.940 --> 00:33:15.865 Nonlinear equations in general are hard to solve. 00:33:15.865 --> 00:33:17.740 When you converted that to a momentum balance 00:33:17.740 --> 00:33:20.110 you got linear equations. 00:33:20.110 --> 00:33:23.590 You got the second derivative of position 00:33:23.590 --> 00:33:27.160 was equal to stiffness over mass times position. 00:33:27.160 --> 00:33:29.290 It was a linear differential equation to solve, 00:33:29.290 --> 00:33:30.300 which is easy to solve. 00:33:30.300 --> 00:33:32.950 So you've got higher order differentials, 00:33:32.950 --> 00:33:36.310 but you linearize the equation so it made it easier to solve, 00:33:36.310 --> 00:33:38.740 in a certain sense. 00:33:38.740 --> 00:33:42.790 But we'd like to sometimes understand what 00:33:42.790 --> 00:33:44.245 we can do with these equations. 00:33:44.245 --> 00:33:45.970 There are certain limiting cases where 00:33:45.970 --> 00:33:48.053 we can do interesting things with these equations. 00:33:48.053 --> 00:33:49.600 I'm going to show you one now. 00:33:49.600 --> 00:33:53.650 So in the case where the Jacobian of this 00:33:53.650 --> 00:33:58.021 function f with respect to the differential of the state 00:33:58.021 --> 00:33:58.520 variable. 00:33:58.520 --> 00:34:00.340 So take the derivative of f with respect 00:34:00.340 --> 00:34:04.080 to x dot and hold time in the state constant. 00:34:04.080 --> 00:34:06.830 So this is a matrix, one that's full rank. 00:34:06.830 --> 00:34:10.330 Then the DA be represented as an equivalent to ODE. 00:34:10.330 --> 00:34:12.190 Here's how you can show that. 00:34:12.190 --> 00:34:16.389 The total differential of f is this Jacobian with respect 00:34:16.389 --> 00:34:21.969 to x dot times the change x dot plus the Jacobian with respect 00:34:21.969 --> 00:34:24.760 to x times the change in x. 00:34:24.760 --> 00:34:28.400 Plus the partial derivative of f with respect to time, 00:34:28.400 --> 00:34:30.230 times the change in time. 00:34:30.230 --> 00:34:32.230 And this total differential has to be equal to 0 00:34:32.230 --> 00:34:37.080 because f is equal to 0 for all states in time. 00:34:37.080 --> 00:34:39.949 So let's divide by dt. 00:34:39.949 --> 00:34:43.050 We'd like to know how f changes with time. 00:34:43.050 --> 00:34:45.139 So total differential of f with time. 00:34:45.139 --> 00:34:48.300 We divide each of these little differentials by time. 00:34:48.300 --> 00:34:53.060 And then we solve for dx dot dt. 00:34:53.060 --> 00:34:55.060 That's going to involve moving this term 00:34:55.060 --> 00:34:56.469 to the other side of the equation 00:34:56.469 --> 00:34:58.340 and inverting this Jacobian. 00:34:58.340 --> 00:35:00.910 So if I can invert the Jacobian then 00:35:00.910 --> 00:35:05.125 I can convert my DAE system into a system of ODEs. 00:35:08.197 --> 00:35:09.530 This equals 0 shouldn't be here. 00:35:09.530 --> 00:35:10.050 I apologize. 00:35:10.050 --> 00:35:13.120 That's another typo. 00:35:13.120 --> 00:35:15.920 I move this term to the other side of the equation 00:35:15.920 --> 00:35:19.580 to where the 0 is and I solve. 00:35:19.580 --> 00:35:24.180 So I went from a first order equation, a first order 00:35:24.180 --> 00:35:27.330 equation in x to a second order equation in x. 00:35:27.330 --> 00:35:30.120 Two differentials of x with time instead of one differential 00:35:30.120 --> 00:35:30.900 of x with time. 00:35:30.900 --> 00:35:32.525 But I could convert it to an equivalent 00:35:32.525 --> 00:35:34.990 system of ordinary differential equations. 00:35:34.990 --> 00:35:36.836 Does that make sense? 00:35:36.836 --> 00:35:39.220 Do you see how that works? 00:35:39.220 --> 00:35:40.776 Any questions? 00:35:40.776 --> 00:35:44.860 Let's go back and look at an example. 00:35:44.860 --> 00:35:45.900 The mass spring system. 00:35:49.430 --> 00:35:53.810 So here's our ODE that we're trying to solve. 00:35:53.810 --> 00:35:56.360 I'm sorry, our DAE that we're trying to solve. 00:35:56.360 --> 00:36:00.440 Can you convert this to an equivalent ordinary 00:36:00.440 --> 00:36:04.030 differential equation using the approach I just described? 00:36:04.030 --> 00:36:06.900 Think you guys can figure that out? 00:36:06.900 --> 00:36:07.400 Try it. 00:36:07.400 --> 00:36:10.555 You can work together. 00:36:10.555 --> 00:37:07.000 [SIDE CONVERSATIONS] 00:37:07.000 --> 00:37:10.070 WILLIAM GREEN, JR.: Guys are either tired 00:37:10.070 --> 00:37:11.880 or friendships have been destroyed 00:37:11.880 --> 00:37:15.150 over the mid-term season. 00:37:15.150 --> 00:37:15.770 Maybe both. 00:37:38.757 --> 00:37:40.590 So let's compute the things we need to know. 00:37:40.590 --> 00:37:44.400 We needed to know the partials of this f with respect 00:37:44.400 --> 00:37:46.180 to the states and with respect to time. 00:37:46.180 --> 00:37:50.190 So partial f respect to x dot, holding all the other variables 00:37:50.190 --> 00:37:51.840 constant, is mx dot. 00:37:51.840 --> 00:37:55.120 Partial f with respect to x is kx. 00:37:55.120 --> 00:37:58.591 Partial f with respect to t is t. 00:37:58.591 --> 00:38:00.090 And I told you before that you could 00:38:00.090 --> 00:38:05.645 write this says dx dot dt is equal to df 00:38:05.645 --> 00:38:10.670 dx dot inverse multiplied by-- 00:38:10.670 --> 00:38:14.264 is equal to minus df dx dot inverse multiplied by df dx. 00:38:14.264 --> 00:38:16.430 This was on the previous slide, which just gives you 00:38:16.430 --> 00:38:19.240 conservation of momentum. 00:38:19.240 --> 00:38:25.390 This gives you the acceleration is equal to the force minus kx 00:38:25.390 --> 00:38:26.870 divided by the mass m. 00:38:31.300 --> 00:38:33.875 So we can solve either, the DAE or the ODE. 00:38:37.890 --> 00:38:40.220 Here's a more complicated example, 00:38:40.220 --> 00:38:43.320 but it works, well, it's more complicated. 00:38:43.320 --> 00:38:46.010 So a lot of times we do molecular dynamics simulations. 00:38:46.010 --> 00:38:48.130 We try to model molecules moving around. 00:38:48.130 --> 00:38:50.232 There we also want to conserve energy. 00:38:50.232 --> 00:38:52.190 There's usually a conserved quantity, something 00:38:52.190 --> 00:38:53.870 we're trying to hold constant. 00:38:53.870 --> 00:38:55.160 So suppose we try to do that. 00:38:55.160 --> 00:38:57.770 The total energy in this system of molecules, 00:38:57.770 --> 00:39:01.610 maybe it's 1/2 mass times the velocity squared. 00:39:01.610 --> 00:39:04.010 Now x is the position of all these molecules. 00:39:04.010 --> 00:39:05.990 Ad x dot is the velocity of all these molecules 00:39:05.990 --> 00:39:09.650 plus the potential energy, which is a function of f, 00:39:09.650 --> 00:39:12.510 a function of their positions. 00:39:12.510 --> 00:39:17.391 So the DAE representing this constraint is this. 00:39:17.391 --> 00:39:18.890 It's the same as for the mass spring 00:39:18.890 --> 00:39:20.931 system, where we've replaced the potential energy 00:39:20.931 --> 00:39:21.980 with a generic one. 00:39:21.980 --> 00:39:24.520 This actually isn't a well-posed DEA. 00:39:24.520 --> 00:39:28.430 We have one equation for, let's see, 00:39:28.430 --> 00:39:31.146 if we have n molecules in here we have three n states. 00:39:31.146 --> 00:39:32.770 We can move around in three dimensions. 00:39:32.770 --> 00:39:35.605 So this is a really ill-posed equation. 00:39:35.605 --> 00:39:37.730 If we try to apply the same procedure, then we say, 00:39:37.730 --> 00:39:40.040 well, this quantity still should be conserved. 00:39:40.040 --> 00:39:43.180 Over time it needs to be equal to zero. 00:39:43.180 --> 00:39:47.630 Then we can try to compute all these differentials. 00:39:47.630 --> 00:39:50.630 But we'll find out that partial f, partial x dot, 00:39:50.630 --> 00:39:53.630 the Jacobian of f with respect to x dot, 00:39:53.630 --> 00:39:56.320 that's just the momentum. 00:39:56.320 --> 00:39:57.605 That's not a matrix. 00:39:57.605 --> 00:40:00.120 That's not invertible. 00:40:00.120 --> 00:40:02.390 So we can't convert this equation 00:40:02.390 --> 00:40:04.600 to an equivalent system of ordinary differential 00:40:04.600 --> 00:40:05.100 equations. 00:40:05.100 --> 00:40:06.558 So it's just something pathological 00:40:06.558 --> 00:40:08.930 for a one-dimensional system. 00:40:08.930 --> 00:40:13.010 Instead, what we find is zero is equal to the velocity 00:40:13.010 --> 00:40:17.740 dotted with the mass times the acceleration plus the gradient 00:40:17.740 --> 00:40:19.040 in the potential. 00:40:19.040 --> 00:40:22.880 So minus the forces acting on the molecules. 00:40:22.880 --> 00:40:24.590 We know that if we were to integrate 00:40:24.590 --> 00:40:27.590 the equations of motion for the molecules exactly, 00:40:27.590 --> 00:40:30.710 the momentum balances for those molecules exactly, 00:40:30.710 --> 00:40:32.900 then this term here would always be equal to zero. 00:40:32.900 --> 00:40:37.090 And we'd satisfy conservation of energy. 00:40:37.090 --> 00:40:40.130 Of course, all solutions of ordinary differential equations 00:40:40.130 --> 00:40:42.170 require approximations. 00:40:42.170 --> 00:40:44.312 So as you move the molecules around, 00:40:44.312 --> 00:40:45.770 their velocities and positions will 00:40:45.770 --> 00:40:48.200 tend to drift away from the solution 00:40:48.200 --> 00:40:51.950 where this is exactly in balance and equal to zero. 00:40:51.950 --> 00:40:54.950 If you desire to do that in the molecular dynamic simulation 00:40:54.950 --> 00:40:56.630 then you better have a method that 00:40:56.630 --> 00:40:59.960 somehow respects this geometric property instead. 00:40:59.960 --> 00:41:02.950 That the velocities are orthogonal to m 00:41:02.950 --> 00:41:05.000 x dot plus grad v. 00:41:05.000 --> 00:41:08.180 And there is a set of ODE-IVP methods that do that. 00:41:08.180 --> 00:41:10.841 They're called symplectic integrators. 00:41:10.841 --> 00:41:12.590 And they integrate the equations of motion 00:41:12.590 --> 00:41:16.220 while exerting some control over the error in the total energy. 00:41:16.220 --> 00:41:19.040 So if you were just the take a ODE45 00:41:19.040 --> 00:41:21.080 and try to solve this system of equations, 00:41:21.080 --> 00:41:24.530 over long times you would find that the energy drifts away 00:41:24.530 --> 00:41:26.340 from the place where it started. 00:41:26.340 --> 00:41:28.670 Even though your positions are being integrated 00:41:28.670 --> 00:41:31.297 with reasonable accuracy, over time 00:41:31.297 --> 00:41:33.380 the energy will drift away from where you started. 00:41:33.380 --> 00:41:36.590 The system will heat up effectively or cool down. 00:41:36.590 --> 00:41:39.200 That's undesirable if you're trying to do modeling. 00:41:39.200 --> 00:41:41.630 So instead, one uses the symplectic integrators, 00:41:41.630 --> 00:41:45.810 which are designed to control the error and the energy 00:41:45.810 --> 00:41:48.336 in addition to integrate the equations of motion. 00:41:48.336 --> 00:41:50.210 So you can look those up and read about them. 00:41:50.210 --> 00:41:51.694 I think they're quite interesting. 00:41:51.694 --> 00:41:53.360 If you're going to do molecular modeling 00:41:53.360 --> 00:41:55.430 you'd like to understand better how they work. 00:41:55.430 --> 00:41:56.929 But actually, conservation of energy 00:41:56.929 --> 00:41:59.870 here doesn't give us a well-posed DAE. 00:41:59.870 --> 00:42:02.060 We actually need the momentum equations 00:42:02.060 --> 00:42:07.680 to formulate or to determine the trajectories of the molecules. 00:42:07.680 --> 00:42:12.120 Let's talk about their numerical simulation. 00:42:12.120 --> 00:42:14.540 So let's think about these semi-explicit DAEs, 00:42:14.540 --> 00:42:16.640 the ones that can be split between differential 00:42:16.640 --> 00:42:17.790 and algebraic states. 00:42:17.790 --> 00:42:21.740 It's kind of useful to look at this problem in particular. 00:42:21.740 --> 00:42:25.700 And you might say, well, let's do the simplest possible thing. 00:42:25.700 --> 00:42:29.000 Let's do the forward Euler approximation. 00:42:29.000 --> 00:42:31.880 So let's take our derivative here and represent it 00:42:31.880 --> 00:42:36.500 as the state evaluated at x plus delta t minus the state at t 00:42:36.500 --> 00:42:38.080 divided by delta t. 00:42:38.080 --> 00:42:40.070 Let's make the right-hand side of our equation 00:42:40.070 --> 00:42:43.340 be evaluated at time t. 00:42:43.340 --> 00:42:46.550 And let's try to step in time, do time marching forward 00:42:46.550 --> 00:42:49.550 to determine the trajectory of x and t. 00:42:49.550 --> 00:42:51.440 How do the states change? 00:42:51.440 --> 00:42:53.210 We see the first equation is an easy one 00:42:53.210 --> 00:42:58.870 to solve because presumably I know x of t and y of t. 00:42:58.870 --> 00:43:01.550 Well, do I know y of t? y of t has to satisfy 00:43:01.550 --> 00:43:02.720 this second equation. 00:43:02.720 --> 00:43:06.200 So first, if I know x of t I need to solve this nonlinear 00:43:06.200 --> 00:43:08.190 equation determine y of t. 00:43:08.190 --> 00:43:09.830 Then I know x of t and y of t and I 00:43:09.830 --> 00:43:14.390 can step forward in time to determine x of t plus delta t. 00:43:14.390 --> 00:43:17.720 So every iteration here with the simplest possible approximation 00:43:17.720 --> 00:43:20.030 for the differential requires solution of a system 00:43:20.030 --> 00:43:21.330 of nonlinear equations. 00:43:24.230 --> 00:43:27.540 We already saw before that a forward Euler approximation 00:43:27.540 --> 00:43:32.640 is not a great one to apply to solutions of IDE-IVPs 00:43:32.640 --> 00:43:35.660 because it has very limited stability properties. 00:43:35.660 --> 00:43:37.630 You need small time steps in order 00:43:37.630 --> 00:43:39.560 to stably integrate in time. 00:43:42.334 --> 00:43:43.750 So even though this seemed like it 00:43:43.750 --> 00:43:46.480 was easy to get a solution for the differential equation, 00:43:46.480 --> 00:43:51.620 we've always got to satisfy this nonlinear equation here. 00:43:51.620 --> 00:43:54.970 So inherently, simulation of DAEs are implicit. 00:43:54.970 --> 00:43:57.700 They require solution of nonlinear equations 00:43:57.700 --> 00:43:59.980 in order to advance in time. 00:43:59.980 --> 00:44:03.590 There's actually no point in using this weakly stable method 00:44:03.590 --> 00:44:04.990 for any very small time steps. 00:44:04.990 --> 00:44:08.200 It makes more sense to choose a naturally implicit method 00:44:08.200 --> 00:44:11.830 for the differential equation, where I can take big time steps 00:44:11.830 --> 00:44:16.376 and still have reasonable stability of the integration. 00:44:16.376 --> 00:44:17.250 Does that make sense? 00:44:21.040 --> 00:44:25.180 Consider now the fully implicit DAE. 00:44:25.180 --> 00:44:28.480 So f of x, x dot, and t is equal to 0 00:44:28.480 --> 00:44:30.565 and have some set of initial conditions. x of 0 00:44:30.565 --> 00:44:32.590 is equal to x 0. 00:44:32.590 --> 00:44:34.960 Here there is no way in general to avoid solving 00:44:34.960 --> 00:44:37.330 systems of nonlinear equations. 00:44:37.330 --> 00:44:39.850 And so one substitutes often backwards 00:44:39.850 --> 00:44:42.670 difference approximations for x dot. 00:44:42.670 --> 00:44:45.355 And then you solve for the next point in time. 00:44:45.355 --> 00:44:47.230 So here's an example with the backward Euler, 00:44:47.230 --> 00:44:48.730 implicit Euler approximations. 00:44:48.730 --> 00:44:51.070 So I replaced the time derivative 00:44:51.070 --> 00:44:55.000 with the difference between the state at point tk 00:44:55.000 --> 00:44:57.160 and the state of point tk minus 1 00:44:57.160 --> 00:45:00.880 divided by the difference between tk and tk minus 1. 00:45:00.880 --> 00:45:03.670 The other state variables are evaluated at tk. 00:45:03.670 --> 00:45:05.130 So x at tk. 00:45:05.130 --> 00:45:06.430 And I evaluate time at tk. 00:45:06.430 --> 00:45:08.650 And I try to satisfy this equation 00:45:08.650 --> 00:45:10.460 and use it to determine x of tk. 00:45:10.460 --> 00:45:13.690 So I solve this system of nonlinear equations 00:45:13.690 --> 00:45:15.160 for x at tk. 00:45:15.160 --> 00:45:18.580 And then I go to the next point in time, tk plus 1. 00:45:18.580 --> 00:45:22.440 And I solve the same equation, but replacing tk with tk plus 1 00:45:22.440 --> 00:45:25.060 and tk minus 1 with tk. 00:45:25.060 --> 00:45:28.060 So at each iteration I solve a system of nonlinear equations. 00:45:28.060 --> 00:45:30.400 No more painful than doing the forward Euler 00:45:30.400 --> 00:45:32.860 approximation I showed you before. 00:45:32.860 --> 00:45:35.590 So that's good. 00:45:35.590 --> 00:45:38.865 But still a lot of work. 00:45:38.865 --> 00:45:40.240 We might wonder how do we control 00:45:40.240 --> 00:45:42.872 the error in these sorts of approximations. 00:45:42.872 --> 00:45:44.830 And that's what we're going to talk about next. 00:45:48.394 --> 00:45:50.310 So here's the equivalent time marching formula 00:45:50.310 --> 00:45:52.770 applied to a semi-explicit DAE. 00:45:55.910 --> 00:46:02.790 So I've got to satisfy the equation 00:46:02.790 --> 00:46:05.850 0 is equal to the derivative of x with respect to time, 00:46:05.850 --> 00:46:09.060 which I approximate with the backward Euler approximation, 00:46:09.060 --> 00:46:12.300 minus f evaluated at tk. 00:46:12.300 --> 00:46:19.410 And the equation g of x at tk, y at tk for this x of tk. 00:46:19.410 --> 00:46:20.670 And actually y of tk as well. 00:46:20.670 --> 00:46:22.800 I'm sorry, I left that off. y of tk as well. 00:46:22.800 --> 00:46:24.890 Got to solve for the differential 00:46:24.890 --> 00:46:28.340 and the algebraic states simultaneously. 00:46:28.340 --> 00:46:30.030 And I use a backwards difference formula 00:46:30.030 --> 00:46:31.877 to try to do this in a stable way. 00:46:36.850 --> 00:46:43.180 One way to think about in particular these sorts of DAEs 00:46:43.180 --> 00:46:46.990 is as exceptionally stiff ordinary differential 00:46:46.990 --> 00:46:47.490 equations. 00:46:47.490 --> 00:46:50.320 Did you guys discuss stiffness? 00:46:50.320 --> 00:46:55.750 So imagine if I replace the 0 with, say, dy dt. 00:46:55.750 --> 00:46:58.090 And I introduced a constant on the right-hand side, 00:46:58.090 --> 00:47:05.485 a multiplicative constant, say, 1 over lambda here. 00:47:05.485 --> 00:47:06.610 Let's say lambda, actually. 00:47:06.610 --> 00:47:08.610 Say I introduce a multiplicative constant lambda 00:47:08.610 --> 00:47:10.940 on the right-hand side here. 00:47:10.940 --> 00:47:15.010 So as lambda becomes very large, the system 00:47:15.010 --> 00:47:17.620 becomes increasingly stiff. 00:47:17.620 --> 00:47:19.600 The dynamics of this equation take 00:47:19.600 --> 00:47:22.790 place over very, very short time scales. 00:47:22.790 --> 00:47:24.820 So short in fact, that eventually, 00:47:24.820 --> 00:47:26.500 if I let lambda diverge to infinity, 00:47:26.500 --> 00:47:29.500 I've just got to satisfy this algebraic equation 00:47:29.500 --> 00:47:31.400 at each point in time. 00:47:31.400 --> 00:47:34.240 So these are exceptionally stiff equations. 00:47:34.240 --> 00:47:36.270 And the methods for solving these equations 00:47:36.270 --> 00:47:39.880 share a lot in common with stiff ODE-IVP solvers. 00:47:45.069 --> 00:47:46.610 So I've got a couple of minutes left. 00:47:46.610 --> 00:47:48.880 I'm going to work through some examples. 00:47:48.880 --> 00:47:51.540 So consider the stirred tank example one. 00:47:51.540 --> 00:47:53.250 dc 2 dt. 00:47:53.250 --> 00:47:57.210 It's related to q over v, flow rate over volume. 00:47:57.210 --> 00:47:58.770 The difference in the concentrations 00:47:58.770 --> 00:48:00.570 coming in and out of the tank. 00:48:00.570 --> 00:48:03.060 And c1 is equal to this input signal gamma. 00:48:03.060 --> 00:48:06.030 Maybe that's a measurement and I'm trying to predict c2. 00:48:06.030 --> 00:48:08.850 So can you apply the backward Euler method to these equations 00:48:08.850 --> 00:48:12.810 to determine c1 and c2? 00:48:12.810 --> 00:48:15.670 What would one step of the backward Euler method look 00:48:15.670 --> 00:48:24.500 like, say, going from time 0 to time delta t, or time t to t 00:48:24.500 --> 00:48:27.260 plus delta t, or tk to tk plus 1, however you 00:48:27.260 --> 00:48:28.748 want to write it? 00:49:22.314 --> 00:49:23.730 So do you get something like this? 00:49:26.480 --> 00:49:28.530 So this equation is easy. 00:49:28.530 --> 00:49:32.360 c1 to tk is gamma of tk. 00:49:32.360 --> 00:49:34.390 I got a substitute up here. 00:49:34.390 --> 00:49:38.570 c2 of tk minus c2 of tk minus 1 over delta t. 00:49:38.570 --> 00:49:41.210 T Just tk minus tk minus 1. 00:49:41.210 --> 00:49:42.800 And then solve for c2 of dk. 00:49:42.800 --> 00:49:44.390 And you get a formula like this. 00:49:49.010 --> 00:49:52.010 This formula was based on an approximation 00:49:52.010 --> 00:49:56.510 for the derivative, which had a leading order 00:49:56.510 --> 00:50:03.140 error that was proportional to tk minus tk minus 1. 00:50:03.140 --> 00:50:07.070 When I go through and I solve for c2 of tk, the leading order 00:50:07.070 --> 00:50:10.600 error in c2, the local truncation error 00:50:10.600 --> 00:50:16.250 is going to be order tk minus tk minus 1 squared. 00:50:16.250 --> 00:50:20.840 So as I shrink the spacing between each of my two time 00:50:20.840 --> 00:50:24.050 points, the local accuracy, the error 00:50:24.050 --> 00:50:27.830 induced in one advance of this time stepping algorithm 00:50:27.830 --> 00:50:29.790 is going to scale like the step size squared. 00:50:29.790 --> 00:50:31.550 So I have the error-- 00:50:31.550 --> 00:50:32.369 I have the spacing. 00:50:32.369 --> 00:50:34.160 The error will go down by a factor of four. 00:50:34.160 --> 00:50:37.667 I reduce the spacing by an order of magnitude. 00:50:37.667 --> 00:50:39.750 The error will go down by two orders of magnitude. 00:50:39.750 --> 00:50:43.605 This is a really desirable sort of error control and a method. 00:50:43.605 --> 00:50:44.480 And we can do better. 00:50:44.480 --> 00:50:48.460 You've seen us do better with other methods. 00:50:48.460 --> 00:50:51.180 Let's look at this equation now. 00:50:51.180 --> 00:50:52.590 This is stirred tank example two. 00:50:52.590 --> 00:50:53.797 I just replace c1 with c2. 00:50:53.797 --> 00:50:55.380 You don't have to write this one down. 00:50:55.380 --> 00:50:56.484 I'll just show it to you. 00:50:56.484 --> 00:50:58.650 The notes are online so you're not missing anything. 00:50:58.650 --> 00:51:00.900 All these things will pop up in the notes online. 00:51:00.900 --> 00:51:05.220 So c2 at tk is equal to gamma of tk. 00:51:05.220 --> 00:51:06.930 I got to come up to my first equation now 00:51:06.930 --> 00:51:08.700 and I have to substitute an approximation 00:51:08.700 --> 00:51:09.810 for the derivative. 00:51:09.810 --> 00:51:13.390 It's this, the backwards difference formula. 00:51:13.390 --> 00:51:17.610 So if find c1 of tk is c2 of tk plus v over q. 00:51:17.610 --> 00:51:21.660 c2 of tk minus c2 of tk minus 1 divided by the time step. 00:51:21.660 --> 00:51:23.640 Remember, this approximation for the derivative 00:51:23.640 --> 00:51:27.292 has a leading order error that goes like tk minus tk minus 1. 00:51:27.292 --> 00:51:32.235 So my approximation for c1, the local approximation for c1, 00:51:32.235 --> 00:51:34.040 is not second order accurate. 00:51:34.040 --> 00:51:38.010 It's only first order accurate. 00:51:38.010 --> 00:51:42.450 I applied the same approximation method to both equations. 00:51:42.450 --> 00:51:45.050 One equation I got something that was second order accurate. 00:51:45.050 --> 00:51:46.425 The next equation I got something 00:51:46.425 --> 00:51:47.840 that was first order accurate. 00:51:47.840 --> 00:51:49.520 They look almost identical. 00:51:49.520 --> 00:51:51.530 It might be hard to see from the start 00:51:51.530 --> 00:51:53.420 why it should work out that way. 00:51:53.420 --> 00:51:56.190 That's how it works out. 00:51:56.190 --> 00:51:57.410 Let's do one more example. 00:52:00.640 --> 00:52:03.680 Here's a system of DAEs. 00:52:03.680 --> 00:52:05.760 c2 dot is equal to c1. 00:52:05.760 --> 00:52:07.900 c3 dot is equal to c2. 00:52:07.900 --> 00:52:10.240 And 0 is equal to c3 minus gamma. 00:52:10.240 --> 00:52:12.250 So gamma, again, is some measurement. 00:52:12.250 --> 00:52:14.660 We solve this equation to determine c3. 00:52:14.660 --> 00:52:17.640 We substitute c3 up here to determine c2. 00:52:17.640 --> 00:52:20.600 We substitute c2 up here to determine c1. 00:52:20.600 --> 00:52:23.965 But we want to find the solution with the backward Euler method. 00:52:27.100 --> 00:52:32.680 And the solution there is going to be c3 of tk is gamma of tk. 00:52:32.680 --> 00:52:35.410 c2 is related to the derivative of c3. 00:52:35.410 --> 00:52:37.180 So I use my backwards difference formula 00:52:37.180 --> 00:52:39.060 for the derivative of c3. 00:52:39.060 --> 00:52:41.620 c1 is related to the derivative of c2. 00:52:41.620 --> 00:52:44.710 So I use my backward difference formula for c2. 00:52:44.710 --> 00:52:46.930 But how accurate are each of these approximations? 00:52:49.840 --> 00:52:52.780 So here I have to determine c2 from an approximation 00:52:52.780 --> 00:52:54.880 for the derivative of c3. 00:52:54.880 --> 00:52:56.920 We know this backwards difference approximation 00:52:56.920 --> 00:52:58.510 has a leading order error proportional 00:52:58.510 --> 00:53:02.350 to delta t, the time step. 00:53:02.350 --> 00:53:07.750 Now I got to approximate c1 with the finite difference backwards 00:53:07.750 --> 00:53:09.922 difference approximation in c2. 00:53:09.922 --> 00:53:12.380 We know this should incur an error proportional to the time 00:53:12.380 --> 00:53:17.350 step delta t, except I don't know c2 at tk exactly. 00:53:17.350 --> 00:53:23.530 Actually, c2 tk, the exact value of c2 at tk 00:53:23.530 --> 00:53:27.830 is this plus something proportional to the time step. 00:53:27.830 --> 00:53:30.560 Which means the exact value of c1 00:53:30.560 --> 00:53:34.430 is this plus something that's order one because I 00:53:34.430 --> 00:53:36.299 carry over the error from my solution 00:53:36.299 --> 00:53:37.340 to the previous equation. 00:53:37.340 --> 00:53:40.477 And I divide it by the time step. 00:53:40.477 --> 00:53:42.560 So I went from having something-- well, let's see. 00:53:42.560 --> 00:53:44.340 The first equation had order delta t 00:53:44.340 --> 00:53:47.040 squared local truncation error. 00:53:47.040 --> 00:53:50.010 As I shrink the spacing, my solution 00:53:50.010 --> 00:53:52.440 gets more and more locally accurate. 00:53:52.440 --> 00:53:57.300 And it seems to do it in a fairly robust fashion. 00:53:57.300 --> 00:53:59.910 The next system of equations we solved 00:53:59.910 --> 00:54:02.170 had a local truncation error proportional to delta t. 00:54:02.170 --> 00:54:03.545 Again, I can shrink delta t and I 00:54:03.545 --> 00:54:05.800 can make my solution more and more accurate. 00:54:05.800 --> 00:54:08.220 That seems like what you want in a numerical algorithm. 00:54:08.220 --> 00:54:11.580 This one, there's no hope. 00:54:11.580 --> 00:54:15.980 Change delta t, who knows what you're going to get? 00:54:15.980 --> 00:54:18.855 And this looks fairly innocuous. 00:54:18.855 --> 00:54:20.990 You just look at it and you say, well, OK. 00:54:20.990 --> 00:54:22.980 Let's solve it and see what happens. 00:54:22.980 --> 00:54:26.040 Well actually, you can't solve this problem accurately 00:54:26.040 --> 00:54:28.322 with the backward Euler method. 00:54:28.322 --> 00:54:30.030 There's something fundamentally different 00:54:30.030 --> 00:54:32.359 about this problem from the previous one. 00:54:32.359 --> 00:54:34.650 There's something fundamentally different about stirred 00:54:34.650 --> 00:54:37.500 tank example two from stirred tank example one. 00:54:37.500 --> 00:54:40.170 Right you might have already perceived that. 00:54:40.170 --> 00:54:42.160 You can see what's going on here. 00:54:42.160 --> 00:54:44.520 c3 is equal to gamma. 00:54:44.520 --> 00:54:46.025 That's the exact solution. 00:54:46.025 --> 00:54:48.390 c2 is equal to gamma dot. 00:54:48.390 --> 00:54:50.700 And c1 is equal to two derivatives 00:54:50.700 --> 00:54:53.860 of gamma, the exact solution to this equation. 00:54:53.860 --> 00:54:59.310 So c1 is incredibly sensitive to changes in gamma. 00:54:59.310 --> 00:55:00.620 And that's peculiar. 00:55:00.620 --> 00:55:03.180 And it has important numerical consequences. 00:55:03.180 --> 00:55:06.380 So I'm going to show you, our next lecture on Wednesday, 00:55:06.380 --> 00:55:09.620 how to formalize an understanding 00:55:09.620 --> 00:55:11.450 of the differences between these problems. 00:55:11.450 --> 00:55:14.990 And we'll talk about some more subtleties associated 00:55:14.990 --> 00:55:17.150 with differential algebraic equations. 00:55:17.150 --> 00:55:18.650 You guys can pick up your quizzes. 00:55:18.650 --> 00:55:23.000 We should do it outside rather than 00:55:23.000 --> 00:55:25.596 inside because the next class has to start.