1 00:00:01,580 --> 00:00:03,920 The following content is provided under a Creative 2 00:00:03,920 --> 00:00:05,340 Commons license. 3 00:00:05,340 --> 00:00:07,550 Your support will help MIT OpenCourseWare 4 00:00:07,550 --> 00:00:11,640 continue to offer high quality educational resources for free. 5 00:00:11,640 --> 00:00:14,180 To make a donation or to view additional materials 6 00:00:14,180 --> 00:00:18,110 from hundreds of MIT courses, visit MIT OpenCourseWare 7 00:00:18,110 --> 00:00:19,340 at ocw.mit.edu. 8 00:00:23,562 --> 00:00:25,270 WILLIAM GREEN, JR.: Shall we get started? 9 00:00:25,270 --> 00:00:27,340 Yes. 10 00:00:27,340 --> 00:00:29,710 Apologies for the notes being posted so late. 11 00:00:29,710 --> 00:00:33,020 I had trouble with the new two factor authentication system, 12 00:00:33,020 --> 00:00:35,980 which I need to get into Stellar and put things up there. 13 00:00:35,980 --> 00:00:39,610 My phone wasn't working so I had to find my iPad someplace 14 00:00:39,610 --> 00:00:42,542 and use that to get credentials to log in. 15 00:00:42,542 --> 00:00:44,500 So if you don't have them don't worry about it. 16 00:00:44,500 --> 00:00:47,540 The full notes are posted online, no missing pieces. 17 00:00:47,540 --> 00:00:49,480 So you can utilize those afterwards 18 00:00:49,480 --> 00:00:51,060 if something is unclear. 19 00:00:51,060 --> 00:00:51,560 All right? 20 00:00:56,250 --> 00:00:58,070 So your quizzes are graded. 21 00:00:58,070 --> 00:01:00,390 The TAs have them. 22 00:01:00,390 --> 00:01:02,610 I think overall, we were really pleased 23 00:01:02,610 --> 00:01:04,349 with the outcome of this quiz. 24 00:01:04,349 --> 00:01:06,360 It is very conceptual in nature and we 25 00:01:06,360 --> 00:01:08,790 thought people really understood a lot of the material 26 00:01:08,790 --> 00:01:09,540 quite well. 27 00:01:09,540 --> 00:01:12,440 And you can see that from the distribution of scores. 28 00:01:12,440 --> 00:01:15,650 So the mean on the quiz was about 77, 29 00:01:15,650 --> 00:01:17,250 a standard deviation of 12. 30 00:01:17,250 --> 00:01:19,930 I think everyone did really quite well. 31 00:01:19,930 --> 00:01:23,580 And on a number of the problems, perfect solutions were posted. 32 00:01:23,580 --> 00:01:26,670 So I'll take a second here. 33 00:01:26,670 --> 00:01:28,800 If there are any questions about the material that 34 00:01:28,800 --> 00:01:31,464 was on the quiz, or the format of it 35 00:01:31,464 --> 00:01:32,880 that you'd like to have addressed, 36 00:01:32,880 --> 00:01:35,351 I'm happy to answer them. 37 00:01:35,351 --> 00:01:35,850 No? 38 00:01:38,580 --> 00:01:41,440 Does it seem like it took too much time 39 00:01:41,440 --> 00:01:44,080 to complete all three questions in the two hours allotted? 40 00:01:44,080 --> 00:01:49,510 Or was it about the right length of material for the time slot? 41 00:01:49,510 --> 00:01:51,030 So we vetted it with the two TAs. 42 00:01:51,030 --> 00:01:53,720 And it took each of them less than an hour to complete it. 43 00:01:53,720 --> 00:01:55,969 So that was a good sign that we thought 44 00:01:55,969 --> 00:01:57,760 you guys should be able to get through most 45 00:01:57,760 --> 00:01:59,140 of the material in two hours. 46 00:02:05,560 --> 00:02:08,012 So let's recap. 47 00:02:08,012 --> 00:02:10,470 We're going to move on today to a topic called differential 48 00:02:10,470 --> 00:02:11,650 algebraic equations. 49 00:02:11,650 --> 00:02:14,310 But before doing that, there are a couple of concepts 50 00:02:14,310 --> 00:02:16,260 I want to review from numerical integration, 51 00:02:16,260 --> 00:02:19,669 or actually concepts that weren't covered when we talked 52 00:02:19,669 --> 00:02:21,210 about numerical integration on Friday 53 00:02:21,210 --> 00:02:22,480 that I think are important. 54 00:02:22,480 --> 00:02:23,855 And then I want to review briefly 55 00:02:23,855 --> 00:02:26,190 implicit methods for ODE-IVPs because those 56 00:02:26,190 --> 00:02:28,710 are going to be important for differential algebraic 57 00:02:28,710 --> 00:02:29,250 equations. 58 00:02:29,250 --> 00:02:32,550 So for numerical integration we talked 59 00:02:32,550 --> 00:02:35,250 about quadrature of various integrals, 60 00:02:35,250 --> 00:02:36,960 how to develop quadrature formulas. 61 00:02:36,960 --> 00:02:38,130 But there's actually one type of integral 62 00:02:38,130 --> 00:02:39,150 we didn't talk about there, which 63 00:02:39,150 --> 00:02:41,910 is sort of improper integrals, where the limits of integration 64 00:02:41,910 --> 00:02:43,440 are unbounded. 65 00:02:43,440 --> 00:02:46,670 And these come up all the time in engineering problems. 66 00:02:46,670 --> 00:02:49,560 We'd like to know what the, say, integrated 67 00:02:49,560 --> 00:02:53,576 response of some function is over all possible times. 68 00:02:53,576 --> 00:02:55,450 And we may have to evaluate that numerically. 69 00:02:55,450 --> 00:02:57,860 So how do you do those sorts of integrals? 70 00:02:57,860 --> 00:03:00,090 And it turns out the strategy that's most often used 71 00:03:00,090 --> 00:03:03,420 is to divide this domain up into two pieces. 72 00:03:03,420 --> 00:03:05,370 One piece, which is proper integral 73 00:03:05,370 --> 00:03:07,680 over a finite domain, and another piece, 74 00:03:07,680 --> 00:03:09,930 which is an improper integral over an infinite domain. 75 00:03:09,930 --> 00:03:11,513 This first integral, you can handle it 76 00:03:11,513 --> 00:03:15,120 with an ODE-IVP method or some higher order polynomial 77 00:03:15,120 --> 00:03:16,560 interpolation. 78 00:03:16,560 --> 00:03:19,494 But the second one we have to do differently. 79 00:03:19,494 --> 00:03:21,160 And there are couple of ways to do this. 80 00:03:21,160 --> 00:03:23,370 One is to transform variables. 81 00:03:23,370 --> 00:03:27,070 So you map this domain onto a finite domain. 82 00:03:27,070 --> 00:03:28,920 And the other option is to substitute 83 00:03:28,920 --> 00:03:30,420 an asymptotic approximation. 84 00:03:30,420 --> 00:03:34,740 So we say, when, say, tau is large, f of tau 85 00:03:34,740 --> 00:03:37,410 has a characteristic asymptotic approximation. 86 00:03:37,410 --> 00:03:40,200 Maybe it's exponentially decaying. 87 00:03:40,200 --> 00:03:43,080 And we take the integral of that asymptotic approximation. 88 00:03:43,080 --> 00:03:44,910 And the more accurate that approximation 89 00:03:44,910 --> 00:03:48,150 is the better our approximation for this whole 90 00:03:48,150 --> 00:03:49,800 integral will be. 91 00:03:49,800 --> 00:03:52,370 There's another type of interval that we have to do, 92 00:03:52,370 --> 00:03:54,960 where the same idea can be applied. 93 00:03:54,960 --> 00:03:58,485 That's the integral over integrable singularities. 94 00:03:58,485 --> 00:03:59,360 So here's an example. 95 00:04:02,810 --> 00:04:06,300 So say we want to integrate the function cosine of tau 96 00:04:06,300 --> 00:04:10,185 over square root of tau from 0 to some finite positive time 97 00:04:10,185 --> 00:04:12,800 point. 98 00:04:12,800 --> 00:04:15,260 As tau goes to 0, cosine goes to 1 and square root 99 00:04:15,260 --> 00:04:17,430 of tau diverges. 100 00:04:17,430 --> 00:04:20,390 But this integral actually has a finite value. 101 00:04:20,390 --> 00:04:21,980 This 1 over square root of tau, that's 102 00:04:21,980 --> 00:04:25,650 an integrable singularity. 103 00:04:25,650 --> 00:04:26,670 So how do you do this? 104 00:04:26,670 --> 00:04:29,730 Well, you want to split the domain again into two parts. 105 00:04:29,730 --> 00:04:32,840 One part that contains the integrable singularity 106 00:04:32,840 --> 00:04:35,270 and the other part, which excludes it. 107 00:04:35,270 --> 00:04:38,210 And in the part that contains the singularity 108 00:04:38,210 --> 00:04:41,540 we might do an asymptotic expansion of the function, 109 00:04:41,540 --> 00:04:44,510 and integrate each of those asymptotically accurate terms 110 00:04:44,510 --> 00:04:45,740 separately. 111 00:04:45,740 --> 00:04:48,530 So the integral of 1 over square root of tau 112 00:04:48,530 --> 00:04:51,680 is going to give us 2 tau 0 to the 1/2. 113 00:04:51,680 --> 00:04:54,230 The integral of this ratio here is 114 00:04:54,230 --> 00:04:57,840 going to give us minus 1/2 t0 to the 5/2. 115 00:04:57,840 --> 00:05:00,600 And then we have an integral over a finite domain, 116 00:05:00,600 --> 00:05:02,580 which doesn't include the singularity. 117 00:05:02,580 --> 00:05:05,390 And the accuracy of this method can be dictated by two things 118 00:05:05,390 --> 00:05:06,140 now. 119 00:05:06,140 --> 00:05:08,370 One is how accurately can we do this integral. 120 00:05:08,370 --> 00:05:13,190 And the other is how accurate is our asymptotic expansion here. 121 00:05:13,190 --> 00:05:16,130 So if we want higher accuracy we need more terms. 122 00:05:16,130 --> 00:05:17,810 We might need to be selective about how 123 00:05:17,810 --> 00:05:20,000 we choose the end point for this domain 124 00:05:20,000 --> 00:05:21,489 in order to minimize the error. 125 00:05:21,489 --> 00:05:23,030 There are lots of methods to do this. 126 00:05:23,030 --> 00:05:23,490 Sam. 127 00:05:23,490 --> 00:05:25,031 AUDIENCE: Will you talk a little more 128 00:05:25,031 --> 00:05:28,360 about what typical integral [INAUDIBLE] is? 129 00:05:28,360 --> 00:05:31,210 WILLIAM GREEN, JR.: Yeah. 130 00:05:31,210 --> 00:05:33,850 So you know these things. 131 00:05:33,850 --> 00:05:40,260 If I try to integrate 1 over x, as x goes to 0 132 00:05:40,260 --> 00:05:42,310 this thing will always diverge. 133 00:05:42,310 --> 00:05:46,180 You know the antiderivative of 1 over x is going to be the log. 134 00:05:46,180 --> 00:05:47,890 And the log diverges. 135 00:05:47,890 --> 00:05:53,950 So power is weaker than 1 over x, like 1 over x to the 1/2. 136 00:05:53,950 --> 00:05:55,290 Don't diverge like the log. 137 00:05:55,290 --> 00:05:56,830 They actually go to a finite value. 138 00:05:56,830 --> 00:05:59,770 The log is like the most weakly singular function there is. 139 00:05:59,770 --> 00:06:02,800 And integrals over weaker powers of 1 over x actually 140 00:06:02,800 --> 00:06:04,994 don't produce a singularity anymore. 141 00:06:04,994 --> 00:06:06,910 They're what we call integrable singularities. 142 00:06:06,910 --> 00:06:11,200 The function itself diverges but its integral does not. 143 00:06:11,200 --> 00:06:15,100 So you usually have this asymptotic power law type 144 00:06:15,100 --> 00:06:16,270 behavior. 145 00:06:16,270 --> 00:06:18,140 OK. 146 00:06:18,140 --> 00:06:20,100 So they come up all the time. 147 00:06:20,100 --> 00:06:23,910 We see these things in places like, say, 148 00:06:23,910 --> 00:06:27,460 squeezing a thin film of fluid between two plates. 149 00:06:27,460 --> 00:06:31,567 How long does it take for these two plates to come together? 150 00:06:31,567 --> 00:06:33,900 You might say, well, as the plates get closer and closer 151 00:06:33,900 --> 00:06:36,737 together the pressure starts to diverge in the gap 152 00:06:36,737 --> 00:06:38,070 and they'll never come together. 153 00:06:38,070 --> 00:06:40,950 But maybe depending on the geometry of the plates 154 00:06:40,950 --> 00:06:42,750 there can be a sort of finite time 155 00:06:42,750 --> 00:06:44,190 at which they'll come together. 156 00:06:44,190 --> 00:06:45,273 Depends on their geometry. 157 00:06:48,557 --> 00:06:49,390 So that's one topic. 158 00:06:49,390 --> 00:06:51,473 I think that's useful because these things come up 159 00:06:51,473 --> 00:06:52,030 all the time. 160 00:06:52,030 --> 00:06:55,330 And you can't actually handle those cases with quadrature 161 00:06:55,330 --> 00:06:56,800 by itself. 162 00:06:56,800 --> 00:06:57,910 The function diverges. 163 00:06:57,910 --> 00:06:59,680 There's no polynomial interpolation 164 00:06:59,680 --> 00:07:01,150 that's suitable to match it. 165 00:07:01,150 --> 00:07:02,777 So the quadrature can't handle it. 166 00:07:02,777 --> 00:07:04,360 And if you're trying to integrate over 167 00:07:04,360 --> 00:07:07,210 an infinite domain, well, good luck. 168 00:07:07,210 --> 00:07:09,580 You're never going to be able to fit enough points in 169 00:07:09,580 --> 00:07:14,440 to know that you've accurately integrated out to infinity. 170 00:07:14,440 --> 00:07:17,260 The other thing I want to recap here, 171 00:07:17,260 --> 00:07:19,750 and this is a problem for you to try. 172 00:07:19,750 --> 00:07:23,890 So here we have a simple first order ordinary differential 173 00:07:23,890 --> 00:07:24,390 equation. 174 00:07:24,390 --> 00:07:27,080 We want to use the implicit Euler method to solve it. 175 00:07:27,080 --> 00:07:29,840 So can you give a closed form solution 176 00:07:29,840 --> 00:07:33,260 for one step of the implicit Euler method, 177 00:07:33,260 --> 00:07:36,620 or for n steps of the implicit Euler method? 178 00:07:36,620 --> 00:07:38,300 What result does that produce? 179 00:07:38,300 --> 00:07:40,520 Can you compute this? 180 00:07:40,520 --> 00:07:48,556 Do you guys talk about backward Euler methods for for ODE-IVPs? 181 00:07:48,556 --> 00:07:51,180 Backwards difference methods in general and the backwards Euler 182 00:07:51,180 --> 00:07:54,600 method specifically? 183 00:07:54,600 --> 00:07:55,940 No. 184 00:07:55,940 --> 00:07:57,150 OK. 185 00:07:57,150 --> 00:08:01,770 Well, clearly I should have been here. 186 00:08:01,770 --> 00:08:03,150 So what's the strategy? 187 00:08:03,150 --> 00:08:05,100 We have to approximate the derivative. 188 00:08:05,100 --> 00:08:06,960 Yes? 189 00:08:06,960 --> 00:08:08,880 So the approximation for the derivative we use 190 00:08:08,880 --> 00:08:11,160 is a finite difference approximation. 191 00:08:11,160 --> 00:08:14,940 So we write dx dt as x at a point k 192 00:08:14,940 --> 00:08:18,870 plus 1 minus x at the point k divided by delta t. 193 00:08:18,870 --> 00:08:20,700 Backwards difference. 194 00:08:20,700 --> 00:08:24,930 We evaluate the right-hand side of the ODE-IVP at k 195 00:08:24,930 --> 00:08:27,940 plus 1, xk plus 1. 196 00:08:27,940 --> 00:08:30,360 And then we solve for xk plus 1. 197 00:08:30,360 --> 00:08:32,159 We just substitute in here. 198 00:08:32,159 --> 00:08:33,690 Backwards difference approximation 199 00:08:33,690 --> 00:08:35,320 for the differential. 200 00:08:35,320 --> 00:08:37,979 And xk plus 1 for the right-hand side. 201 00:08:37,979 --> 00:08:41,010 And we find that xk plus 1 will be 1 over 1 202 00:08:41,010 --> 00:08:46,950 minus delta t times this lambda times xk. 203 00:08:46,950 --> 00:08:51,680 And if we iterate this from k equals 0 up to some finite k, 204 00:08:51,680 --> 00:08:54,240 it's like multiplying by powers of 1 over 1 205 00:08:54,240 --> 00:08:55,980 minus delta t times lambda. 206 00:08:55,980 --> 00:08:57,044 Yes? 207 00:08:57,044 --> 00:08:57,960 Does this makes sense? 208 00:09:01,240 --> 00:09:05,910 So our backwards Euler solution to this fundamental problem 209 00:09:05,910 --> 00:09:08,386 is going to look like this. 210 00:09:08,386 --> 00:09:10,260 And a lot of times we ask about the stability 211 00:09:10,260 --> 00:09:12,240 of these sorts of solutions. 212 00:09:12,240 --> 00:09:16,850 So if for example, the solution to this equation 213 00:09:16,850 --> 00:09:20,000 was supposed to decay to zero, we 214 00:09:20,000 --> 00:09:21,770 would expect our numerical method 215 00:09:21,770 --> 00:09:23,890 to also yield results which decay to zero. 216 00:09:23,890 --> 00:09:27,110 We would call that stable. 217 00:09:27,110 --> 00:09:29,660 So we all know under what conditions 218 00:09:29,660 --> 00:09:34,662 does xk, for example, go to 0 as I change the product delta 219 00:09:34,662 --> 00:09:36,350 t times lambda. 220 00:09:36,350 --> 00:09:40,850 Lambda could be a real number or it can be a complex number too. 221 00:09:40,850 --> 00:09:42,360 Doesn't really matter. 222 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? 223 00:09:46,830 --> 00:09:49,340 Well what needs to be true is this quantity in parentheses 224 00:09:49,340 --> 00:09:50,960 needs to be smaller than 1 because I'm 225 00:09:50,960 --> 00:09:52,400 taking lots of products of things 226 00:09:52,400 --> 00:09:53,483 that are smaller than one. 227 00:09:53,483 --> 00:09:57,770 That will continue to shrink xk from x 0 until it goes to zero. 228 00:09:57,770 --> 00:10:01,790 If this is smaller than 1, then this quantity down here 229 00:10:01,790 --> 00:10:04,610 needs to be bigger than 1 in absolute value. 230 00:10:04,610 --> 00:10:07,010 And then I'll get solutions that slowly decay to zero. 231 00:10:07,010 --> 00:10:10,490 So if I ask, for what values of delta t times lambda 232 00:10:10,490 --> 00:10:13,770 is the absolute value of this thing bigger than 1, 233 00:10:13,770 --> 00:10:17,245 and I allow lambda to be a complex number, for example? 234 00:10:17,245 --> 00:10:20,110 I'll find out that's true when 1 minus delta t 235 00:10:20,110 --> 00:10:22,490 times the real part of lambda squared 236 00:10:22,490 --> 00:10:26,630 plus delta t times the imaginary part of lambda squared 237 00:10:26,630 --> 00:10:27,530 is bigger than 1. 238 00:10:30,100 --> 00:10:33,440 And if I plot the real part of delta t times 239 00:10:33,440 --> 00:10:37,630 lambda versus the imaginary part of delta times lambda, 240 00:10:37,630 --> 00:10:41,160 this inequality is represented by this pink zone 241 00:10:41,160 --> 00:10:42,470 on the outside of the circle. 242 00:10:42,470 --> 00:10:44,510 So any values of delta t times lambda 243 00:10:44,510 --> 00:10:49,460 that live in this pink area will give stable integration. 244 00:10:49,460 --> 00:10:54,500 And the solution xk, well, the k to zero as k goes to infinity. 245 00:10:57,020 --> 00:10:59,530 So it makes sense? 246 00:10:59,530 --> 00:11:02,010 This is a little funny though. 247 00:11:02,010 --> 00:11:03,780 What does the solution to this equation 248 00:11:03,780 --> 00:11:08,012 do when lambda is bigger than 1? 249 00:11:08,012 --> 00:11:09,470 When it's bigger than 0, excuse me. 250 00:11:09,470 --> 00:11:10,886 When lambda is bigger than 0, what 251 00:11:10,886 --> 00:11:14,860 does the solution to this ODE-IVP do? 252 00:11:14,860 --> 00:11:15,560 Does it decay? 253 00:11:15,560 --> 00:11:16,640 AUDIENCE: [INAUDIBLE] 254 00:11:16,640 --> 00:11:17,973 WILLIAM GREEN, JR.: It blows up. 255 00:11:17,973 --> 00:11:19,100 It grows exponentially. 256 00:11:19,100 --> 00:11:21,920 But the numerical method here, when 257 00:11:21,920 --> 00:11:25,070 I have lambda is bigger than zero, 258 00:11:25,070 --> 00:11:28,659 can be stable and produce decaying solutions. 259 00:11:28,659 --> 00:11:29,825 So this is sort of peculiar. 260 00:11:32,890 --> 00:11:35,290 So we're often concerned when we solve ODE-IVPs 261 00:11:35,290 --> 00:11:37,640 and we try to use these sorts of backwards difference 262 00:11:37,640 --> 00:11:40,880 methods, these implicit methods, to solve them. 263 00:11:40,880 --> 00:11:44,390 We're often concerned with getting stability. 264 00:11:44,390 --> 00:11:46,299 We want to get solutions that decay 265 00:11:46,299 --> 00:11:47,590 when they're supposed to decay. 266 00:11:47,590 --> 00:11:49,340 Sometimes we get solutions that decay even 267 00:11:49,340 --> 00:11:52,690 when the solution is supposed to blow up. 268 00:11:52,690 --> 00:11:55,240 So that's sort of funny. 269 00:11:55,240 --> 00:11:57,835 So the exact solution should look like this. 270 00:11:57,835 --> 00:11:59,710 But there are going to be circumstances where 271 00:11:59,710 --> 00:12:00,940 lambda is bigger than zero. 272 00:12:00,940 --> 00:12:02,530 In fact, lambda is bigger than-- 273 00:12:02,530 --> 00:12:06,340 lambda delta t is bigger than 2 along the real axis, 274 00:12:06,340 --> 00:12:11,060 where the solution will actually blow up rather than decay away. 275 00:12:11,060 --> 00:12:12,970 So stability and accuracy don't always 276 00:12:12,970 --> 00:12:15,760 correlate with each other. 277 00:12:15,760 --> 00:12:18,760 We really have to understand the underlying dynamics 278 00:12:18,760 --> 00:12:21,040 of the problem we're trying to solve 279 00:12:21,040 --> 00:12:24,110 before we choose a numerical method to apply to it. 280 00:12:24,110 --> 00:12:27,280 And this is going to be true of differential algebraic 281 00:12:27,280 --> 00:12:29,740 equations as well. 282 00:12:29,740 --> 00:12:33,310 Are there any questions about this? 283 00:12:33,310 --> 00:12:37,030 There's a discussion of these things in your textbook, 284 00:12:37,030 --> 00:12:39,160 these sorts of stability regions associated 285 00:12:39,160 --> 00:12:41,290 with different backwards difference methods. 286 00:12:41,290 --> 00:12:46,450 The implicit Euler method is one backwards difference formula. 287 00:12:46,450 --> 00:12:48,760 It's the simplest one. 288 00:12:48,760 --> 00:12:52,160 It's the least accurate one that you can derive 289 00:12:52,160 --> 00:12:53,420 but it's an easy one to use. 290 00:12:56,074 --> 00:12:58,490 But now we want to move on to a different sort of problem. 291 00:12:58,490 --> 00:13:01,737 And it's one that we come across in engineering all the time. 292 00:13:01,737 --> 00:13:03,320 These problems are called differential 293 00:13:03,320 --> 00:13:04,670 algebraic equations. 294 00:13:04,670 --> 00:13:08,270 And they're problems of the sort a function of some state 295 00:13:08,270 --> 00:13:12,500 variables, and the derivatives of those state variables 296 00:13:12,500 --> 00:13:15,147 and time is equal to zero. 297 00:13:15,147 --> 00:13:16,730 And there are some initial conditions. 298 00:13:16,730 --> 00:13:19,825 The state variables at that time 0 is equal to some x0. 299 00:13:24,840 --> 00:13:26,490 So this is a vector valued function 300 00:13:26,490 --> 00:13:29,320 of a vector valued state and it's time derivatives. 301 00:13:29,320 --> 00:13:31,650 And we want to understand the dynamics of this system. 302 00:13:31,650 --> 00:13:33,660 How does x vary with time in a way 303 00:13:33,660 --> 00:13:37,500 that's consistent with this governing equation? 304 00:13:37,500 --> 00:13:39,840 Usually we call these well-posed problems 305 00:13:39,840 --> 00:13:42,690 when the dimension of the state variable and the function 306 00:13:42,690 --> 00:13:43,590 are the same. 307 00:13:43,590 --> 00:13:46,730 They both live in the vector space rn. 308 00:13:46,730 --> 00:13:48,375 So we have n different states. 309 00:13:48,375 --> 00:13:50,550 I have n different functions that I have to satisfy. 310 00:13:50,550 --> 00:13:55,090 And I want to understand how x varies with time. 311 00:13:55,090 --> 00:13:57,320 Here's an example. 312 00:13:57,320 --> 00:13:58,200 A stirred tank. 313 00:13:58,200 --> 00:14:00,460 I call this stirred tank, example one. 314 00:14:00,460 --> 00:14:06,684 So into the stirred tank comes some concentration of solute. 315 00:14:06,684 --> 00:14:09,100 Out of the stirred tank comes some different concentration 316 00:14:09,100 --> 00:14:09,891 of the same solute. 317 00:14:09,891 --> 00:14:14,120 And we want to model the dynamics as a function of time. 318 00:14:14,120 --> 00:14:16,900 So we know that a material balance 319 00:14:16,900 --> 00:14:19,390 of the solute on the stirred tank, if it's well mixed, 320 00:14:19,390 --> 00:14:26,320 will tell us that d c2 dt is related to the volumetric flow 321 00:14:26,320 --> 00:14:28,960 rate into the tank divided by the volume of the tank 322 00:14:28,960 --> 00:14:32,140 multiplied by the difference between c1 and c2, 323 00:14:32,140 --> 00:14:37,112 the amount carried in and the amount carried out of the tank. 324 00:14:37,112 --> 00:14:38,890 And let's put a constraint in here too. 325 00:14:38,890 --> 00:14:42,430 Let's suppose we're trying to control this tank in some way. 326 00:14:42,430 --> 00:14:45,820 So we say that c1 of time has got 327 00:14:45,820 --> 00:14:49,190 to be equal to some control signal, gamma of t. 328 00:14:52,390 --> 00:14:55,810 And then we want to understand the dynamics of this system. 329 00:14:55,810 --> 00:14:59,710 At time 0, c2 is some c0, the initial concentration 330 00:14:59,710 --> 00:15:00,570 of the tank. 331 00:15:00,570 --> 00:15:03,490 And c1 is whatever the initial value of this control signal 332 00:15:03,490 --> 00:15:04,960 is. 333 00:15:04,960 --> 00:15:08,350 And the solution is c1 equals gamma of t. 334 00:15:08,350 --> 00:15:11,470 And c2 is, well, the initial concentration's 335 00:15:11,470 --> 00:15:13,090 got to decay out exponentially. 336 00:15:13,090 --> 00:15:15,220 And then I've got to integrate over how 337 00:15:15,220 --> 00:15:17,770 the input is changing in time. 338 00:15:17,770 --> 00:15:20,060 Those also produce exponential decay. 339 00:15:20,060 --> 00:15:22,870 So if the input is a little delta function spike 340 00:15:22,870 --> 00:15:25,330 then I'll get an exponential decay leaving the tank. 341 00:15:25,330 --> 00:15:27,310 If the input is changing dynamically in time, 342 00:15:27,310 --> 00:15:30,550 I need this convolution of the input with this integral. 343 00:15:30,550 --> 00:15:33,130 So you just to solve this system, 344 00:15:33,130 --> 00:15:34,900 this ordinary differential equation 345 00:15:34,900 --> 00:15:37,310 by substituting c1 in there. 346 00:15:37,310 --> 00:15:40,280 But this system of equations here is not just 347 00:15:40,280 --> 00:15:42,040 a differential equation. 348 00:15:42,040 --> 00:15:45,847 It's a differential equation and an algebraic equation. 349 00:15:45,847 --> 00:15:47,680 There's something peculiar about that that's 350 00:15:47,680 --> 00:15:50,140 different than just solving systems 351 00:15:50,140 --> 00:15:52,440 of ordinary differential equations. 352 00:15:52,440 --> 00:15:57,400 Oftentimes in models, we write down all 353 00:15:57,400 --> 00:15:59,720 the governing equations associated with the model. 354 00:15:59,720 --> 00:16:01,636 Some of them may be differential, some of them 355 00:16:01,636 --> 00:16:02,800 may be algebraic. 356 00:16:02,800 --> 00:16:04,960 There are times when we can look at those equations 357 00:16:04,960 --> 00:16:07,480 and see clever substitutions that we can make. 358 00:16:07,480 --> 00:16:10,840 But if you have a hundred equations, or a thousand 359 00:16:10,840 --> 00:16:12,900 equations, or a million equations, 360 00:16:12,900 --> 00:16:14,820 there's no way to do that reliably. 361 00:16:14,820 --> 00:16:18,310 We just kind of wind up with a system of ordinary differential 362 00:16:18,310 --> 00:16:21,730 equations and algebraic equations mixed together. 363 00:16:21,730 --> 00:16:24,777 And we have to solve these reliably. 364 00:16:24,777 --> 00:16:25,735 Here's another example. 365 00:16:28,650 --> 00:16:29,720 Oh, excuse me. 366 00:16:29,720 --> 00:16:31,800 So let me write this in this form. 367 00:16:31,800 --> 00:16:34,110 So we said there should be a vector valued 368 00:16:34,110 --> 00:16:38,460 function of the states and the derivatives of the states. 369 00:16:38,460 --> 00:16:41,490 The states are the concentration c1 and c2. 370 00:16:41,490 --> 00:16:45,090 And the derivative here is dc 2 dt. 371 00:16:45,090 --> 00:16:46,720 So here's my vector valued function. 372 00:16:46,720 --> 00:16:48,900 The first element is this, and it's equal to 0. 373 00:16:48,900 --> 00:16:51,660 The second element is this, and it's equal to 0. 374 00:16:51,660 --> 00:16:54,900 Here's my initial condition over the states. 375 00:16:54,900 --> 00:16:57,930 So I can just transform that list of equations 376 00:16:57,930 --> 00:17:00,450 to the sort of canonical form for a differential 377 00:17:00,450 --> 00:17:01,420 algebraic equation. 378 00:17:04,690 --> 00:17:07,619 Here's a separate example. 379 00:17:07,619 --> 00:17:12,790 Suppose instead, I'm trying to either control 380 00:17:12,790 --> 00:17:16,450 or make a measurement of the outlet concentrations c2. 381 00:17:16,450 --> 00:17:18,220 And we call that gamma of t. 382 00:17:22,060 --> 00:17:26,329 And I want to know now what is c2 and what is c1? 383 00:17:26,329 --> 00:17:29,340 So I've got to solve this system of equations for c2 and c1. 384 00:17:29,340 --> 00:17:31,740 The solution for c2 is easy. 385 00:17:31,740 --> 00:17:32,360 I know that. 386 00:17:32,360 --> 00:17:35,060 c2 is gamma, by definition. 387 00:17:35,060 --> 00:17:37,430 I got to substitute this into the first equation 388 00:17:37,430 --> 00:17:38,420 and solve for c1. 389 00:17:38,420 --> 00:17:43,199 And I see, c1 is gamma plus v over q gamma dot. 390 00:17:43,199 --> 00:17:45,240 And I have to have initial conditions to go along 391 00:17:45,240 --> 00:17:47,992 with this. 392 00:17:47,992 --> 00:17:49,200 There's something funny here. 393 00:17:49,200 --> 00:17:52,800 We can see the initial condition for c2 has to be gamma 0. 394 00:17:52,800 --> 00:17:57,710 The initial condition for c1, what does that have to be? 395 00:17:57,710 --> 00:18:01,430 Well, it needs to be consistent with the solution here. 396 00:18:01,430 --> 00:18:03,800 There were no free variables when 397 00:18:03,800 --> 00:18:05,484 I solved this equation for c1. 398 00:18:05,484 --> 00:18:07,400 It isn't like solving a differential equation, 399 00:18:07,400 --> 00:18:10,220 having a free coefficient to specify. 400 00:18:10,220 --> 00:18:13,130 So it better be the case that this c0 here 401 00:18:13,130 --> 00:18:18,680 is the same as gamma 0 plus v over q gamma dot 0. 402 00:18:18,680 --> 00:18:21,410 Somehow I have to know this input signal to prescribe 403 00:18:21,410 --> 00:18:23,030 the initial conditions for c1. 404 00:18:23,030 --> 00:18:28,805 So that's peculiar and different from just ODE-IVPs. 405 00:18:28,805 --> 00:18:29,680 You see this picture? 406 00:18:29,680 --> 00:18:32,940 You see how this goes together? 407 00:18:32,940 --> 00:18:34,630 The solution is funny too. 408 00:18:34,630 --> 00:18:36,880 Suppose we are trying to use this system of equations 409 00:18:36,880 --> 00:18:39,040 to do the following. 410 00:18:39,040 --> 00:18:43,720 We are trying to measure c2 and use the measurement of c2 411 00:18:43,720 --> 00:18:46,030 to predict c1. 412 00:18:46,030 --> 00:18:48,580 So gamma is the measurement, and we're 413 00:18:48,580 --> 00:18:52,060 trying to solve this system of differential algebraic 414 00:18:52,060 --> 00:18:55,120 equations to predict what c1 is. 415 00:18:55,120 --> 00:18:58,300 All measurements incur numerical error. 416 00:18:58,300 --> 00:19:01,560 So even though our signal for gamma that we measure 417 00:19:01,560 --> 00:19:06,670 may be continuous, it's bouncing around wildly. 418 00:19:06,670 --> 00:19:08,290 It's not a constant signal. 419 00:19:08,290 --> 00:19:10,900 It's going to move around a great deal. 420 00:19:10,900 --> 00:19:13,960 And c1, our predicted value for c1, 421 00:19:13,960 --> 00:19:15,670 depends on the derivatives of gamma, 422 00:19:15,670 --> 00:19:18,490 which means c1 is going to be incredibly sensitive to how 423 00:19:18,490 --> 00:19:21,440 that measurement is made. 424 00:19:21,440 --> 00:19:22,395 So it's peculiar. 425 00:19:22,395 --> 00:19:24,520 It means that there's going to be a lot of problems 426 00:19:24,520 --> 00:19:28,840 with stably solving these equations because the solutions 427 00:19:28,840 --> 00:19:34,390 can admit cases where they're not integrals of the input 428 00:19:34,390 --> 00:19:37,990 but they're actually related to derivatives of the input. 429 00:19:37,990 --> 00:19:39,960 Look at the solution back here again. 430 00:19:39,960 --> 00:19:43,140 Oop, sorry. 431 00:19:43,140 --> 00:19:45,960 The solution here was related to an integral of the input. 432 00:19:45,960 --> 00:19:48,030 Suppose I was making a measurement of c1 433 00:19:48,030 --> 00:19:50,940 in trying to predict c2 instead? 434 00:19:50,940 --> 00:19:53,480 My prediction for c2 would be related to the integral 435 00:19:53,480 --> 00:19:54,480 of that measured signal. 436 00:19:54,480 --> 00:19:56,850 And we know integrals smooth signals out. 437 00:19:56,850 --> 00:19:59,340 So it's not going to be so sensitive to the measurement. 438 00:19:59,340 --> 00:20:01,740 So one way of doing this seems really sensible. 439 00:20:01,740 --> 00:20:06,540 And the other way of doing it seems a little bizarre. 440 00:20:06,540 --> 00:20:09,317 Well, we can construct these equations however we want. 441 00:20:09,317 --> 00:20:11,400 So there's something about how we formulate models 442 00:20:11,400 --> 00:20:13,140 that are going to make them either easily 443 00:20:13,140 --> 00:20:18,850 solvable, as DAEs, or quite difficult to solve. 444 00:20:18,850 --> 00:20:23,190 So these pop up all the time in engineering problems. 445 00:20:23,190 --> 00:20:25,530 They pop up in dynamic simulations, 446 00:20:25,530 --> 00:20:27,900 where we have some sort of underlying conservation 447 00:20:27,900 --> 00:20:31,950 principle, or constraints, or equilibria. 448 00:20:31,950 --> 00:20:33,450 So if we have to conserve something, 449 00:20:33,450 --> 00:20:38,550 like total energy, total mass, total momentum, total particle 450 00:20:38,550 --> 00:20:41,760 number, total atomic species, total charge, 451 00:20:41,760 --> 00:20:44,040 there will always be an algebraic constraint that 452 00:20:44,040 --> 00:20:48,210 tells us the total mass has to stay this, 453 00:20:48,210 --> 00:20:51,500 or the total number of atoms of a certain type has to say this. 454 00:20:51,500 --> 00:20:52,230 Yeah, Ian. 455 00:20:52,230 --> 00:20:56,950 AUDIENCE: Just to be clear, were you saying, on the stirred tank 456 00:20:56,950 --> 00:20:59,722 example, that those were the same problem 457 00:20:59,722 --> 00:21:02,119 with different solution approaches? 458 00:21:02,119 --> 00:21:03,660 WILLIAM GREEN, JR.: So good question. 459 00:21:03,660 --> 00:21:04,285 Let me go back. 460 00:21:07,260 --> 00:21:10,960 So they are fundamentally different problems. 461 00:21:10,960 --> 00:21:17,450 So in one case I specified the value of c2 to be this gamma. 462 00:21:17,450 --> 00:21:21,450 I tried to make a measurement of c2 and then predict c1. 463 00:21:21,450 --> 00:21:24,180 In the other case, I specified c1. 464 00:21:24,180 --> 00:21:28,100 I tried to measure c1 and predict c2. 465 00:21:28,100 --> 00:21:29,850 So one case I measured the input and tried 466 00:21:29,850 --> 00:21:31,869 to predict the output to the tank. 467 00:21:31,869 --> 00:21:33,660 And in the other case I measured the output 468 00:21:33,660 --> 00:21:36,600 and tried to use it to predict the input. 469 00:21:36,600 --> 00:21:37,273 Yeah? 470 00:21:37,273 --> 00:21:39,125 AUDIENCE: So physically the same system 471 00:21:39,125 --> 00:21:40,790 with a choice by the operator. 472 00:21:40,790 --> 00:21:41,790 WILLIAM GREEN, JR.: Yes. 473 00:21:41,790 --> 00:21:42,569 That's right. 474 00:21:42,569 --> 00:21:43,110 That's right. 475 00:21:43,110 --> 00:21:45,270 So I formulated the model differently. 476 00:21:45,270 --> 00:21:49,080 The same underlying system but I formulated 477 00:21:49,080 --> 00:21:50,700 the model differently. 478 00:21:50,700 --> 00:21:52,950 And it led to completely different behavior 479 00:21:52,950 --> 00:21:55,770 in the underlying solution to the differential algebraic 480 00:21:55,770 --> 00:21:57,520 equations. 481 00:21:57,520 --> 00:21:59,160 So you've already solved problems 482 00:21:59,160 --> 00:22:00,360 that involve conservation. 483 00:22:00,360 --> 00:22:04,770 You've done things like ODE-IVPs on batch reactors, 484 00:22:04,770 --> 00:22:07,320 where the total amount of material in the reactor 485 00:22:07,320 --> 00:22:09,015 had to remain constant. 486 00:22:09,015 --> 00:22:10,890 Instead, you probably solved for the dynamics 487 00:22:10,890 --> 00:22:14,130 of each of the components undergoing the reaction. 488 00:22:14,130 --> 00:22:18,540 You might have checked to see if at each point in time 489 00:22:18,540 --> 00:22:21,200 the total mass in the reactor stayed the same. 490 00:22:21,200 --> 00:22:23,754 Because you solved all these differential equations, 491 00:22:23,754 --> 00:22:25,420 they were interconnected with each other 492 00:22:25,420 --> 00:22:27,450 but they all incurred some numerical error. 493 00:22:27,450 --> 00:22:29,370 It's possible that numerical error 494 00:22:29,370 --> 00:22:33,150 accrues in a way that actually has you losing mass 495 00:22:33,150 --> 00:22:35,430 from your system. 496 00:22:35,430 --> 00:22:37,230 So there may be benefits to actually trying 497 00:22:37,230 --> 00:22:39,360 to solve these sorts of systems of equations 498 00:22:39,360 --> 00:22:42,150 with, say, one less differential equation 499 00:22:42,150 --> 00:22:45,390 for one of the components and instead an algebraic constraint 500 00:22:45,390 --> 00:22:48,330 governing the total mass or total number of moles 501 00:22:48,330 --> 00:22:49,050 in the system. 502 00:22:53,010 --> 00:22:55,320 So these pop up all the time. 503 00:22:55,320 --> 00:22:58,080 You see them in models of reaction networks, 504 00:22:58,080 --> 00:23:00,540 where we try to use a pseudo steady state approximation. 505 00:23:00,540 --> 00:23:03,960 We say some reactions go so fast that they equilibrate. 506 00:23:03,960 --> 00:23:06,120 So they're not governed by differential equations 507 00:23:06,120 --> 00:23:11,790 but by algebraic equations for equilibria. 508 00:23:11,790 --> 00:23:14,520 They can pop up in models of control, where we neglect 509 00:23:14,520 --> 00:23:16,680 the controller dynamics. 510 00:23:16,680 --> 00:23:19,440 So you say I'm going to try to control this process 511 00:23:19,440 --> 00:23:23,010 and I get to instantaneously turn on or off this valve 512 00:23:23,010 --> 00:23:24,930 in response to a measurement. 513 00:23:24,930 --> 00:23:28,472 Actually, controllers have inherent dynamics in them. 514 00:23:28,472 --> 00:23:30,180 But if I don't put those dynamics in then 515 00:23:30,180 --> 00:23:32,760 I may get an algebraic equation instead of a differential 516 00:23:32,760 --> 00:23:34,755 equation for how the control process occurs. 517 00:23:39,690 --> 00:23:44,250 There are different types of DAEs that we talk about. 518 00:23:44,250 --> 00:23:48,640 So there's one type that's called semi-explicit DAEs. 519 00:23:48,640 --> 00:23:50,980 And in these types of differential algebraic 520 00:23:50,980 --> 00:23:53,710 equations we can write them in the form 521 00:23:53,710 --> 00:24:00,310 M dx dt is equal to f of x and t, some initial condition. 522 00:24:00,310 --> 00:24:03,410 You say that this looks like an ODE-IVP. 523 00:24:03,410 --> 00:24:06,430 M is called the mass matrix. 524 00:24:06,430 --> 00:24:10,990 And it may or may not be the case that M is invertible. 525 00:24:10,990 --> 00:24:16,150 M may or may not be a full rank matrix. 526 00:24:16,150 --> 00:24:20,600 So let's look at that stirred tank example one. 527 00:24:20,600 --> 00:24:24,890 This was the underlying equation, f of x and dx dt 528 00:24:24,890 --> 00:24:25,790 and t. 529 00:24:25,790 --> 00:24:27,860 Can write it in semi-explicit form. 530 00:24:27,860 --> 00:24:35,590 If c has two components, c1 and c2, then dc 2 dt-- oops, 531 00:24:35,590 --> 00:24:36,640 this is a typo here. 532 00:24:36,640 --> 00:24:37,810 I really apologize. 533 00:24:37,810 --> 00:24:42,180 This should be a 0, and this element should be a 1. 534 00:24:42,180 --> 00:24:43,230 I'll correct that online. 535 00:24:43,230 --> 00:24:46,200 I apologize for that. 536 00:24:46,200 --> 00:24:48,709 So this should be this matrix multiplied 537 00:24:48,709 --> 00:24:50,250 with this differential should give me 538 00:24:50,250 --> 00:24:52,750 dc 2 dt, this term here. 539 00:24:52,750 --> 00:24:56,410 And a 0 for the second equation. 540 00:24:56,410 --> 00:25:02,450 And then I've got q over v multiplying c1 and minus c2. 541 00:25:02,450 --> 00:25:05,010 So that gives me the first line here. 542 00:25:05,010 --> 00:25:09,850 I've got minus c1 coming on the second equation here. 543 00:25:09,850 --> 00:25:12,820 And those balance with a 0 and a gamma. 544 00:25:12,820 --> 00:25:16,000 So the lower line reads like this equation here 545 00:25:16,000 --> 00:25:17,650 and the upper line-- 546 00:25:17,650 --> 00:25:19,660 fix this typo-- replace 0 with one, 547 00:25:19,660 --> 00:25:21,790 reads like the upper equation here. 548 00:25:21,790 --> 00:25:24,950 This is the semi-explicit form of the differential equation. 549 00:25:24,950 --> 00:25:27,970 You can see this mass matrix is singular and has 550 00:25:27,970 --> 00:25:29,200 a row that's all zeros. 551 00:25:29,200 --> 00:25:31,000 There's no way to invert this matrix 552 00:25:31,000 --> 00:25:34,540 and convert it into a system of ODE-IVPs, which you already 553 00:25:34,540 --> 00:25:37,360 know how to solve. 554 00:25:37,360 --> 00:25:41,080 If the mass matrix is full rank, which can happen, 555 00:25:41,080 --> 00:25:43,570 you can formulate your model in such a way 556 00:25:43,570 --> 00:25:45,980 that it takes on this semi-explicit form. 557 00:25:45,980 --> 00:25:48,530 But if the mass matrix is full rank and invertible, 558 00:25:48,530 --> 00:25:50,730 then you just invert the mass matrix. 559 00:25:50,730 --> 00:25:53,260 And now you can apply all your ODE-IVP techniques 560 00:25:53,260 --> 00:25:56,710 to solve this thing. 561 00:25:56,710 --> 00:25:58,960 The mass matrix doesn't need to be constant. 562 00:25:58,960 --> 00:26:01,720 Could depend on the concentration. 563 00:26:01,720 --> 00:26:06,397 It can depend on the states as well in this form. 564 00:26:06,397 --> 00:26:08,480 I just chose an example where that isn't the case. 565 00:26:08,480 --> 00:26:10,345 But in general, it can depend on the states. 566 00:26:10,345 --> 00:26:12,662 If it depends on the states and we invert it, 567 00:26:12,662 --> 00:26:14,620 then we need to be sure that the mass matrix is 568 00:26:14,620 --> 00:26:18,100 invertible for all values of the states, which may or may not 569 00:26:18,100 --> 00:26:19,395 be easy to ensure. 570 00:26:24,700 --> 00:26:28,010 So here's what I'd like you to try to do. 571 00:26:28,010 --> 00:26:30,850 So here's that stirred tank example two. 572 00:26:30,850 --> 00:26:34,240 The only difference between this and the previous one 573 00:26:34,240 --> 00:26:38,130 is that c2 now is set equal to gamma instead of c1. 574 00:26:38,130 --> 00:26:40,430 I want you to try to write it in semi-explicit form. 575 00:26:40,430 --> 00:26:42,790 A sort of test of understanding. 576 00:26:42,790 --> 00:26:46,700 Can you define the mass matrix? 577 00:26:46,700 --> 00:26:48,380 Can you write out the right-hand side 578 00:26:48,380 --> 00:26:50,175 of the equation in semi-explicit form? 579 00:26:50,175 --> 00:26:51,800 Take a second, work with your neighbor. 580 00:26:51,800 --> 00:26:53,508 See if you can figure out how to do that. 581 00:28:04,658 --> 00:28:06,020 AUDIENCE: Is this clear? 582 00:28:06,020 --> 00:28:07,296 OK to do this? 583 00:28:12,920 --> 00:28:15,420 WILLIAM GREEN, JR.: You can look back on the previous slide. 584 00:28:15,420 --> 00:28:18,690 Actually, the semi-exclusive form is going to be the same. 585 00:28:18,690 --> 00:28:20,910 This is a 0, that's a 1. 586 00:28:20,910 --> 00:28:23,850 The only difference is going to be c2 is now 587 00:28:23,850 --> 00:28:25,970 balanced with gamma, the input. 588 00:28:25,970 --> 00:28:27,880 So this minus 1 is got to shift over. 589 00:28:27,880 --> 00:28:29,460 So I switched the 0 for the minus 1, 590 00:28:29,460 --> 00:28:32,040 you got the second, the semi-explicit form. 591 00:28:38,270 --> 00:28:43,140 There's another way that these equations can pop up. 592 00:28:43,140 --> 00:28:44,880 So we have this mass matrix. 593 00:28:44,880 --> 00:28:47,790 And it might be the case that M is 594 00:28:47,790 --> 00:28:51,450 diagonal over many of the states and then 595 00:28:51,450 --> 00:28:52,789 0 for all the other states. 596 00:28:52,789 --> 00:28:55,080 So we might be able to write the equation in this form. 597 00:28:55,080 --> 00:28:56,204 This comes up all the time. 598 00:28:56,204 --> 00:28:58,260 So we might be able to write dx dt is 599 00:28:58,260 --> 00:29:01,110 a function of x, y, and t. 600 00:29:01,110 --> 00:29:06,950 And 0 is equal to another function of x, y, and t. 601 00:29:06,950 --> 00:29:09,300 And we have some initial conditions 602 00:29:09,300 --> 00:29:13,410 prescribed for the x's, of which we're 603 00:29:13,410 --> 00:29:16,560 taking the differential in one of these equations. 604 00:29:16,560 --> 00:29:19,890 And we've also got to satisfy those initial conditions 605 00:29:19,890 --> 00:29:23,400 for these variables y associated with the second nonlinear 606 00:29:23,400 --> 00:29:24,390 equation. 607 00:29:24,390 --> 00:29:27,750 The x's are called the differential states 608 00:29:27,750 --> 00:29:30,144 because they're involved in equations where they're 609 00:29:30,144 --> 00:29:32,310 difference with respect to time, or the differential 610 00:29:32,310 --> 00:29:34,260 with respect to time is taken. 611 00:29:34,260 --> 00:29:37,740 y are called the algebraic states 612 00:29:37,740 --> 00:29:41,280 because they only appear as themselves and not 613 00:29:41,280 --> 00:29:44,470 as their time differential in any of these equations. 614 00:29:44,470 --> 00:29:47,590 So we might want to solve for both the differential 615 00:29:47,590 --> 00:29:50,040 and the algebraic states simultaneously. 616 00:29:50,040 --> 00:29:52,710 We have to know these algebraic states to determine 617 00:29:52,710 --> 00:29:54,024 the differential states. 618 00:29:54,024 --> 00:29:56,190 We have to know the differential states to determine 619 00:29:56,190 --> 00:29:57,120 the algebraic states. 620 00:30:00,771 --> 00:30:01,740 AUDIENCE: Professor? 621 00:30:01,740 --> 00:30:02,739 WILLIAM GREEN, JR.: Yes. 622 00:30:02,739 --> 00:30:07,695 AUDIENCE: [INAUDIBLE] 623 00:30:07,695 --> 00:30:09,070 WILLIAM GREEN, JR.: That's right. 624 00:30:09,070 --> 00:30:13,370 So here the x that I have is not all of the states. 625 00:30:13,370 --> 00:30:14,690 It's not all of the unknowns. 626 00:30:14,690 --> 00:30:17,390 It's only the set of the unknowns of which I 627 00:30:17,390 --> 00:30:18,576 have to take a differential. 628 00:30:18,576 --> 00:30:19,700 That's right, that's right. 629 00:30:25,040 --> 00:30:29,030 So we've sort of tooled down the kinds of equations 630 00:30:29,030 --> 00:30:30,810 we might want to look at. 631 00:30:30,810 --> 00:30:33,860 We've gone from the generic differential algebraic equation 632 00:30:33,860 --> 00:30:38,390 to a semi-explicit form to the splitting, this differential 633 00:30:38,390 --> 00:30:39,320 algebraic splitting. 634 00:30:39,320 --> 00:30:43,400 Many problems can be formulated in this fashion. 635 00:30:43,400 --> 00:30:45,251 This one's the easiest one to think 636 00:30:45,251 --> 00:30:47,000 about so that's the one that we'll discuss 637 00:30:47,000 --> 00:30:48,530 for most of the lecture. 638 00:30:51,740 --> 00:30:53,720 So let's look at some examples. 639 00:30:53,720 --> 00:30:57,340 I think examples are interesting to think about. 640 00:30:57,340 --> 00:30:58,420 So here's one. 641 00:30:58,420 --> 00:31:02,140 Think about a mass spring system. 642 00:31:02,140 --> 00:31:04,360 There's an algebraic equation that 643 00:31:04,360 --> 00:31:07,270 governs a conserved quantity, the total energy 644 00:31:07,270 --> 00:31:09,370 of this mass spring system. 645 00:31:09,370 --> 00:31:13,060 So the kinetic energy, 1/2 m times the velocity of the mass 646 00:31:13,060 --> 00:31:14,200 squared. 647 00:31:14,200 --> 00:31:18,444 Plus the energy stored in the spring, 1/2 kx squared 648 00:31:18,444 --> 00:31:20,110 has got to be equal to the total energy. 649 00:31:20,110 --> 00:31:22,670 And the total energy can change in time. 650 00:31:22,670 --> 00:31:27,010 So we have a differential algebraic equation. 651 00:31:27,010 --> 00:31:29,860 We have a differential state because we're taking 652 00:31:29,860 --> 00:31:31,420 the time derivative of x. 653 00:31:31,420 --> 00:31:33,610 And we'd like to be able to determine 654 00:31:33,610 --> 00:31:36,040 x as a function of time. 655 00:31:36,040 --> 00:31:38,920 So f of x, x dot and t. 656 00:31:38,920 --> 00:31:42,250 That's this equation, 1/2 mv squared plus one half kx 657 00:31:42,250 --> 00:31:45,254 squared minus E is equal to 0. 658 00:31:45,254 --> 00:31:47,170 We'd like to solve this differential algebraic 659 00:31:47,170 --> 00:31:48,070 equation. 660 00:31:48,070 --> 00:31:49,300 It's well-posed. 661 00:31:49,300 --> 00:31:54,580 We have one equation for one unknown, one state. 662 00:31:54,580 --> 00:31:58,260 The equation has a solution, which 663 00:31:58,260 --> 00:32:00,520 is a times cosine omega t. 664 00:32:00,520 --> 00:32:02,020 It's not the only possible solution. 665 00:32:02,020 --> 00:32:04,290 It's a solution to this equation. 666 00:32:04,290 --> 00:32:06,660 We can substitute that solution into the equation 667 00:32:06,660 --> 00:32:10,560 and determine, for example what this oscillation frequency 668 00:32:10,560 --> 00:32:11,490 omega has to be. 669 00:32:11,490 --> 00:32:12,810 It's square root of k over m. 670 00:32:12,810 --> 00:32:14,250 Or what the amplitude has to be. 671 00:32:14,250 --> 00:32:17,130 It's the square root of E over k. 672 00:32:17,130 --> 00:32:20,496 So you've seen this from physics before. 673 00:32:20,496 --> 00:32:21,870 We've already solved differential 674 00:32:21,870 --> 00:32:24,030 algebraic equations. 675 00:32:24,030 --> 00:32:26,520 The thing you saw before though, was actually 676 00:32:26,520 --> 00:32:28,500 conservation of momentum. 677 00:32:28,500 --> 00:32:32,190 You converted this entirely to a second order 678 00:32:32,190 --> 00:32:34,380 differential equation in x and solved it. 679 00:32:34,380 --> 00:32:37,980 But in principle, we could have tried to find various solutions 680 00:32:37,980 --> 00:32:39,150 to this equation instead. 681 00:32:39,150 --> 00:32:41,802 We would have said, well, we've observed masses and springs 682 00:32:41,802 --> 00:32:43,260 and we know they oscillate in time. 683 00:32:43,260 --> 00:32:45,450 So let's propose some oscillatory solutions 684 00:32:45,450 --> 00:32:48,300 and see what values of the amplitude and omega we 685 00:32:48,300 --> 00:32:50,356 need to satisfy this equation. 686 00:32:50,356 --> 00:32:52,230 It would have been a perfectly acceptable way 687 00:32:52,230 --> 00:32:53,280 of solving this problem. 688 00:33:01,010 --> 00:33:03,940 The thing is, this mass spring problem, 689 00:33:03,940 --> 00:33:05,980 the differential algebraic equation 690 00:33:05,980 --> 00:33:09,880 contains nonlinearities with respect to the states. 691 00:33:09,880 --> 00:33:11,470 So it's nonlinear in the velocity 692 00:33:11,470 --> 00:33:12,940 and it's nonlinear in the position. 693 00:33:12,940 --> 00:33:15,865 Nonlinear equations in general are hard to solve. 694 00:33:15,865 --> 00:33:17,740 When you converted that to a momentum balance 695 00:33:17,740 --> 00:33:20,110 you got linear equations. 696 00:33:20,110 --> 00:33:23,590 You got the second derivative of position 697 00:33:23,590 --> 00:33:27,160 was equal to stiffness over mass times position. 698 00:33:27,160 --> 00:33:29,290 It was a linear differential equation to solve, 699 00:33:29,290 --> 00:33:30,300 which is easy to solve. 700 00:33:30,300 --> 00:33:32,950 So you've got higher order differentials, 701 00:33:32,950 --> 00:33:36,310 but you linearize the equation so it made it easier to solve, 702 00:33:36,310 --> 00:33:38,740 in a certain sense. 703 00:33:38,740 --> 00:33:42,790 But we'd like to sometimes understand what 704 00:33:42,790 --> 00:33:44,245 we can do with these equations. 705 00:33:44,245 --> 00:33:45,970 There are certain limiting cases where 706 00:33:45,970 --> 00:33:48,053 we can do interesting things with these equations. 707 00:33:48,053 --> 00:33:49,600 I'm going to show you one now. 708 00:33:49,600 --> 00:33:53,650 So in the case where the Jacobian of this 709 00:33:53,650 --> 00:33:58,021 function f with respect to the differential of the state 710 00:33:58,021 --> 00:33:58,520 variable. 711 00:33:58,520 --> 00:34:00,340 So take the derivative of f with respect 712 00:34:00,340 --> 00:34:04,080 to x dot and hold time in the state constant. 713 00:34:04,080 --> 00:34:06,830 So this is a matrix, one that's full rank. 714 00:34:06,830 --> 00:34:10,330 Then the DA be represented as an equivalent to ODE. 715 00:34:10,330 --> 00:34:12,190 Here's how you can show that. 716 00:34:12,190 --> 00:34:16,389 The total differential of f is this Jacobian with respect 717 00:34:16,389 --> 00:34:21,969 to x dot times the change x dot plus the Jacobian with respect 718 00:34:21,969 --> 00:34:24,760 to x times the change in x. 719 00:34:24,760 --> 00:34:28,400 Plus the partial derivative of f with respect to time, 720 00:34:28,400 --> 00:34:30,230 times the change in time. 721 00:34:30,230 --> 00:34:32,230 And this total differential has to be equal to 0 722 00:34:32,230 --> 00:34:37,080 because f is equal to 0 for all states in time. 723 00:34:37,080 --> 00:34:39,949 So let's divide by dt. 724 00:34:39,949 --> 00:34:43,050 We'd like to know how f changes with time. 725 00:34:43,050 --> 00:34:45,139 So total differential of f with time. 726 00:34:45,139 --> 00:34:48,300 We divide each of these little differentials by time. 727 00:34:48,300 --> 00:34:53,060 And then we solve for dx dot dt. 728 00:34:53,060 --> 00:34:55,060 That's going to involve moving this term 729 00:34:55,060 --> 00:34:56,469 to the other side of the equation 730 00:34:56,469 --> 00:34:58,340 and inverting this Jacobian. 731 00:34:58,340 --> 00:35:00,910 So if I can invert the Jacobian then 732 00:35:00,910 --> 00:35:05,125 I can convert my DAE system into a system of ODEs. 733 00:35:08,197 --> 00:35:09,530 This equals 0 shouldn't be here. 734 00:35:09,530 --> 00:35:10,050 I apologize. 735 00:35:10,050 --> 00:35:13,120 That's another typo. 736 00:35:13,120 --> 00:35:15,920 I move this term to the other side of the equation 737 00:35:15,920 --> 00:35:19,580 to where the 0 is and I solve. 738 00:35:19,580 --> 00:35:24,180 So I went from a first order equation, a first order 739 00:35:24,180 --> 00:35:27,330 equation in x to a second order equation in x. 740 00:35:27,330 --> 00:35:30,120 Two differentials of x with time instead of one differential 741 00:35:30,120 --> 00:35:30,900 of x with time. 742 00:35:30,900 --> 00:35:32,525 But I could convert it to an equivalent 743 00:35:32,525 --> 00:35:34,990 system of ordinary differential equations. 744 00:35:34,990 --> 00:35:36,836 Does that make sense? 745 00:35:36,836 --> 00:35:39,220 Do you see how that works? 746 00:35:39,220 --> 00:35:40,776 Any questions? 747 00:35:40,776 --> 00:35:44,860 Let's go back and look at an example. 748 00:35:44,860 --> 00:35:45,900 The mass spring system. 749 00:35:49,430 --> 00:35:53,810 So here's our ODE that we're trying to solve. 750 00:35:53,810 --> 00:35:56,360 I'm sorry, our DAE that we're trying to solve. 751 00:35:56,360 --> 00:36:00,440 Can you convert this to an equivalent ordinary 752 00:36:00,440 --> 00:36:04,030 differential equation using the approach I just described? 753 00:36:04,030 --> 00:36:06,900 Think you guys can figure that out? 754 00:36:06,900 --> 00:36:07,400 Try it. 755 00:36:07,400 --> 00:36:10,555 You can work together. 756 00:36:10,555 --> 00:37:07,000 [SIDE CONVERSATIONS] 757 00:37:07,000 --> 00:37:10,070 WILLIAM GREEN, JR.: Guys are either tired 758 00:37:10,070 --> 00:37:11,880 or friendships have been destroyed 759 00:37:11,880 --> 00:37:15,150 over the mid-term season. 760 00:37:15,150 --> 00:37:15,770 Maybe both. 761 00:37:38,757 --> 00:37:40,590 So let's compute the things we need to know. 762 00:37:40,590 --> 00:37:44,400 We needed to know the partials of this f with respect 763 00:37:44,400 --> 00:37:46,180 to the states and with respect to time. 764 00:37:46,180 --> 00:37:50,190 So partial f respect to x dot, holding all the other variables 765 00:37:50,190 --> 00:37:51,840 constant, is mx dot. 766 00:37:51,840 --> 00:37:55,120 Partial f with respect to x is kx. 767 00:37:55,120 --> 00:37:58,591 Partial f with respect to t is t. 768 00:37:58,591 --> 00:38:00,090 And I told you before that you could 769 00:38:00,090 --> 00:38:05,645 write this says dx dot dt is equal to df 770 00:38:05,645 --> 00:38:10,670 dx dot inverse multiplied by-- 771 00:38:10,670 --> 00:38:14,264 is equal to minus df dx dot inverse multiplied by df dx. 772 00:38:14,264 --> 00:38:16,430 This was on the previous slide, which just gives you 773 00:38:16,430 --> 00:38:19,240 conservation of momentum. 774 00:38:19,240 --> 00:38:25,390 This gives you the acceleration is equal to the force minus kx 775 00:38:25,390 --> 00:38:26,870 divided by the mass m. 776 00:38:31,300 --> 00:38:33,875 So we can solve either, the DAE or the ODE. 777 00:38:37,890 --> 00:38:40,220 Here's a more complicated example, 778 00:38:40,220 --> 00:38:43,320 but it works, well, it's more complicated. 779 00:38:43,320 --> 00:38:46,010 So a lot of times we do molecular dynamics simulations. 780 00:38:46,010 --> 00:38:48,130 We try to model molecules moving around. 781 00:38:48,130 --> 00:38:50,232 There we also want to conserve energy. 782 00:38:50,232 --> 00:38:52,190 There's usually a conserved quantity, something 783 00:38:52,190 --> 00:38:53,870 we're trying to hold constant. 784 00:38:53,870 --> 00:38:55,160 So suppose we try to do that. 785 00:38:55,160 --> 00:38:57,770 The total energy in this system of molecules, 786 00:38:57,770 --> 00:39:01,610 maybe it's 1/2 mass times the velocity squared. 787 00:39:01,610 --> 00:39:04,010 Now x is the position of all these molecules. 788 00:39:04,010 --> 00:39:05,990 Ad x dot is the velocity of all these molecules 789 00:39:05,990 --> 00:39:09,650 plus the potential energy, which is a function of f, 790 00:39:09,650 --> 00:39:12,510 a function of their positions. 791 00:39:12,510 --> 00:39:17,391 So the DAE representing this constraint is this. 792 00:39:17,391 --> 00:39:18,890 It's the same as for the mass spring 793 00:39:18,890 --> 00:39:20,931 system, where we've replaced the potential energy 794 00:39:20,931 --> 00:39:21,980 with a generic one. 795 00:39:21,980 --> 00:39:24,520 This actually isn't a well-posed DEA. 796 00:39:24,520 --> 00:39:28,430 We have one equation for, let's see, 797 00:39:28,430 --> 00:39:31,146 if we have n molecules in here we have three n states. 798 00:39:31,146 --> 00:39:32,770 We can move around in three dimensions. 799 00:39:32,770 --> 00:39:35,605 So this is a really ill-posed equation. 800 00:39:35,605 --> 00:39:37,730 If we try to apply the same procedure, then we say, 801 00:39:37,730 --> 00:39:40,040 well, this quantity still should be conserved. 802 00:39:40,040 --> 00:39:43,180 Over time it needs to be equal to zero. 803 00:39:43,180 --> 00:39:47,630 Then we can try to compute all these differentials. 804 00:39:47,630 --> 00:39:50,630 But we'll find out that partial f, partial x dot, 805 00:39:50,630 --> 00:39:53,630 the Jacobian of f with respect to x dot, 806 00:39:53,630 --> 00:39:56,320 that's just the momentum. 807 00:39:56,320 --> 00:39:57,605 That's not a matrix. 808 00:39:57,605 --> 00:40:00,120 That's not invertible. 809 00:40:00,120 --> 00:40:02,390 So we can't convert this equation 810 00:40:02,390 --> 00:40:04,600 to an equivalent system of ordinary differential 811 00:40:04,600 --> 00:40:05,100 equations. 812 00:40:05,100 --> 00:40:06,558 So it's just something pathological 813 00:40:06,558 --> 00:40:08,930 for a one-dimensional system. 814 00:40:08,930 --> 00:40:13,010 Instead, what we find is zero is equal to the velocity 815 00:40:13,010 --> 00:40:17,740 dotted with the mass times the acceleration plus the gradient 816 00:40:17,740 --> 00:40:19,040 in the potential. 817 00:40:19,040 --> 00:40:22,880 So minus the forces acting on the molecules. 818 00:40:22,880 --> 00:40:24,590 We know that if we were to integrate 819 00:40:24,590 --> 00:40:27,590 the equations of motion for the molecules exactly, 820 00:40:27,590 --> 00:40:30,710 the momentum balances for those molecules exactly, 821 00:40:30,710 --> 00:40:32,900 then this term here would always be equal to zero. 822 00:40:32,900 --> 00:40:37,090 And we'd satisfy conservation of energy. 823 00:40:37,090 --> 00:40:40,130 Of course, all solutions of ordinary differential equations 824 00:40:40,130 --> 00:40:42,170 require approximations. 825 00:40:42,170 --> 00:40:44,312 So as you move the molecules around, 826 00:40:44,312 --> 00:40:45,770 their velocities and positions will 827 00:40:45,770 --> 00:40:48,200 tend to drift away from the solution 828 00:40:48,200 --> 00:40:51,950 where this is exactly in balance and equal to zero. 829 00:40:51,950 --> 00:40:54,950 If you desire to do that in the molecular dynamic simulation 830 00:40:54,950 --> 00:40:56,630 then you better have a method that 831 00:40:56,630 --> 00:40:59,960 somehow respects this geometric property instead. 832 00:40:59,960 --> 00:41:02,950 That the velocities are orthogonal to m 833 00:41:02,950 --> 00:41:05,000 x dot plus grad v. 834 00:41:05,000 --> 00:41:08,180 And there is a set of ODE-IVP methods that do that. 835 00:41:08,180 --> 00:41:10,841 They're called symplectic integrators. 836 00:41:10,841 --> 00:41:12,590 And they integrate the equations of motion 837 00:41:12,590 --> 00:41:16,220 while exerting some control over the error in the total energy. 838 00:41:16,220 --> 00:41:19,040 So if you were just the take a ODE45 839 00:41:19,040 --> 00:41:21,080 and try to solve this system of equations, 840 00:41:21,080 --> 00:41:24,530 over long times you would find that the energy drifts away 841 00:41:24,530 --> 00:41:26,340 from the place where it started. 842 00:41:26,340 --> 00:41:28,670 Even though your positions are being integrated 843 00:41:28,670 --> 00:41:31,297 with reasonable accuracy, over time 844 00:41:31,297 --> 00:41:33,380 the energy will drift away from where you started. 845 00:41:33,380 --> 00:41:36,590 The system will heat up effectively or cool down. 846 00:41:36,590 --> 00:41:39,200 That's undesirable if you're trying to do modeling. 847 00:41:39,200 --> 00:41:41,630 So instead, one uses the symplectic integrators, 848 00:41:41,630 --> 00:41:45,810 which are designed to control the error and the energy 849 00:41:45,810 --> 00:41:48,336 in addition to integrate the equations of motion. 850 00:41:48,336 --> 00:41:50,210 So you can look those up and read about them. 851 00:41:50,210 --> 00:41:51,694 I think they're quite interesting. 852 00:41:51,694 --> 00:41:53,360 If you're going to do molecular modeling 853 00:41:53,360 --> 00:41:55,430 you'd like to understand better how they work. 854 00:41:55,430 --> 00:41:56,929 But actually, conservation of energy 855 00:41:56,929 --> 00:41:59,870 here doesn't give us a well-posed DAE. 856 00:41:59,870 --> 00:42:02,060 We actually need the momentum equations 857 00:42:02,060 --> 00:42:07,680 to formulate or to determine the trajectories of the molecules. 858 00:42:07,680 --> 00:42:12,120 Let's talk about their numerical simulation. 859 00:42:12,120 --> 00:42:14,540 So let's think about these semi-explicit DAEs, 860 00:42:14,540 --> 00:42:16,640 the ones that can be split between differential 861 00:42:16,640 --> 00:42:17,790 and algebraic states. 862 00:42:17,790 --> 00:42:21,740 It's kind of useful to look at this problem in particular. 863 00:42:21,740 --> 00:42:25,700 And you might say, well, let's do the simplest possible thing. 864 00:42:25,700 --> 00:42:29,000 Let's do the forward Euler approximation. 865 00:42:29,000 --> 00:42:31,880 So let's take our derivative here and represent it 866 00:42:31,880 --> 00:42:36,500 as the state evaluated at x plus delta t minus the state at t 867 00:42:36,500 --> 00:42:38,080 divided by delta t. 868 00:42:38,080 --> 00:42:40,070 Let's make the right-hand side of our equation 869 00:42:40,070 --> 00:42:43,340 be evaluated at time t. 870 00:42:43,340 --> 00:42:46,550 And let's try to step in time, do time marching forward 871 00:42:46,550 --> 00:42:49,550 to determine the trajectory of x and t. 872 00:42:49,550 --> 00:42:51,440 How do the states change? 873 00:42:51,440 --> 00:42:53,210 We see the first equation is an easy one 874 00:42:53,210 --> 00:42:58,870 to solve because presumably I know x of t and y of t. 875 00:42:58,870 --> 00:43:01,550 Well, do I know y of t? y of t has to satisfy 876 00:43:01,550 --> 00:43:02,720 this second equation. 877 00:43:02,720 --> 00:43:06,200 So first, if I know x of t I need to solve this nonlinear 878 00:43:06,200 --> 00:43:08,190 equation determine y of t. 879 00:43:08,190 --> 00:43:09,830 Then I know x of t and y of t and I 880 00:43:09,830 --> 00:43:14,390 can step forward in time to determine x of t plus delta t. 881 00:43:14,390 --> 00:43:17,720 So every iteration here with the simplest possible approximation 882 00:43:17,720 --> 00:43:20,030 for the differential requires solution of a system 883 00:43:20,030 --> 00:43:21,330 of nonlinear equations. 884 00:43:24,230 --> 00:43:27,540 We already saw before that a forward Euler approximation 885 00:43:27,540 --> 00:43:32,640 is not a great one to apply to solutions of IDE-IVPs 886 00:43:32,640 --> 00:43:35,660 because it has very limited stability properties. 887 00:43:35,660 --> 00:43:37,630 You need small time steps in order 888 00:43:37,630 --> 00:43:39,560 to stably integrate in time. 889 00:43:42,334 --> 00:43:43,750 So even though this seemed like it 890 00:43:43,750 --> 00:43:46,480 was easy to get a solution for the differential equation, 891 00:43:46,480 --> 00:43:51,620 we've always got to satisfy this nonlinear equation here. 892 00:43:51,620 --> 00:43:54,970 So inherently, simulation of DAEs are implicit. 893 00:43:54,970 --> 00:43:57,700 They require solution of nonlinear equations 894 00:43:57,700 --> 00:43:59,980 in order to advance in time. 895 00:43:59,980 --> 00:44:03,590 There's actually no point in using this weakly stable method 896 00:44:03,590 --> 00:44:04,990 for any very small time steps. 897 00:44:04,990 --> 00:44:08,200 It makes more sense to choose a naturally implicit method 898 00:44:08,200 --> 00:44:11,830 for the differential equation, where I can take big time steps 899 00:44:11,830 --> 00:44:16,376 and still have reasonable stability of the integration. 900 00:44:16,376 --> 00:44:17,250 Does that make sense? 901 00:44:21,040 --> 00:44:25,180 Consider now the fully implicit DAE. 902 00:44:25,180 --> 00:44:28,480 So f of x, x dot, and t is equal to 0 903 00:44:28,480 --> 00:44:30,565 and have some set of initial conditions. x of 0 904 00:44:30,565 --> 00:44:32,590 is equal to x 0. 905 00:44:32,590 --> 00:44:34,960 Here there is no way in general to avoid solving 906 00:44:34,960 --> 00:44:37,330 systems of nonlinear equations. 907 00:44:37,330 --> 00:44:39,850 And so one substitutes often backwards 908 00:44:39,850 --> 00:44:42,670 difference approximations for x dot. 909 00:44:42,670 --> 00:44:45,355 And then you solve for the next point in time. 910 00:44:45,355 --> 00:44:47,230 So here's an example with the backward Euler, 911 00:44:47,230 --> 00:44:48,730 implicit Euler approximations. 912 00:44:48,730 --> 00:44:51,070 So I replaced the time derivative 913 00:44:51,070 --> 00:44:55,000 with the difference between the state at point tk 914 00:44:55,000 --> 00:44:57,160 and the state of point tk minus 1 915 00:44:57,160 --> 00:45:00,880 divided by the difference between tk and tk minus 1. 916 00:45:00,880 --> 00:45:03,670 The other state variables are evaluated at tk. 917 00:45:03,670 --> 00:45:05,130 So x at tk. 918 00:45:05,130 --> 00:45:06,430 And I evaluate time at tk. 919 00:45:06,430 --> 00:45:08,650 And I try to satisfy this equation 920 00:45:08,650 --> 00:45:10,460 and use it to determine x of tk. 921 00:45:10,460 --> 00:45:13,690 So I solve this system of nonlinear equations 922 00:45:13,690 --> 00:45:15,160 for x at tk. 923 00:45:15,160 --> 00:45:18,580 And then I go to the next point in time, tk plus 1. 924 00:45:18,580 --> 00:45:22,440 And I solve the same equation, but replacing tk with tk plus 1 925 00:45:22,440 --> 00:45:25,060 and tk minus 1 with tk. 926 00:45:25,060 --> 00:45:28,060 So at each iteration I solve a system of nonlinear equations. 927 00:45:28,060 --> 00:45:30,400 No more painful than doing the forward Euler 928 00:45:30,400 --> 00:45:32,860 approximation I showed you before. 929 00:45:32,860 --> 00:45:35,590 So that's good. 930 00:45:35,590 --> 00:45:38,865 But still a lot of work. 931 00:45:38,865 --> 00:45:40,240 We might wonder how do we control 932 00:45:40,240 --> 00:45:42,872 the error in these sorts of approximations. 933 00:45:42,872 --> 00:45:44,830 And that's what we're going to talk about next. 934 00:45:48,394 --> 00:45:50,310 So here's the equivalent time marching formula 935 00:45:50,310 --> 00:45:52,770 applied to a semi-explicit DAE. 936 00:45:55,910 --> 00:46:02,790 So I've got to satisfy the equation 937 00:46:02,790 --> 00:46:05,850 0 is equal to the derivative of x with respect to time, 938 00:46:05,850 --> 00:46:09,060 which I approximate with the backward Euler approximation, 939 00:46:09,060 --> 00:46:12,300 minus f evaluated at tk. 940 00:46:12,300 --> 00:46:19,410 And the equation g of x at tk, y at tk for this x of tk. 941 00:46:19,410 --> 00:46:20,670 And actually y of tk as well. 942 00:46:20,670 --> 00:46:22,800 I'm sorry, I left that off. y of tk as well. 943 00:46:22,800 --> 00:46:24,890 Got to solve for the differential 944 00:46:24,890 --> 00:46:28,340 and the algebraic states simultaneously. 945 00:46:28,340 --> 00:46:30,030 And I use a backwards difference formula 946 00:46:30,030 --> 00:46:31,877 to try to do this in a stable way. 947 00:46:36,850 --> 00:46:43,180 One way to think about in particular these sorts of DAEs 948 00:46:43,180 --> 00:46:46,990 is as exceptionally stiff ordinary differential 949 00:46:46,990 --> 00:46:47,490 equations. 950 00:46:47,490 --> 00:46:50,320 Did you guys discuss stiffness? 951 00:46:50,320 --> 00:46:55,750 So imagine if I replace the 0 with, say, dy dt. 952 00:46:55,750 --> 00:46:58,090 And I introduced a constant on the right-hand side, 953 00:46:58,090 --> 00:47:05,485 a multiplicative constant, say, 1 over lambda here. 954 00:47:05,485 --> 00:47:06,610 Let's say lambda, actually. 955 00:47:06,610 --> 00:47:08,610 Say I introduce a multiplicative constant lambda 956 00:47:08,610 --> 00:47:10,940 on the right-hand side here. 957 00:47:10,940 --> 00:47:15,010 So as lambda becomes very large, the system 958 00:47:15,010 --> 00:47:17,620 becomes increasingly stiff. 959 00:47:17,620 --> 00:47:19,600 The dynamics of this equation take 960 00:47:19,600 --> 00:47:22,790 place over very, very short time scales. 961 00:47:22,790 --> 00:47:24,820 So short in fact, that eventually, 962 00:47:24,820 --> 00:47:26,500 if I let lambda diverge to infinity, 963 00:47:26,500 --> 00:47:29,500 I've just got to satisfy this algebraic equation 964 00:47:29,500 --> 00:47:31,400 at each point in time. 965 00:47:31,400 --> 00:47:34,240 So these are exceptionally stiff equations. 966 00:47:34,240 --> 00:47:36,270 And the methods for solving these equations 967 00:47:36,270 --> 00:47:39,880 share a lot in common with stiff ODE-IVP solvers. 968 00:47:45,069 --> 00:47:46,610 So I've got a couple of minutes left. 969 00:47:46,610 --> 00:47:48,880 I'm going to work through some examples. 970 00:47:48,880 --> 00:47:51,540 So consider the stirred tank example one. 971 00:47:51,540 --> 00:47:53,250 dc 2 dt. 972 00:47:53,250 --> 00:47:57,210 It's related to q over v, flow rate over volume. 973 00:47:57,210 --> 00:47:58,770 The difference in the concentrations 974 00:47:58,770 --> 00:48:00,570 coming in and out of the tank. 975 00:48:00,570 --> 00:48:03,060 And c1 is equal to this input signal gamma. 976 00:48:03,060 --> 00:48:06,030 Maybe that's a measurement and I'm trying to predict c2. 977 00:48:06,030 --> 00:48:08,850 So can you apply the backward Euler method to these equations 978 00:48:08,850 --> 00:48:12,810 to determine c1 and c2? 979 00:48:12,810 --> 00:48:15,670 What would one step of the backward Euler method look 980 00:48:15,670 --> 00:48:24,500 like, say, going from time 0 to time delta t, or time t to t 981 00:48:24,500 --> 00:48:27,260 plus delta t, or tk to tk plus 1, however you 982 00:48:27,260 --> 00:48:28,748 want to write it? 983 00:49:22,314 --> 00:49:23,730 So do you get something like this? 984 00:49:26,480 --> 00:49:28,530 So this equation is easy. 985 00:49:28,530 --> 00:49:32,360 c1 to tk is gamma of tk. 986 00:49:32,360 --> 00:49:34,390 I got a substitute up here. 987 00:49:34,390 --> 00:49:38,570 c2 of tk minus c2 of tk minus 1 over delta t. 988 00:49:38,570 --> 00:49:41,210 T Just tk minus tk minus 1. 989 00:49:41,210 --> 00:49:42,800 And then solve for c2 of dk. 990 00:49:42,800 --> 00:49:44,390 And you get a formula like this. 991 00:49:49,010 --> 00:49:52,010 This formula was based on an approximation 992 00:49:52,010 --> 00:49:56,510 for the derivative, which had a leading order 993 00:49:56,510 --> 00:50:03,140 error that was proportional to tk minus tk minus 1. 994 00:50:03,140 --> 00:50:07,070 When I go through and I solve for c2 of tk, the leading order 995 00:50:07,070 --> 00:50:10,600 error in c2, the local truncation error 996 00:50:10,600 --> 00:50:16,250 is going to be order tk minus tk minus 1 squared. 997 00:50:16,250 --> 00:50:20,840 So as I shrink the spacing between each of my two time 998 00:50:20,840 --> 00:50:24,050 points, the local accuracy, the error 999 00:50:24,050 --> 00:50:27,830 induced in one advance of this time stepping algorithm 1000 00:50:27,830 --> 00:50:29,790 is going to scale like the step size squared. 1001 00:50:29,790 --> 00:50:31,550 So I have the error-- 1002 00:50:31,550 --> 00:50:32,369 I have the spacing. 1003 00:50:32,369 --> 00:50:34,160 The error will go down by a factor of four. 1004 00:50:34,160 --> 00:50:37,667 I reduce the spacing by an order of magnitude. 1005 00:50:37,667 --> 00:50:39,750 The error will go down by two orders of magnitude. 1006 00:50:39,750 --> 00:50:43,605 This is a really desirable sort of error control and a method. 1007 00:50:43,605 --> 00:50:44,480 And we can do better. 1008 00:50:44,480 --> 00:50:48,460 You've seen us do better with other methods. 1009 00:50:48,460 --> 00:50:51,180 Let's look at this equation now. 1010 00:50:51,180 --> 00:50:52,590 This is stirred tank example two. 1011 00:50:52,590 --> 00:50:53,797 I just replace c1 with c2. 1012 00:50:53,797 --> 00:50:55,380 You don't have to write this one down. 1013 00:50:55,380 --> 00:50:56,484 I'll just show it to you. 1014 00:50:56,484 --> 00:50:58,650 The notes are online so you're not missing anything. 1015 00:50:58,650 --> 00:51:00,900 All these things will pop up in the notes online. 1016 00:51:00,900 --> 00:51:05,220 So c2 at tk is equal to gamma of tk. 1017 00:51:05,220 --> 00:51:06,930 I got to come up to my first equation now 1018 00:51:06,930 --> 00:51:08,700 and I have to substitute an approximation 1019 00:51:08,700 --> 00:51:09,810 for the derivative. 1020 00:51:09,810 --> 00:51:13,390 It's this, the backwards difference formula. 1021 00:51:13,390 --> 00:51:17,610 So if find c1 of tk is c2 of tk plus v over q. 1022 00:51:17,610 --> 00:51:21,660 c2 of tk minus c2 of tk minus 1 divided by the time step. 1023 00:51:21,660 --> 00:51:23,640 Remember, this approximation for the derivative 1024 00:51:23,640 --> 00:51:27,292 has a leading order error that goes like tk minus tk minus 1. 1025 00:51:27,292 --> 00:51:32,235 So my approximation for c1, the local approximation for c1, 1026 00:51:32,235 --> 00:51:34,040 is not second order accurate. 1027 00:51:34,040 --> 00:51:38,010 It's only first order accurate. 1028 00:51:38,010 --> 00:51:42,450 I applied the same approximation method to both equations. 1029 00:51:42,450 --> 00:51:45,050 One equation I got something that was second order accurate. 1030 00:51:45,050 --> 00:51:46,425 The next equation I got something 1031 00:51:46,425 --> 00:51:47,840 that was first order accurate. 1032 00:51:47,840 --> 00:51:49,520 They look almost identical. 1033 00:51:49,520 --> 00:51:51,530 It might be hard to see from the start 1034 00:51:51,530 --> 00:51:53,420 why it should work out that way. 1035 00:51:53,420 --> 00:51:56,190 That's how it works out. 1036 00:51:56,190 --> 00:51:57,410 Let's do one more example. 1037 00:52:00,640 --> 00:52:03,680 Here's a system of DAEs. 1038 00:52:03,680 --> 00:52:05,760 c2 dot is equal to c1. 1039 00:52:05,760 --> 00:52:07,900 c3 dot is equal to c2. 1040 00:52:07,900 --> 00:52:10,240 And 0 is equal to c3 minus gamma. 1041 00:52:10,240 --> 00:52:12,250 So gamma, again, is some measurement. 1042 00:52:12,250 --> 00:52:14,660 We solve this equation to determine c3. 1043 00:52:14,660 --> 00:52:17,640 We substitute c3 up here to determine c2. 1044 00:52:17,640 --> 00:52:20,600 We substitute c2 up here to determine c1. 1045 00:52:20,600 --> 00:52:23,965 But we want to find the solution with the backward Euler method. 1046 00:52:27,100 --> 00:52:32,680 And the solution there is going to be c3 of tk is gamma of tk. 1047 00:52:32,680 --> 00:52:35,410 c2 is related to the derivative of c3. 1048 00:52:35,410 --> 00:52:37,180 So I use my backwards difference formula 1049 00:52:37,180 --> 00:52:39,060 for the derivative of c3. 1050 00:52:39,060 --> 00:52:41,620 c1 is related to the derivative of c2. 1051 00:52:41,620 --> 00:52:44,710 So I use my backward difference formula for c2. 1052 00:52:44,710 --> 00:52:46,930 But how accurate are each of these approximations? 1053 00:52:49,840 --> 00:52:52,780 So here I have to determine c2 from an approximation 1054 00:52:52,780 --> 00:52:54,880 for the derivative of c3. 1055 00:52:54,880 --> 00:52:56,920 We know this backwards difference approximation 1056 00:52:56,920 --> 00:52:58,510 has a leading order error proportional 1057 00:52:58,510 --> 00:53:02,350 to delta t, the time step. 1058 00:53:02,350 --> 00:53:07,750 Now I got to approximate c1 with the finite difference backwards 1059 00:53:07,750 --> 00:53:09,922 difference approximation in c2. 1060 00:53:09,922 --> 00:53:12,380 We know this should incur an error proportional to the time 1061 00:53:12,380 --> 00:53:17,350 step delta t, except I don't know c2 at tk exactly. 1062 00:53:17,350 --> 00:53:23,530 Actually, c2 tk, the exact value of c2 at tk 1063 00:53:23,530 --> 00:53:27,830 is this plus something proportional to the time step. 1064 00:53:27,830 --> 00:53:30,560 Which means the exact value of c1 1065 00:53:30,560 --> 00:53:34,430 is this plus something that's order one because I 1066 00:53:34,430 --> 00:53:36,299 carry over the error from my solution 1067 00:53:36,299 --> 00:53:37,340 to the previous equation. 1068 00:53:37,340 --> 00:53:40,477 And I divide it by the time step. 1069 00:53:40,477 --> 00:53:42,560 So I went from having something-- well, let's see. 1070 00:53:42,560 --> 00:53:44,340 The first equation had order delta t 1071 00:53:44,340 --> 00:53:47,040 squared local truncation error. 1072 00:53:47,040 --> 00:53:50,010 As I shrink the spacing, my solution 1073 00:53:50,010 --> 00:53:52,440 gets more and more locally accurate. 1074 00:53:52,440 --> 00:53:57,300 And it seems to do it in a fairly robust fashion. 1075 00:53:57,300 --> 00:53:59,910 The next system of equations we solved 1076 00:53:59,910 --> 00:54:02,170 had a local truncation error proportional to delta t. 1077 00:54:02,170 --> 00:54:03,545 Again, I can shrink delta t and I 1078 00:54:03,545 --> 00:54:05,800 can make my solution more and more accurate. 1079 00:54:05,800 --> 00:54:08,220 That seems like what you want in a numerical algorithm. 1080 00:54:08,220 --> 00:54:11,580 This one, there's no hope. 1081 00:54:11,580 --> 00:54:15,980 Change delta t, who knows what you're going to get? 1082 00:54:15,980 --> 00:54:18,855 And this looks fairly innocuous. 1083 00:54:18,855 --> 00:54:20,990 You just look at it and you say, well, OK. 1084 00:54:20,990 --> 00:54:22,980 Let's solve it and see what happens. 1085 00:54:22,980 --> 00:54:26,040 Well actually, you can't solve this problem accurately 1086 00:54:26,040 --> 00:54:28,322 with the backward Euler method. 1087 00:54:28,322 --> 00:54:30,030 There's something fundamentally different 1088 00:54:30,030 --> 00:54:32,359 about this problem from the previous one. 1089 00:54:32,359 --> 00:54:34,650 There's something fundamentally different about stirred 1090 00:54:34,650 --> 00:54:37,500 tank example two from stirred tank example one. 1091 00:54:37,500 --> 00:54:40,170 Right you might have already perceived that. 1092 00:54:40,170 --> 00:54:42,160 You can see what's going on here. 1093 00:54:42,160 --> 00:54:44,520 c3 is equal to gamma. 1094 00:54:44,520 --> 00:54:46,025 That's the exact solution. 1095 00:54:46,025 --> 00:54:48,390 c2 is equal to gamma dot. 1096 00:54:48,390 --> 00:54:50,700 And c1 is equal to two derivatives 1097 00:54:50,700 --> 00:54:53,860 of gamma, the exact solution to this equation. 1098 00:54:53,860 --> 00:54:59,310 So c1 is incredibly sensitive to changes in gamma. 1099 00:54:59,310 --> 00:55:00,620 And that's peculiar. 1100 00:55:00,620 --> 00:55:03,180 And it has important numerical consequences. 1101 00:55:03,180 --> 00:55:06,380 So I'm going to show you, our next lecture on Wednesday, 1102 00:55:06,380 --> 00:55:09,620 how to formalize an understanding 1103 00:55:09,620 --> 00:55:11,450 of the differences between these problems. 1104 00:55:11,450 --> 00:55:14,990 And we'll talk about some more subtleties associated 1105 00:55:14,990 --> 00:55:17,150 with differential algebraic equations. 1106 00:55:17,150 --> 00:55:18,650 You guys can pick up your quizzes. 1107 00:55:18,650 --> 00:55:23,000 We should do it outside rather than 1108 00:55:23,000 --> 00:55:25,596 inside because the next class has to start.