1 00:00:00,030 --> 00:00:02,400 The following content is provided under a Creative 2 00:00:02,400 --> 00:00:03,770 Commons license. 3 00:00:03,770 --> 00:00:06,020 Your support will help MIT OpenCourseWare 4 00:00:06,020 --> 00:00:10,100 continue to offer high quality educational resources for free. 5 00:00:10,100 --> 00:00:12,670 To make a donation or to view additional materials 6 00:00:12,670 --> 00:00:16,580 from hundreds of MIT courses, visit MIT OpenCourseWare 7 00:00:16,580 --> 00:00:17,235 at ocw.mit.edu. 8 00:00:26,194 --> 00:00:27,360 QIQI WANG: Hello, everybody. 9 00:00:29,940 --> 00:00:37,000 So I guess yesterday was the a really serious reading due. 10 00:00:37,000 --> 00:00:38,294 How is that going? 11 00:00:38,294 --> 00:00:39,460 How do you like the reading? 12 00:00:39,460 --> 00:00:40,668 AUDIENCE: What's consistency? 13 00:00:40,668 --> 00:00:41,450 QIQI WANG: Hm? 14 00:00:41,450 --> 00:00:43,130 AUDIENCE: What's consistency? 15 00:00:43,130 --> 00:00:45,590 QIQI WANG: What's consistency? 16 00:00:45,590 --> 00:00:52,100 Consistency is the scheme being at least first order accurate. 17 00:00:52,100 --> 00:00:52,600 OK? 18 00:00:52,600 --> 00:00:56,350 So the it's local-- consistently is related 19 00:00:56,350 --> 00:00:58,390 to local order of accuracy. 20 00:00:58,390 --> 00:01:03,600 So if the scheme is at least first order local accuracy, 21 00:01:03,600 --> 00:01:04,806 then it is consistent. 22 00:01:04,806 --> 00:01:05,722 AUDIENCE: [INAUDIBLE]. 23 00:01:08,352 --> 00:01:10,060 QIQI WANG: It means the truncating error. 24 00:01:10,060 --> 00:01:12,345 If you do the truncating error analysis, 25 00:01:12,345 --> 00:01:14,700 the data series gives you a truncating error 26 00:01:14,700 --> 00:01:20,410 of at least delta t squared if you don't divide by delta t. 27 00:01:20,410 --> 00:01:23,530 At least delta t if you divide by delta t. 28 00:01:23,530 --> 00:01:26,680 So that's what it means by being consistent. 29 00:01:26,680 --> 00:01:28,730 I mean the word consistent basically 30 00:01:28,730 --> 00:01:32,935 means if you look at the scheme, it 31 00:01:32,935 --> 00:01:35,810 is somewhat consistent to the differential equation. 32 00:01:35,810 --> 00:01:40,290 So the scheme on one hand, the differential equation 33 00:01:40,290 --> 00:01:41,150 on the other hand. 34 00:01:41,150 --> 00:01:43,724 They are consistent, right? 35 00:01:43,724 --> 00:01:46,172 AUDIENCE: So that would mean that the scheme is 36 00:01:46,172 --> 00:01:49,417 self-satisfied with the [INAUDIBLE]? 37 00:01:49,417 --> 00:01:51,250 QIQI WANG: The scheme here-- so the question 38 00:01:51,250 --> 00:01:54,170 is does that mean the scheme satisfies the differential 39 00:01:54,170 --> 00:01:54,670 equation? 40 00:01:54,670 --> 00:01:56,500 The answer is the scheme is never 41 00:01:56,500 --> 00:01:59,090 going to satisfy the differential equation exactly. 42 00:01:59,090 --> 00:02:02,580 But as you make delta t smaller and smaller, 43 00:02:02,580 --> 00:02:05,060 it is going to satisfy the differential equation better 44 00:02:05,060 --> 00:02:06,640 and better. 45 00:02:06,640 --> 00:02:07,140 OK. 46 00:02:07,140 --> 00:02:10,644 That's what it means by being consistent. 47 00:02:10,644 --> 00:02:12,185 AUDIENCE: So if it wasn't consistent, 48 00:02:12,185 --> 00:02:14,085 it would make the [INAUDIBLE]. 49 00:02:17,614 --> 00:02:18,239 QIQI WANG: Yes. 50 00:02:18,239 --> 00:02:18,739 So, right. 51 00:02:18,739 --> 00:02:23,090 If a scheme is not consistent, what does that mean? 52 00:02:23,090 --> 00:02:26,020 That means if you make the delta t smaller and smaller, 53 00:02:26,020 --> 00:02:28,956 the difference between the discrete scheme 54 00:02:28,956 --> 00:02:30,580 and the continual differential equation 55 00:02:30,580 --> 00:02:32,510 is not going to become smaller. 56 00:02:32,510 --> 00:02:33,620 It may become larger. 57 00:02:33,620 --> 00:02:34,725 It may remain constant. 58 00:02:34,725 --> 00:02:37,053 It may go around but it doesn't go to zero. 59 00:02:37,053 --> 00:02:39,136 AUDIENCE: How does that differ from the definition 60 00:02:39,136 --> 00:02:40,266 of stability? 61 00:02:40,266 --> 00:02:42,570 QIQI WANG: How does that differ from the definition 62 00:02:42,570 --> 00:02:44,060 of stability? 63 00:02:44,060 --> 00:02:46,840 We're going to see there are-- when you talk about stability, 64 00:02:46,840 --> 00:02:49,120 there are two types of stability. 65 00:02:49,120 --> 00:02:55,720 They are different in the sense that a consistency concerns 66 00:02:55,720 --> 00:02:57,964 the Taylor series analysis. 67 00:02:57,964 --> 00:03:02,620 You can figure out consistency by doing data series analysis. 68 00:03:02,620 --> 00:03:07,250 Stability has nothing to do with data series. 69 00:03:07,250 --> 00:03:07,750 Right? 70 00:03:07,750 --> 00:03:10,435 We are going to see when you check stability, 71 00:03:10,435 --> 00:03:13,860 you're plugging in a half solution that 72 00:03:13,860 --> 00:03:20,310 grows exponentially and seeing can the exponential 73 00:03:20,310 --> 00:03:22,710 be a growing exponential? 74 00:03:22,710 --> 00:03:27,150 Or the exponential has to be a decaying exponential? 75 00:03:27,150 --> 00:03:29,290 If all the exponential solutions have 76 00:03:29,290 --> 00:03:38,320 to be decaying exponentials, then the scheme is stable. 77 00:03:38,320 --> 00:03:44,910 If there is one mode that can be a growing exponential, then 78 00:03:44,910 --> 00:03:45,948 the scheme is unstable. 79 00:03:49,010 --> 00:03:49,510 OK. 80 00:03:49,510 --> 00:03:53,380 I clarified a question about consistency. 81 00:03:53,380 --> 00:03:56,894 That is basically the scheme is at least first order accurate. 82 00:03:56,894 --> 00:04:00,490 What is the difference between consistency and stability? 83 00:04:00,490 --> 00:04:03,390 Consistency is by data series analysis. 84 00:04:03,390 --> 00:04:07,390 Stability is by plugging exponentially growing 85 00:04:07,390 --> 00:04:08,826 solutions. 86 00:04:08,826 --> 00:04:10,780 So basically asking does the scheme 87 00:04:10,780 --> 00:04:13,400 allow exponentially growing solutions? 88 00:04:18,970 --> 00:04:21,290 Any other questions that arise in the reading? 89 00:04:21,290 --> 00:04:25,200 If one person has a concern, probably your classmate 90 00:04:25,200 --> 00:04:26,080 also has it. 91 00:04:26,080 --> 00:04:28,366 So please raise it. 92 00:04:28,366 --> 00:04:30,824 AUDIENCE: Maybe one other way to think about the definition 93 00:04:30,824 --> 00:04:41,800 of stability is analysis. 94 00:04:41,800 --> 00:04:46,450 [INAUDIBLE] the solution behaves mostly [INAUDIBLE]. 95 00:04:46,450 --> 00:04:48,912 When we think about stability, meaning that we're 96 00:04:48,912 --> 00:04:49,912 taking the relationship. 97 00:04:49,912 --> 00:04:59,872 We're looking at how [INAUDIBLE] analysis. 98 00:04:59,872 --> 00:05:01,366 Other things [INAUDIBLE]. 99 00:05:06,844 --> 00:05:09,334 The other is [INAUDIBLE]. 100 00:05:09,334 --> 00:05:11,824 Because as w goes to 0, [INAUDIBLE]. 101 00:05:16,804 --> 00:05:22,771 Because I think more and more [INAUDIBLE] really [INAUDIBLE]. 102 00:05:22,771 --> 00:05:25,312 You could probably think there's a few different [INAUDIBLE]. 103 00:05:25,312 --> 00:05:27,256 You're going to need both of those in order 104 00:05:27,256 --> 00:05:28,714 to analyze [INAUDIBLE]. 105 00:05:34,905 --> 00:05:35,530 QIQI WANG: Yes. 106 00:05:35,530 --> 00:05:36,213 Thank you. 107 00:05:36,213 --> 00:05:39,640 Thank you for the [INAUDIBLE]. 108 00:05:39,640 --> 00:05:42,290 Any other questions? 109 00:05:42,290 --> 00:05:43,160 No? 110 00:05:43,160 --> 00:05:46,940 Then I basically want to very quickly review and recap 111 00:05:46,940 --> 00:05:48,960 what we did in the previous lectures 112 00:05:48,960 --> 00:05:52,580 and also what we read before today. 113 00:05:52,580 --> 00:05:54,850 So we did local order of accuracy. 114 00:05:54,850 --> 00:05:59,920 That is basically by plugging Taylor series 115 00:05:59,920 --> 00:06:04,600 into the numerical scheme and figuring out how, 116 00:06:04,600 --> 00:06:05,850 which is the truncating error. 117 00:06:05,850 --> 00:06:07,349 Is it the left hand side [INAUDIBLE] 118 00:06:07,349 --> 00:06:11,520 right hand side of the numerical scheme? 119 00:06:11,520 --> 00:06:14,190 When you plug in both the differential equation 120 00:06:14,190 --> 00:06:19,720 and the data series, does it go-- what order of accuracy 121 00:06:19,720 --> 00:06:21,210 does it go? 122 00:06:21,210 --> 00:06:24,844 So as delta t goes to 0, how fast 123 00:06:24,844 --> 00:06:29,600 does my inconsistency between the differential equation 124 00:06:29,600 --> 00:06:33,170 and my numerical scheme go to 0? 125 00:06:33,170 --> 00:06:34,280 OK. 126 00:06:34,280 --> 00:06:37,860 And the local order of accuracy is closely related 127 00:06:37,860 --> 00:06:39,380 to consistency. 128 00:06:39,380 --> 00:06:45,695 Actually the definition of consistency is by looking at 129 00:06:45,695 --> 00:06:49,660 is the scheme at least a first order accurate or not. 130 00:06:49,660 --> 00:06:50,900 OK. 131 00:06:50,900 --> 00:06:52,130 Global accuracy. 132 00:06:52,130 --> 00:06:55,230 You should have read about global accuracy. 133 00:06:55,230 --> 00:06:57,390 Could somebody answer me the question 134 00:06:57,390 --> 00:07:00,020 of what is the difference between global accuracy 135 00:07:00,020 --> 00:07:01,795 and local accuracy? 136 00:07:01,795 --> 00:07:02,670 Or are they the same? 137 00:07:02,670 --> 00:07:03,292 Yes? 138 00:07:03,292 --> 00:07:04,167 AUDIENCE: [INAUDIBLE] 139 00:07:18,660 --> 00:07:20,150 QIQI WANG: So the global accuracy 140 00:07:20,150 --> 00:07:26,243 is the long term accuracy while the local accuracy is not 141 00:07:26,243 --> 00:07:26,826 the long term? 142 00:07:30,714 --> 00:07:32,307 AUDIENCE: Local accuracy. 143 00:07:32,307 --> 00:07:33,890 QIQI WANG: The local order of accuracy 144 00:07:33,890 --> 00:07:38,750 you can figure out by looking at only one time step, right? 145 00:07:38,750 --> 00:07:41,350 The global order of accuracy, you 146 00:07:41,350 --> 00:07:44,660 have to look at really by taking the number of time steps 147 00:07:44,660 --> 00:07:47,110 to infinity. 148 00:07:47,110 --> 00:07:52,310 And as Professor Willcox said, if you look at the solution U 149 00:07:52,310 --> 00:07:56,620 at a particular time-- let's say U at t equal to 1-- 150 00:07:56,620 --> 00:08:00,470 and compare the numerical solution-- let me say V at t 151 00:08:00,470 --> 00:08:02,260 equal to 1-- right? 152 00:08:02,260 --> 00:08:04,670 If you look at the difference between them, 153 00:08:04,670 --> 00:08:09,330 let's just take the absolute value of the difference. 154 00:08:09,330 --> 00:08:13,360 That difference is really looking 155 00:08:13,360 --> 00:08:18,970 at how much error has the solution accumulated 156 00:08:18,970 --> 00:08:24,045 over one time step all the way from t equals 0 to t equals 1? 157 00:08:24,045 --> 00:08:26,780 Right? 158 00:08:26,780 --> 00:08:33,460 And if I look at that area as I take delta t to infinity, 159 00:08:33,460 --> 00:08:37,470 I'm actually accumulating more time steps, right? 160 00:08:37,470 --> 00:08:40,850 So if delta t is to 1 I'm accumulating 10 time steps. 161 00:08:40,850 --> 00:08:42,659 If delta t is [INAUDIBLE], then I'm 162 00:08:42,659 --> 00:08:45,285 looking at the total area that has been accumulated 163 00:08:45,285 --> 00:08:47,430 over 100 time steps. 164 00:08:47,430 --> 00:08:51,604 Now we can for sure see that this 165 00:08:51,604 --> 00:08:55,930 puts a more stringent requirement on how 166 00:08:55,930 --> 00:08:57,380 my scheme behaves. 167 00:08:57,380 --> 00:09:00,400 Because when I only look at local order of accuracy 168 00:09:00,400 --> 00:09:04,070 and just look at one time step, now I'm looking at not only 169 00:09:04,070 --> 00:09:08,190 how much area have I produced in one time step 170 00:09:08,190 --> 00:09:12,810 but also how much has the area grown? 171 00:09:12,810 --> 00:09:14,722 I mean, does the area made in one time 172 00:09:14,722 --> 00:09:17,000 step actually only matter locally? 173 00:09:17,000 --> 00:09:19,930 Or is this area actually amplified 174 00:09:19,930 --> 00:09:22,430 by the later time steps? 175 00:09:22,430 --> 00:09:26,830 So it puts additional requirements on a scheme. 176 00:09:26,830 --> 00:09:33,430 A scheme is globally accurate only if it is 177 00:09:33,430 --> 00:09:36,980 both locally accurate and what? 178 00:09:43,620 --> 00:09:44,540 I hear two answers. 179 00:09:44,540 --> 00:09:46,000 One says stable. 180 00:09:46,000 --> 00:09:50,020 Another is consistent. 181 00:09:50,020 --> 00:09:51,106 Who votes for stable? 182 00:09:51,106 --> 00:09:52,970 And who votes for consistent? 183 00:09:52,970 --> 00:09:55,620 Stable? 184 00:09:55,620 --> 00:09:57,740 One, two, three, four. 185 00:09:57,740 --> 00:10:00,030 Consistent? 186 00:10:00,030 --> 00:10:02,640 One, two, three, four, five. 187 00:10:02,640 --> 00:10:04,280 OK. 188 00:10:04,280 --> 00:10:06,830 I said a scheme is globally accurate 189 00:10:06,830 --> 00:10:09,600 only if it's locally accurate and what? 190 00:10:09,600 --> 00:10:12,200 So if we the answer is consistent, 191 00:10:12,200 --> 00:10:14,000 then it doesn't add anything, right? 192 00:10:14,000 --> 00:10:18,430 Being consistent means being locally accurate. 193 00:10:18,430 --> 00:10:21,210 OK? 194 00:10:21,210 --> 00:10:23,760 A scheme is locally accurate and consistent means 195 00:10:23,760 --> 00:10:24,885 it's just locally accurate. 196 00:10:24,885 --> 00:10:26,470 It doesn't say anything additional 197 00:10:26,470 --> 00:10:28,480 to being locally accurate. 198 00:10:28,480 --> 00:10:30,370 That's the definition of being consistent. 199 00:10:30,370 --> 00:10:33,935 It is just to be locally accurate. 200 00:10:33,935 --> 00:10:37,930 A scheme is globally accurate only if it's locally accurate 201 00:10:37,930 --> 00:10:40,980 and is stable. 202 00:10:40,980 --> 00:10:42,640 Actually zero stable. 203 00:10:42,640 --> 00:10:44,330 And what does zero stable mean? 204 00:10:44,330 --> 00:10:46,241 You should have just read about it. 205 00:10:49,910 --> 00:10:52,520 I'm just going to say I mean you have read about it. 206 00:10:52,520 --> 00:10:55,320 I'm just going to say another way of understanding 207 00:10:55,320 --> 00:10:59,890 zero accuracy is that when you solve the differential equation 208 00:10:59,890 --> 00:11:05,660 du/dt equal to zero, the scheme is going to be stable. 209 00:11:05,660 --> 00:11:06,360 OK? 210 00:11:06,360 --> 00:11:10,016 That's what it means by a scheme being zero stable. 211 00:11:13,860 --> 00:11:17,020 In this sense, it's a very weak requirement. 212 00:11:17,020 --> 00:11:18,710 Right? 213 00:11:18,710 --> 00:11:20,840 You only need to be stable when solving 214 00:11:20,840 --> 00:11:24,650 this trivial differential equation. 215 00:11:24,650 --> 00:11:25,150 Right? 216 00:11:25,150 --> 00:11:29,060 What is the solution for a differential equation? 217 00:11:29,060 --> 00:11:30,060 AUDIENCE: Confidence. 218 00:11:30,060 --> 00:11:31,910 QIQI WANG: Confidence, right. 219 00:11:31,910 --> 00:11:34,240 I mean, if you're so very [INAUDIBLE], 220 00:11:34,240 --> 00:11:36,880 even your scheme can't even solve that differential 221 00:11:36,880 --> 00:11:40,880 equation, what good is the scheme? 222 00:11:40,880 --> 00:11:41,445 Right? 223 00:11:41,445 --> 00:11:41,945 OK. 224 00:11:45,440 --> 00:11:46,966 So this is, yes. 225 00:11:46,966 --> 00:11:53,080 It is an additional requirement on top of local accuracy. 226 00:11:53,080 --> 00:11:55,090 But it is a pretty weak requirement. 227 00:11:55,090 --> 00:11:56,370 All right? 228 00:11:56,370 --> 00:12:00,450 That is why I think it is such a great result that we 229 00:12:00,450 --> 00:12:03,670 can say being locally accurate [INAUDIBLE]. 230 00:12:06,444 --> 00:12:09,470 And being zero stable, which is such a weak requirement, 231 00:12:09,470 --> 00:12:14,930 you can immediately have the conclusion of global accuracy. 232 00:12:14,930 --> 00:12:16,160 Right? 233 00:12:16,160 --> 00:12:22,324 So what theorem or result gives you that conclusion? 234 00:12:22,324 --> 00:12:23,240 AUDIENCE: [INAUDIBLE]. 235 00:12:32,060 --> 00:12:33,530 QIQI WANG: What is it called? 236 00:12:33,530 --> 00:12:35,500 Dahlquist Equivalence Theorem. 237 00:12:35,500 --> 00:12:36,620 Yes. 238 00:12:36,620 --> 00:12:37,370 Right. 239 00:12:37,370 --> 00:12:39,420 So that theorem basically says that these 240 00:12:39,420 --> 00:12:45,140 are locally accurate, which only looks at one time step. 241 00:12:45,140 --> 00:12:48,022 And if you have are zero stable-- which 242 00:12:48,022 --> 00:12:53,020 means you can solve this trivial differential equation-- then 243 00:12:53,020 --> 00:12:55,040 your scheme is globally accurate which 244 00:12:55,040 --> 00:12:57,420 means you can solve the differential 245 00:12:57,420 --> 00:13:01,246 equation for a fixed amount of time 246 00:13:01,246 --> 00:13:03,490 as it takes delta t to go to zero, 247 00:13:03,490 --> 00:13:08,250 which also means the amount of time steps goes to infinity. 248 00:13:08,250 --> 00:13:13,180 Your global solution becomes more and more accurate. 249 00:13:13,180 --> 00:13:14,610 OK. 250 00:13:14,610 --> 00:13:17,150 And in addition to this, zero stability 251 00:13:17,150 --> 00:13:20,440 gives you the result that the local order of accuracy 252 00:13:20,440 --> 00:13:25,760 is the same as the global order of accuracy, right? 253 00:13:25,760 --> 00:13:27,800 So if you can do Taylor series analysis, 254 00:13:27,800 --> 00:13:29,960 you can not only know the local order of accuracy. 255 00:13:29,960 --> 00:13:32,395 You also know the global order of accuracy. 256 00:13:35,450 --> 00:13:36,120 OK. 257 00:13:36,120 --> 00:13:42,930 Eigenvalue stability is a step beyond zero stability. 258 00:13:42,930 --> 00:13:46,760 This tells you that your scheme has 259 00:13:46,760 --> 00:13:48,970 to be not only solving this equation 260 00:13:48,970 --> 00:13:54,760 but also du/dt equal to lambda U. Right? 261 00:13:54,760 --> 00:13:57,340 And then this is also a pretty easy differential equation. 262 00:13:57,340 --> 00:14:01,890 You have analytical solutions to that. 263 00:14:01,890 --> 00:14:06,598 The zero stability, you can think of as a special case. 264 00:14:06,598 --> 00:14:10,510 It's a much weaker requirement than Eigenvalue stability. 265 00:14:10,510 --> 00:14:13,113 Eigenvalue stability means you need 266 00:14:13,113 --> 00:14:16,490 to solve now nontrivial differential equations. 267 00:14:16,490 --> 00:14:16,990 OK. 268 00:14:16,990 --> 00:14:18,970 So today what we're going to do is 269 00:14:18,970 --> 00:14:22,640 we are going to exercise this local order of accuracy, 270 00:14:22,640 --> 00:14:24,600 global order of accuracy, and the zero 271 00:14:24,600 --> 00:14:26,860 stability, and Eigenvalue stability 272 00:14:26,860 --> 00:14:31,890 by looking at a pretty simple problem. 273 00:14:31,890 --> 00:14:34,450 It's a ballistic trajectory prediction problem. 274 00:14:34,450 --> 00:14:39,225 I hope most of you brought your computer. 275 00:14:39,225 --> 00:14:41,680 Does everybody have your computer with you? 276 00:14:44,340 --> 00:14:48,830 And do you have MATLAB installed? 277 00:14:48,830 --> 00:14:50,040 Yes? 278 00:14:50,040 --> 00:14:52,430 OK. 279 00:14:52,430 --> 00:14:54,260 So what I want you to do is, I want 280 00:14:54,260 --> 00:14:59,810 you to do mathematical modeling which is deriving some ODEs. 281 00:14:59,810 --> 00:15:01,390 I'm not going to give you any ODEs. 282 00:15:01,390 --> 00:15:02,950 You need to derive it. 283 00:15:02,950 --> 00:15:07,180 And you're going to solve the ODE using forward Euler. 284 00:15:07,180 --> 00:15:09,720 Solve the ODE using midpoint. 285 00:15:09,720 --> 00:15:13,940 And solve the ODE with one of the homework, 286 00:15:13,940 --> 00:15:17,890 the first homework question you did, 287 00:15:17,890 --> 00:15:23,530 which is the most accurate two-step explicit scheme. 288 00:15:23,530 --> 00:15:25,740 Anybody who can tell me what that scheme is? 289 00:15:29,080 --> 00:15:39,509 V of m plus 1 equal to minus 4 V m plus 5 V m minus 1. 290 00:15:39,509 --> 00:15:40,425 AUDIENCE: [INAUDIBLE]. 291 00:15:44,150 --> 00:15:48,120 QIQI WANG: Plus 4 delta t F of V m 292 00:15:48,120 --> 00:15:52,310 plus 2 delta t F of Vm minus 1. 293 00:15:52,310 --> 00:15:53,690 So I mean, yeah. 294 00:15:53,690 --> 00:15:57,270 It is the most accurate two-step scheme. 295 00:15:57,270 --> 00:15:59,620 But why doesn't it have a name? 296 00:15:59,620 --> 00:16:03,150 Something like this should have a name, right? 297 00:16:03,150 --> 00:16:07,170 Let's try to see why it doesn't have a name by applying this 298 00:16:07,170 --> 00:16:12,020 onto the ballistic trajectory prediction problem. 299 00:16:12,020 --> 00:16:12,950 All right. 300 00:16:12,950 --> 00:16:16,850 And let's look at accuracy, stability, convergence 301 00:16:16,850 --> 00:16:18,885 and how they relate to each other 302 00:16:18,885 --> 00:16:20,260 through the Dahlquist Equivalence 303 00:16:20,260 --> 00:16:24,881 Theorem using this example. 304 00:16:24,881 --> 00:16:25,380 OK. 305 00:16:25,380 --> 00:16:27,740 Ballistic trajectory predictions. 306 00:16:27,740 --> 00:16:28,910 Here is some history. 307 00:16:28,910 --> 00:16:32,140 It is really the motivating factor 308 00:16:32,140 --> 00:16:35,210 of developing the first real computer, 309 00:16:35,210 --> 00:16:37,830 the first real electrical computer is what they call it. 310 00:16:37,830 --> 00:16:41,290 They are really designed to help engineers 311 00:16:41,290 --> 00:16:43,740 solve this kind of problem. 312 00:16:43,740 --> 00:16:48,320 And here we are going to be programming a MATLAB to solve 313 00:16:48,320 --> 00:16:50,920 such a prediction problem. 314 00:16:50,920 --> 00:16:58,737 So we are thinking if the motor fires a small 3 kilogram and 10 315 00:16:58,737 --> 00:17:02,760 centimeters in diameter spherical cannonball, 316 00:17:02,760 --> 00:17:06,030 so it is fired over here. 317 00:17:06,030 --> 00:17:07,690 At sea level in standard atmosphere, 318 00:17:07,690 --> 00:17:16,119 you can look for the air density or whatever from the internet. 319 00:17:16,119 --> 00:17:18,359 And so I want to derive the differential equation. 320 00:17:18,359 --> 00:17:24,230 That is going to allow me to predict the impact if I give 321 00:17:24,230 --> 00:17:27,369 you the initial velocity, the x velocity and y 322 00:17:27,369 --> 00:17:33,060 velocity of this thing. 323 00:17:33,060 --> 00:17:36,570 Here is some additional data. 324 00:17:36,570 --> 00:17:39,290 Aerodynamic drag has to be considered. 325 00:17:39,290 --> 00:17:43,440 The cannon ball, the force on the cannonball 326 00:17:43,440 --> 00:17:46,460 is gravity and drag. 327 00:17:46,460 --> 00:17:46,980 Right? 328 00:17:46,980 --> 00:17:49,890 Those are the only two forces. 329 00:17:49,890 --> 00:17:54,740 And let's assume it's fired at a subsonic speed. 330 00:17:54,740 --> 00:17:57,000 Subsonic Cd, that coefficient. 331 00:17:57,000 --> 00:18:00,100 Let's use 0.5. 332 00:18:00,100 --> 00:18:03,186 And yeah, actually we give you the air density. 333 00:18:03,186 --> 00:18:05,680 So you don't have to find it out. 334 00:18:08,960 --> 00:18:10,935 So what's the differential equation? 335 00:18:10,935 --> 00:18:14,072 The differential equation such that you can [INAUDIBLE] 336 00:18:14,072 --> 00:18:17,050 in MATLAB. 337 00:18:17,050 --> 00:18:18,630 For deriving the equation, you can 338 00:18:18,630 --> 00:18:21,760 work in groups of two or three. 339 00:18:21,760 --> 00:18:24,992 Writing the code has to be done individually. 340 00:18:24,992 --> 00:18:26,470 All right? 341 00:18:26,470 --> 00:18:29,955 So if you're sitting alone, please tend to gather. 342 00:18:29,955 --> 00:18:30,881 Form groups. 343 00:18:30,881 --> 00:18:33,440 And continue working out the differential equations 344 00:18:33,440 --> 00:18:36,310 together. 345 00:18:36,310 --> 00:18:41,630 Any questions before we dive into the lab? 346 00:18:41,630 --> 00:18:44,055 AUDIENCE: [INAUDIBLE]. 347 00:18:44,055 --> 00:18:45,995 QIQI WANG: Uh huh. 348 00:18:45,995 --> 00:18:46,970 AUDIENCE: [INAUDIBLE] 349 00:18:46,970 --> 00:18:49,200 QIQI WANG: What is the angle? 350 00:18:49,200 --> 00:18:50,880 This is [INAUDIBLE]. 351 00:18:50,880 --> 00:18:53,280 And I'm going to tell you that. 352 00:18:53,280 --> 00:18:56,150 Let's leave it as the parameter [INAUDIBLE]. 353 00:18:56,150 --> 00:19:01,970 And I'm going to tell you where the [INAUDIBLE] what happens 354 00:19:01,970 --> 00:19:03,090 [INAUDIBLE]. 355 00:19:03,090 --> 00:19:03,590 OK. 356 00:19:03,590 --> 00:19:06,590 Lot of you have figured out the ODE. 357 00:19:06,590 --> 00:19:08,500 So I'm just going to write it down 358 00:19:08,500 --> 00:19:11,700 here so that everybody is on the same page when 359 00:19:11,700 --> 00:19:13,990 you go into implementation. 360 00:19:13,990 --> 00:19:14,490 OK. 361 00:19:18,700 --> 00:19:22,510 What happens is that the force-- if you look at the cannon ball, 362 00:19:22,510 --> 00:19:31,130 if the velocity is Ux and Uy, it has 363 00:19:31,130 --> 00:19:33,820 a force of gravity, m times g. 364 00:19:33,820 --> 00:19:38,410 It has a force of drag, right? 365 00:19:38,410 --> 00:19:42,300 And the magnitude of the drag is equal to Cv 366 00:19:42,300 --> 00:19:45,290 times half of rho U squared. 367 00:19:45,290 --> 00:19:48,190 U squared is Ux squared plus Uy squared. 368 00:19:51,380 --> 00:19:54,200 And the drag also has an angle. 369 00:19:54,200 --> 00:19:59,750 The angle is proportional to the angle of Ux and Uy 370 00:19:59,750 --> 00:20:02,210 but in the opposite direction. 371 00:20:02,210 --> 00:20:05,300 So if you write down the differential equations, 372 00:20:05,300 --> 00:20:05,880 it's dUx/dt. 373 00:20:11,150 --> 00:20:14,110 The force on the x component only 374 00:20:14,110 --> 00:20:16,060 has a component of drag on it. 375 00:20:16,060 --> 00:20:16,990 Right? 376 00:20:16,990 --> 00:20:24,600 So it is minus D times Ux divided by square root 377 00:20:24,600 --> 00:20:27,510 of Ux squared plus Uy squared. 378 00:20:27,510 --> 00:20:28,780 This is the angle. 379 00:20:28,780 --> 00:20:34,560 And as you plug it in, it is-- oh, and times the area. 380 00:20:34,560 --> 00:20:38,480 When you plug it in, it is equal to minus 1/2 381 00:20:38,480 --> 00:20:48,320 of Cd times row S times square root of Ux squared plus Uy 382 00:20:48,320 --> 00:20:50,890 squared times Ux. 383 00:20:50,890 --> 00:21:06,702 And dUy dt is equal to first of all-- 384 00:21:06,702 --> 00:21:08,290 Wait. 385 00:21:08,290 --> 00:21:11,160 I think I'm missing an m here. 386 00:21:11,160 --> 00:21:12,770 m times d Uy/dt. 387 00:21:12,770 --> 00:21:15,480 That is the acceleration on the y 388 00:21:15,480 --> 00:21:18,060 direction, which is also equal to the force in the y 389 00:21:18,060 --> 00:21:18,990 direction. 390 00:21:18,990 --> 00:21:25,760 It is equal to minus m times g minus a similar drag 391 00:21:25,760 --> 00:21:30,850 component which has the same Cd Row S over two. 392 00:21:30,850 --> 00:21:37,090 It has the same square root of Ux squared plus Uy squared. 393 00:21:37,090 --> 00:21:38,940 But the direction is in the y. 394 00:21:38,940 --> 00:21:42,540 So it is proportional to the Uy. 395 00:21:42,540 --> 00:21:43,040 All right. 396 00:21:43,040 --> 00:21:45,160 So these are the two differential equations. 397 00:21:45,160 --> 00:21:48,480 And when you want to compute the trajectory, 398 00:21:48,480 --> 00:21:52,090 you also need dx/dt equals what? 399 00:21:55,250 --> 00:22:00,030 x and y is the location of the cannon ball. 400 00:22:00,030 --> 00:22:04,845 Ux have and dy/dt equal to Uy. 401 00:22:04,845 --> 00:22:08,170 So these are the four differential equations 402 00:22:08,170 --> 00:22:11,490 you need to implement and solve. 403 00:22:11,490 --> 00:22:13,648 So let's try fourth order first. 404 00:22:16,940 --> 00:22:21,800 And if you checked your email over the last half an hour, 405 00:22:21,800 --> 00:22:25,026 you're going to see I shared a Dropbox folder with you. 406 00:22:25,026 --> 00:22:27,050 How many people got that? 407 00:22:27,050 --> 00:22:28,440 OK. 408 00:22:28,440 --> 00:22:31,700 So if you are happy with it and if you 409 00:22:31,700 --> 00:22:33,210 have Dropbox on your computer, you 410 00:22:33,210 --> 00:22:37,080 can make the subdirectory in the shared folder 411 00:22:37,080 --> 00:22:39,860 and you can pull in that directory 412 00:22:39,860 --> 00:22:45,000 so I can take a look from here after you put it up. 413 00:22:50,310 --> 00:22:52,430 So if you look at your email, you 414 00:22:52,430 --> 00:22:59,004 are going to see a folder called 1690shared2014 415 00:22:59,004 --> 00:22:59,920 being shared with you. 416 00:22:59,920 --> 00:23:02,490 If you accept it, everybody is going 417 00:23:02,490 --> 00:23:05,060 to see the same content here. 418 00:23:05,060 --> 00:23:09,980 And I see people have already uploaded code over here. 419 00:23:14,530 --> 00:23:15,030 Great. 420 00:23:15,030 --> 00:23:16,220 Thank you very much. 421 00:23:16,220 --> 00:23:17,725 I'm going to open your model. 422 00:23:17,725 --> 00:23:18,641 AUDIENCE: [INAUDIBLE]. 423 00:23:21,530 --> 00:23:22,720 QIQI WANG: That's a what? 424 00:23:22,720 --> 00:23:23,680 AUDIENCE: Let me upload my [INAUDIBLE]. 425 00:23:23,680 --> 00:23:24,190 QIQI WANG: Oh OK. 426 00:23:24,190 --> 00:23:24,690 Sure. 427 00:23:24,690 --> 00:23:25,270 Please. 428 00:23:25,270 --> 00:23:26,590 AUDIENCE: [INAUDIBLE]. 429 00:23:26,590 --> 00:23:27,140 QIQI WANG: Are you going to also-- 430 00:23:27,140 --> 00:23:27,681 AUDIENCE: OK. 431 00:23:30,550 --> 00:23:33,440 QIQI WANG: Oops. 432 00:23:33,440 --> 00:23:33,940 OK. 433 00:23:33,940 --> 00:23:34,740 I'm going to wait. 434 00:23:44,350 --> 00:23:44,940 Oh OK. 435 00:23:44,940 --> 00:23:48,320 So somebody must have not had Dropbox before. 436 00:23:48,320 --> 00:23:52,268 Because if I shared with you and you actually-- 437 00:23:52,268 --> 00:23:53,774 AUDIENCE: [INAUDIBLE]. 438 00:23:53,774 --> 00:23:54,440 QIQI WANG: Yeah. 439 00:23:54,440 --> 00:23:58,444 [LAUGHTER] 440 00:23:58,444 --> 00:23:59,360 AUDIENCE: [INAUDIBLE]. 441 00:24:06,684 --> 00:24:08,350 QIQI WANG: I get more space if I shared. 442 00:24:08,350 --> 00:24:13,150 And as a result of my sharing, somebody who didn't use Dropbox 443 00:24:13,150 --> 00:24:14,380 before now uses Dropbox. 444 00:24:14,380 --> 00:24:15,674 AUDIENCE: Oh. 445 00:24:15,674 --> 00:24:18,596 QIQI WANG: Yes. 446 00:24:18,596 --> 00:24:20,544 AUDIENCE: So you don't get access to anything? 447 00:24:27,362 --> 00:24:29,190 QIQI WANG: All right. 448 00:24:29,190 --> 00:24:29,690 OK. 449 00:24:29,690 --> 00:24:34,650 So I have codes here. 450 00:24:34,650 --> 00:24:39,804 And let me first run it to see if it also runs in my computer. 451 00:24:39,804 --> 00:24:41,130 Mmm. 452 00:24:41,130 --> 00:24:45,587 AUDIENCE: What is [INAUDIBLE]? 453 00:24:45,587 --> 00:24:46,170 QIQI WANG: Oh. 454 00:24:46,170 --> 00:24:46,940 I see, I see. 455 00:24:46,940 --> 00:24:49,760 OK. 456 00:24:49,760 --> 00:24:55,490 So let's go through the [INAUDIBLE] 457 00:24:55,490 --> 00:25:01,684 here and explain to us what is your code doing. 458 00:25:01,684 --> 00:25:02,350 AUDIENCE: Right. 459 00:25:09,796 --> 00:25:11,170 So first, I put in some constants 460 00:25:11,170 --> 00:25:15,130 that we need-- the diameter, the coefficient drag. 461 00:25:15,130 --> 00:25:16,440 I apologize for any errors. 462 00:25:16,440 --> 00:25:17,954 I'm definitely not immune to errors. 463 00:25:17,954 --> 00:25:19,120 I don't think there are any. 464 00:25:19,120 --> 00:25:20,730 But we calculate drag. 465 00:25:20,730 --> 00:25:23,340 And then I calculate the cross-sectional area, 466 00:25:23,340 --> 00:25:24,955 rho, gravity. 467 00:25:24,955 --> 00:25:26,830 And then I figure out the initial conditions. 468 00:25:26,830 --> 00:25:29,740 So I didn't actually hear if you gave us an output. 469 00:25:29,740 --> 00:25:31,111 So I just put in pi over 4. 470 00:25:31,111 --> 00:25:31,570 QIQI WANG: OK great. 471 00:25:31,570 --> 00:25:33,710 AUDIENCE: But I could have put anything in there. 472 00:25:33,710 --> 00:25:38,060 So I calculated the initial x and y velocities here. 473 00:25:38,060 --> 00:25:40,430 And then I create vectors with plenty of space 474 00:25:40,430 --> 00:25:44,244 to hold all my [INAUDIBLE] and set the initial conditions 475 00:25:44,244 --> 00:25:44,910 in the velocity. 476 00:25:44,910 --> 00:25:46,160 And then I create a time step. 477 00:25:46,160 --> 00:25:47,940 And then this is the integration. 478 00:25:47,940 --> 00:25:52,430 So what I say is the position at 1 times 7, 479 00:25:52,430 --> 00:25:56,100 the future x position is the current x position 480 00:25:56,100 --> 00:25:59,680 times et times the x velocity. 481 00:25:59,680 --> 00:26:01,040 And same thing for y. 482 00:26:01,040 --> 00:26:03,520 And I have to calculate the changes to the velocity of both 483 00:26:03,520 --> 00:26:04,450 x and y. 484 00:26:04,450 --> 00:26:07,795 And then I add those, again multiplying by the time step. 485 00:26:07,795 --> 00:26:11,962 And then here it hits the ground again, I think. 486 00:26:11,962 --> 00:26:14,532 And then I plot. 487 00:26:14,532 --> 00:26:15,240 QIQI WANG: Great. 488 00:26:15,240 --> 00:26:17,112 Thank you. 489 00:26:17,112 --> 00:26:18,788 Any questions about this [INAUDIBLE]? 490 00:26:23,468 --> 00:26:25,060 Is it perfectly clear, everything? 491 00:26:25,060 --> 00:26:25,976 AUDIENCE: [INAUDIBLE]. 492 00:26:27,704 --> 00:26:28,370 QIQI WANG: Here. 493 00:26:28,370 --> 00:26:33,420 Let me run it again to see if it actually still works. 494 00:26:33,420 --> 00:26:35,956 I am going to click Run. 495 00:26:35,956 --> 00:26:36,910 It works. 496 00:26:36,910 --> 00:26:39,160 AUDIENCE: And there's some residuals here [INAUDIBLE]. 497 00:26:48,280 --> 00:26:51,880 QIQI WANG: So thank you. 498 00:26:51,880 --> 00:26:55,630 And let me take a look at another code that 499 00:26:55,630 --> 00:27:00,310 is by same thing. 500 00:27:00,310 --> 00:27:01,560 Did you update? 501 00:27:01,560 --> 00:27:03,620 AUDIENCE: Yeah. 502 00:27:03,620 --> 00:27:05,480 QIQI WANG: Let me close it. 503 00:27:05,480 --> 00:27:05,980 So-- 504 00:27:10,780 --> 00:27:12,680 AUDIENCE: Probably ballistic. 505 00:27:12,680 --> 00:27:14,891 QIQI WANG: Ballistic. 506 00:27:14,891 --> 00:27:15,516 AUDIENCE: Yeah. 507 00:27:15,516 --> 00:27:17,187 That's been upgraded. 508 00:27:17,187 --> 00:27:17,770 QIQI WANG: OK. 509 00:27:21,562 --> 00:27:22,510 AUDIENCE: Try again. 510 00:27:25,360 --> 00:27:26,100 QIQI WANG: No. 511 00:27:26,100 --> 00:27:26,600 OK. 512 00:27:26,600 --> 00:27:32,090 I got another code from-- all right. 513 00:27:32,090 --> 00:27:32,990 Thanks. 514 00:27:32,990 --> 00:27:35,506 And which one should I look at first? 515 00:27:35,506 --> 00:27:37,721 AUDIENCE: [INAUDIBLE]. 516 00:27:37,721 --> 00:27:40,000 QIQI WANG: [INAUDIBLE]. 517 00:27:40,000 --> 00:27:42,680 Do you mind coming here and explaining what you did? 518 00:27:42,680 --> 00:27:49,490 I think it's a little bit different from what-- 519 00:27:49,490 --> 00:27:52,500 so you organized the code in a different way. 520 00:27:52,500 --> 00:27:56,440 I think this is more consistent with what 521 00:27:56,440 --> 00:27:58,675 I did in the last section. 522 00:27:58,675 --> 00:27:59,300 AUDIENCE: Yeah. 523 00:28:02,754 --> 00:28:03,920 QIQI WANG: Of the other one. 524 00:28:03,920 --> 00:28:06,560 Yeah, the other one. 525 00:28:06,560 --> 00:28:11,968 AUDIENCE: So what I did was to find all my initial values, 526 00:28:11,968 --> 00:28:16,234 I chose a state of 45 degrees. 527 00:28:16,234 --> 00:28:19,071 And this is all [INAUDIBLE], actually. 528 00:28:19,071 --> 00:28:21,276 [LAUGHTER] 529 00:28:21,276 --> 00:28:26,320 I defined my rate on each side as 1,000 elements. 530 00:28:26,320 --> 00:28:29,454 And then I chose a dt of 0.01. 531 00:28:29,454 --> 00:28:31,340 And then I made an initial vector 532 00:28:31,340 --> 00:28:37,470 which has my x0, y0, and the dx0 and dy0. 533 00:28:37,470 --> 00:28:39,170 And then we go into the while loop 534 00:28:39,170 --> 00:28:42,440 to do our forward Euler method. 535 00:28:42,440 --> 00:28:52,890 And then to calculate our next element in our forward Euler 536 00:28:52,890 --> 00:28:57,490 method, I go to this function called trajectory that I made. 537 00:28:57,490 --> 00:28:59,900 QIQI WANG: So trajectory t is a function you made. 538 00:28:59,900 --> 00:29:01,250 So let me open this also. 539 00:29:03,810 --> 00:29:04,495 All right. 540 00:29:04,495 --> 00:29:05,120 AUDIENCE: Yeah. 541 00:29:05,120 --> 00:29:08,610 And then here, you can see the equation [INAUDIBLE]. 542 00:29:11,140 --> 00:29:13,016 And so what this trajectory does is 543 00:29:13,016 --> 00:29:18,630 it takes in the position and velocity. 544 00:29:18,630 --> 00:29:21,005 It applies them to the [INAUDIBLE] 545 00:29:21,005 --> 00:29:27,490 we derived, and then outputs two velocities and two [INAUDIBLE], 546 00:29:27,490 --> 00:29:30,855 which we then multiply by a change in our time step 547 00:29:30,855 --> 00:29:35,322 to solve for the new position 548 00:29:35,322 --> 00:29:37,030 QIQI WANG: We have any questions on this? 549 00:29:39,970 --> 00:29:46,850 And in the other file, I believe I need to change this to this. 550 00:29:46,850 --> 00:29:48,600 Because I think you changed the name 551 00:29:48,600 --> 00:29:52,534 of this which you also have to change this too. 552 00:29:52,534 --> 00:29:54,904 So that way, [INAUDIBLE]. 553 00:29:54,904 --> 00:29:58,130 So [INAUDIBLE], you can do the derivative that you did. 554 00:29:58,130 --> 00:30:00,230 AUDIENCE: After we did the derivative, 555 00:30:00,230 --> 00:30:04,618 we have to multiply it by our time step 556 00:30:04,618 --> 00:30:08,380 to go up an order to change it from velocity to speed 557 00:30:08,380 --> 00:30:10,100 and position. 558 00:30:10,100 --> 00:30:13,792 We add that to our old position and get our new position. 559 00:30:13,792 --> 00:30:17,410 And then we store that in our position vectors. 560 00:30:17,410 --> 00:30:20,530 And then we [INAUDIBLE] the position. 561 00:30:20,530 --> 00:30:21,580 QIQI WANG: OK. 562 00:30:21,580 --> 00:30:23,560 Let's try to see if it works. 563 00:30:29,480 --> 00:30:29,980 OK. 564 00:30:29,980 --> 00:30:30,980 Let me close everything. 565 00:30:33,562 --> 00:30:36,420 Let's click Run. 566 00:30:36,420 --> 00:30:36,920 Yeah. 567 00:30:36,920 --> 00:30:40,251 We get a [INAUDIBLE]. 568 00:30:40,251 --> 00:30:42,860 AUDIENCE: [INAUDIBLE] 569 00:30:42,860 --> 00:30:43,760 QIQI WANG: Yeah. 570 00:30:43,760 --> 00:30:44,620 [LAUGHTER] 571 00:30:44,620 --> 00:30:50,190 And we do see that [INAUDIBLE], right? 572 00:30:50,190 --> 00:30:57,885 Because the angle here is much deeper than the other angles. 573 00:30:57,885 --> 00:31:03,340 [INAUDIBLE] slower than when it was shut off. 574 00:31:03,340 --> 00:31:03,940 OK. 575 00:31:03,940 --> 00:31:07,760 Now my question to you is, what if I want to challenge you 576 00:31:07,760 --> 00:31:09,380 with a question? 577 00:31:09,380 --> 00:31:12,170 What if I want to change this to midpoint? 578 00:31:12,170 --> 00:31:15,390 From last lecture, we saw that at the same time step, 579 00:31:15,390 --> 00:31:19,860 midpoint was much more accurate than forward Euler. 580 00:31:19,860 --> 00:31:21,180 Right? 581 00:31:21,180 --> 00:31:25,920 So what should we do to change this to midpoint? 582 00:31:25,920 --> 00:31:27,160 I'm asking the whole class. 583 00:31:34,315 --> 00:31:35,810 Anybody have any idea? 584 00:31:40,730 --> 00:31:42,054 Yes? 585 00:31:42,054 --> 00:31:42,970 AUDIENCE: [INAUDIBLE]. 586 00:32:00,010 --> 00:32:00,773 QIQI WANG: Good. 587 00:32:00,773 --> 00:32:03,671 AUDIENCE: [INAUDIBLE]. 588 00:32:03,671 --> 00:32:04,637 QIQI WANG: Thank you. 589 00:32:04,637 --> 00:32:07,052 AUDIENCE: [INAUDIBLE]. 590 00:32:07,052 --> 00:32:08,501 QIQI WANG: [INAUDIBLE]. 591 00:32:08,501 --> 00:32:10,440 OK. 592 00:32:10,440 --> 00:32:11,450 OK. 593 00:32:11,450 --> 00:32:13,360 So, right. 594 00:32:13,360 --> 00:32:15,140 Now the question is what if we want 595 00:32:15,140 --> 00:32:18,210 to change this to midpoint? 596 00:32:18,210 --> 00:32:21,650 If we change this to midpoint, the first thing 597 00:32:21,650 --> 00:32:27,890 is that midpoint is a how many step scheme? 598 00:32:27,890 --> 00:32:29,140 It's a two-step scheme. 599 00:32:29,140 --> 00:32:31,270 Forward Euler is a one step scheme. 600 00:32:31,270 --> 00:32:35,515 So in terms of implementation, what do we need to change? 601 00:32:35,515 --> 00:32:41,230 What do we need to ignore when we switch from a one step 602 00:32:41,230 --> 00:32:45,700 scheme to a two-step scheme? 603 00:32:45,700 --> 00:32:49,585 Why do we call it a two step scheme? 604 00:32:49,585 --> 00:32:52,040 AUDIENCE: [INAUDIBLE]. 605 00:32:52,040 --> 00:32:56,900 QIQI WANG: Because we need to store two time steps. 606 00:32:56,900 --> 00:33:00,580 Here we only need to have that. 607 00:33:00,580 --> 00:33:04,170 That is one previous time step. 608 00:33:04,170 --> 00:33:08,370 Now we have to have two previous time steps. 609 00:33:08,370 --> 00:33:09,800 Right? 610 00:33:09,800 --> 00:33:12,950 So we need this. 611 00:33:12,950 --> 00:33:15,470 And we need one step of forward Euler 612 00:33:15,470 --> 00:33:19,010 to actually get to the next time step. 613 00:33:19,010 --> 00:33:27,111 So I'm just going to say a one step of forward Euler. 614 00:33:27,111 --> 00:33:27,610 OK. 615 00:33:27,610 --> 00:33:29,790 So in doing one step of forward Euler, 616 00:33:29,790 --> 00:33:34,780 I'm just going to copy this and copy this. 617 00:33:34,780 --> 00:33:35,280 Right? 618 00:33:39,210 --> 00:33:42,780 Here, I'm just going to pass my vect 0. 619 00:33:42,780 --> 00:33:47,920 And the vect equal to vect 0 plus Vt times d vect dt. 620 00:33:47,920 --> 00:33:50,220 Any questions on this one step of forward Euler? 621 00:33:55,320 --> 00:33:55,820 OK. 622 00:33:55,820 --> 00:33:59,850 Now we are using midpoint. 623 00:33:59,850 --> 00:34:01,300 OK. 624 00:34:01,300 --> 00:34:08,092 In mid-point, now what we want to change is these two lines. 625 00:34:08,092 --> 00:34:08,909 Right? 626 00:34:08,909 --> 00:34:11,580 Let me scroll it down. 627 00:34:11,580 --> 00:34:13,170 We need to change these two lines. 628 00:34:13,170 --> 00:34:15,569 Can somebody tell me how do I change lines 629 00:34:15,569 --> 00:34:18,200 when [INAUDIBLE] to make this midpoint instead 630 00:34:18,200 --> 00:34:19,311 of forward Euler? 631 00:34:27,159 --> 00:34:27,867 Yeah? 632 00:34:27,867 --> 00:34:30,157 AUDIENCE: You have the vect. 633 00:34:30,157 --> 00:34:38,826 And you add it to 2 times U vect dt and your time step. 634 00:34:38,826 --> 00:34:39,409 QIQI WANG: OK. 635 00:34:39,409 --> 00:34:43,949 So first of all, in midpoint, we have 2 times 636 00:34:43,949 --> 00:34:48,070 delta t times F of v. Right? 637 00:34:48,070 --> 00:34:50,024 So we have a 2 times. 638 00:34:50,024 --> 00:34:51,315 What else do we need to change? 639 00:34:53,910 --> 00:34:56,699 AUDIENCE: [INAUDIBLE]. 640 00:34:56,699 --> 00:34:58,372 QIQI WANG: For which one? 641 00:35:01,236 --> 00:35:02,611 AUDIENCE: In the previous scheme? 642 00:35:06,040 --> 00:35:07,030 QIQI WANG: Yeah. 643 00:35:07,030 --> 00:35:08,640 You are right. 644 00:35:08,640 --> 00:35:21,130 When I do this, I need to also somehow store the old vect. 645 00:35:21,130 --> 00:35:28,210 So I need to do vect of 0. 646 00:35:28,210 --> 00:35:30,850 Let me actually do this. 647 00:35:30,850 --> 00:35:36,830 Let me actually call this vect 00 as the previous time step. 648 00:35:36,830 --> 00:35:43,650 And the vect 0 as-- so, at any point, that 00? 649 00:35:46,190 --> 00:35:49,890 Let me call it k minus 1. 650 00:35:49,890 --> 00:35:50,520 OK. 651 00:35:50,520 --> 00:35:51,020 vect 0. 652 00:35:51,020 --> 00:35:52,470 Let me call it vect k. 653 00:36:02,290 --> 00:36:03,530 All right. 654 00:36:03,530 --> 00:36:08,240 So this is one step of forward Euler going from k minus 1 655 00:36:08,240 --> 00:36:10,410 to k. 656 00:36:10,410 --> 00:36:11,220 All right. 657 00:36:11,220 --> 00:36:18,456 Now in midpoint, where do I compute the derivative? 658 00:36:18,456 --> 00:36:19,490 x time step k? 659 00:36:19,490 --> 00:36:20,520 Or k minus 1? 660 00:36:23,748 --> 00:36:25,200 K, right. 661 00:36:25,200 --> 00:36:27,090 So that's k. 662 00:36:27,090 --> 00:36:32,710 So this d vect/dt is the time derivative at step k. 663 00:36:32,710 --> 00:36:41,750 Now my vect is equal to vect k minus 1 plus 2 times 664 00:36:41,750 --> 00:36:47,590 delta t times F of v. Right? 665 00:36:47,590 --> 00:36:58,080 This is because in midpoint-- so forward Euler is V at k plus 1 666 00:36:58,080 --> 00:37:02,540 equal to Vk plus f of Vk. 667 00:37:02,540 --> 00:37:10,110 Midpoint, we have Vk plus 1 equal to be Vk minus 1 plus 2. 668 00:37:10,110 --> 00:37:11,410 There is a delta t here. 669 00:37:11,410 --> 00:37:14,761 2 times f Vk times delta t. 670 00:37:14,761 --> 00:37:15,260 Right? 671 00:37:15,260 --> 00:37:20,280 This is what we are implementing right now. 672 00:37:20,280 --> 00:37:21,100 All right. 673 00:37:21,100 --> 00:37:24,770 So we have this. 674 00:37:24,770 --> 00:37:28,500 And everything else I think is same. 675 00:37:30,664 --> 00:37:31,580 AUDIENCE: [INAUDIBLE]. 676 00:37:36,852 --> 00:37:37,560 QIQI WANG: Right. 677 00:37:37,560 --> 00:37:39,700 Then I need to also update. 678 00:37:39,700 --> 00:37:44,820 Because I need to set vect of k equal to vect. 679 00:37:44,820 --> 00:37:46,210 Right? 680 00:37:46,210 --> 00:37:52,475 And before I do that, I need to have vect k minus 1 681 00:37:52,475 --> 00:37:56,714 is equal to vect k. 682 00:37:56,714 --> 00:37:59,650 Right? 683 00:37:59,650 --> 00:38:03,119 So I'm done. 684 00:38:03,119 --> 00:38:03,910 Should we run this? 685 00:38:06,790 --> 00:38:09,520 Is it clear like if when you guys have 686 00:38:09,520 --> 00:38:11,570 fourth order implemented on your computer 687 00:38:11,570 --> 00:38:15,680 how to change your implementation into midpoint? 688 00:38:15,680 --> 00:38:17,010 Good? 689 00:38:17,010 --> 00:38:18,030 Let's run it. 690 00:38:23,160 --> 00:38:24,354 Let me close out. 691 00:38:27,030 --> 00:38:28,030 Now I'm going to run it. 692 00:38:34,862 --> 00:38:38,787 Is it the same as before? 693 00:38:38,787 --> 00:38:39,338 No. 694 00:38:39,338 --> 00:38:40,254 AUDIENCE: [INAUDIBLE]. 695 00:38:43,680 --> 00:38:44,547 QIQI WANG: Oh. 696 00:38:44,547 --> 00:38:46,746 AUDIENCE: [INAUDIBLE]. 697 00:38:46,746 --> 00:38:47,329 QIQI WANG: Oh. 698 00:38:47,329 --> 00:38:48,220 Is this something to-- 699 00:38:48,220 --> 00:38:49,270 AUDIENCE: Yeah. [INAUDIBLE]. 700 00:38:49,270 --> 00:38:50,120 QIQI WANG: You mean remotely? 701 00:38:50,120 --> 00:38:51,286 AUDIENCE: Yeah. [INAUDIBLE]. 702 00:38:51,286 --> 00:38:56,367 [LAUGHTER] 703 00:38:56,367 --> 00:38:56,950 QIQI WANG: OK. 704 00:38:56,950 --> 00:38:57,449 OK. 705 00:38:57,449 --> 00:38:59,680 What did you do? 706 00:38:59,680 --> 00:39:02,917 AUDIENCE: I didn't [INAUDIBLE]. 707 00:39:02,917 --> 00:39:03,500 QIQI WANG: Oh. 708 00:39:03,500 --> 00:39:04,660 OK. 709 00:39:04,660 --> 00:39:05,180 OK. 710 00:39:05,180 --> 00:39:11,074 So let me just make it 10 more time steps. 711 00:39:11,074 --> 00:39:12,977 AUDIENCE: [INAUDIBLE]. 712 00:39:12,977 --> 00:39:13,810 QIQI WANG: Oh, yeah. 713 00:39:13,810 --> 00:39:14,310 Right. 714 00:39:18,352 --> 00:39:19,268 AUDIENCE: [INAUDIBLE]. 715 00:39:23,677 --> 00:39:24,302 QIQI WANG: Yes. 716 00:39:27,920 --> 00:39:28,640 OK. 717 00:39:28,640 --> 00:39:29,410 Let me run it. 718 00:39:32,494 --> 00:39:33,410 AUDIENCE: [INAUDIBLE]. 719 00:39:38,910 --> 00:39:43,180 QIQI WANG: Let me-- I think this is too much. 720 00:39:43,180 --> 00:39:50,780 Let me change it to-- I think I just do this much. 721 00:39:50,780 --> 00:39:52,120 It may be enough. 722 00:39:52,120 --> 00:39:53,700 Let me see it. 723 00:39:53,700 --> 00:39:58,580 And when you plot it, let me make sure to do x is equal. 724 00:39:58,580 --> 00:39:59,970 Do you know what this means? 725 00:39:59,970 --> 00:40:02,120 This just means that the x-axis and y-axis 726 00:40:02,120 --> 00:40:04,560 are going to be on scale. 727 00:40:04,560 --> 00:40:05,430 Right? 728 00:40:05,430 --> 00:40:06,500 OK. 729 00:40:06,500 --> 00:40:07,340 That's good. 730 00:40:07,340 --> 00:40:10,360 So midpoint also works. 731 00:40:10,360 --> 00:40:16,440 So first of all, I want to use the remaining few minutes, 732 00:40:16,440 --> 00:40:19,700 let's try to implement the most accurate scheme 733 00:40:19,700 --> 00:40:24,830 we got to see-- we already have a two-step scheme. 734 00:40:24,830 --> 00:40:26,100 Right? 735 00:40:26,100 --> 00:40:30,480 And it should be pretty easy to change to this scheme 736 00:40:30,480 --> 00:40:33,330 we have over here. 737 00:40:33,330 --> 00:40:43,100 So instead of here, so we have a d vect dt at step k, right? 738 00:40:43,100 --> 00:40:47,850 We also need a d vect dt at step k minus 1. 739 00:40:47,850 --> 00:40:53,690 That is just called in the same function at step k minus 1. 740 00:40:53,690 --> 00:40:56,130 Right? 741 00:40:56,130 --> 00:40:57,110 You see what I'm doing? 742 00:40:57,110 --> 00:40:59,820 I'm computing two time derivatives. 743 00:40:59,820 --> 00:41:07,470 Now in this update, I have a minus 4 times Vk 744 00:41:07,470 --> 00:41:17,780 and the plus 5 times Vk minus 1 plus 5 times this, right? 745 00:41:17,780 --> 00:41:18,280 OK. 746 00:41:18,280 --> 00:41:26,450 And then the second line is plus 4 delta t f of Vk 747 00:41:26,450 --> 00:41:30,590 plus 2 delta t f of Vk minus 1. 748 00:41:30,590 --> 00:41:42,485 So this is 4 plus 2 times dt times d vect dt k minus one. 749 00:41:42,485 --> 00:41:42,985 Right? 750 00:41:45,600 --> 00:41:50,830 So that completes our change from midpoint 751 00:41:50,830 --> 00:41:55,270 to an even more accurate two-step scheme. 752 00:41:55,270 --> 00:41:58,210 And by more accurate, I mean global or local order 753 00:41:58,210 --> 00:42:00,744 of accuracy? 754 00:42:00,744 --> 00:42:01,410 AUDIENCE: Local. 755 00:42:01,410 --> 00:42:02,409 QIQI WANG: Local, right? 756 00:42:02,409 --> 00:42:06,950 Because we did this using Taylor series analysis. 757 00:42:06,950 --> 00:42:11,490 So let's see if it translates into global order of accuracy. 758 00:42:11,490 --> 00:42:15,974 Let me click, let me close this first. 759 00:42:15,974 --> 00:42:17,190 And let's run it. 760 00:42:22,547 --> 00:42:25,469 You see this? 761 00:42:25,469 --> 00:42:27,904 You see this? 762 00:42:27,904 --> 00:42:31,800 What does this mean, 10 to the 162? 763 00:42:31,800 --> 00:42:34,048 That's where my cannonball bounced. 764 00:42:38,016 --> 00:42:38,720 What happens? 765 00:42:43,440 --> 00:42:44,007 What happens? 766 00:42:47,010 --> 00:42:49,710 Let's look at our scheme solution 767 00:42:49,710 --> 00:42:51,430 not for this many time steps. 768 00:42:51,430 --> 00:42:54,250 But let's just go 19 time steps. 769 00:42:54,250 --> 00:42:57,400 And let's put a circle at every time step 770 00:42:57,400 --> 00:42:58,777 to see what actually happens. 771 00:43:01,365 --> 00:43:01,864 OK. 772 00:43:06,040 --> 00:43:10,790 After just 19 time steps, we've got 10 to the 22. 773 00:43:10,790 --> 00:43:13,290 That's amazing. 774 00:43:13,290 --> 00:43:20,500 And we basically just shoot out and diverge, right? 775 00:43:20,500 --> 00:43:23,210 So why do you think that is the case? 776 00:43:23,210 --> 00:43:25,450 Do we have any global order of accuracy? 777 00:43:28,150 --> 00:43:31,185 I mean, we can do that by changing the time 778 00:43:31,185 --> 00:43:33,940 steps to be even smaller. 779 00:43:33,940 --> 00:43:36,530 And we would change-- don't change it to this. 780 00:43:36,530 --> 00:43:39,450 Let me change it to twice smaller. 781 00:43:39,450 --> 00:43:44,620 And the time steps to also twice as much. 782 00:43:44,620 --> 00:43:47,440 That is a good way of assessing the global order of accuracy, 783 00:43:47,440 --> 00:43:48,270 right? 784 00:43:48,270 --> 00:43:51,580 It will decrease the time step by a factor of two, 785 00:43:51,580 --> 00:43:54,240 increase the number of time steps by a factor of two, 786 00:43:54,240 --> 00:43:55,226 right? 787 00:43:55,226 --> 00:43:59,880 To see if the solution gets more accurate. 788 00:43:59,880 --> 00:44:01,060 OK. 789 00:44:01,060 --> 00:44:02,537 It gets wilder. 790 00:44:02,537 --> 00:44:04,120 Instead of 10 to the 100 or something, 791 00:44:04,120 --> 00:44:06,760 it's 10 to the 200 and something. 792 00:44:06,760 --> 00:44:07,972 All right. 793 00:44:07,972 --> 00:44:10,400 So something is wrong. 794 00:44:10,400 --> 00:44:11,000 What is wrong? 795 00:44:11,000 --> 00:44:12,600 Can somebody just shout out? 796 00:44:16,512 --> 00:44:18,300 What we reviewed in the beginning? 797 00:44:22,620 --> 00:44:25,490 Our scheme is not zero stable. 798 00:44:25,490 --> 00:44:31,290 It can't even solve this differential equation. 799 00:44:31,290 --> 00:44:32,210 Let me repeat. 800 00:44:32,210 --> 00:44:37,360 The most accurate two-step scheme can't even 801 00:44:37,360 --> 00:44:41,160 solve du/dt equal to 0. 802 00:44:41,160 --> 00:44:42,660 You want me to prove that? 803 00:44:42,660 --> 00:44:46,110 I'm going to prove that analytically. 804 00:44:46,110 --> 00:44:46,810 OK. 805 00:44:46,810 --> 00:44:48,740 We have the following scheme. 806 00:44:48,740 --> 00:44:51,770 I'll see if I can do it in five minutes. 807 00:44:51,770 --> 00:44:57,380 It's equal to minus 4 Vk plus 5 Vk minus 1. 808 00:44:57,380 --> 00:44:59,580 And the rest of them on the right hand side 809 00:44:59,580 --> 00:45:01,470 are f of V, right? 810 00:45:01,470 --> 00:45:06,580 If I solve the differential equation du/dt equal to f of U 811 00:45:06,580 --> 00:45:09,060 equal to 0, then f of U is equal to 0. 812 00:45:09,060 --> 00:45:11,570 I have nothing after this. 813 00:45:11,570 --> 00:45:13,820 And this is my scheme. 814 00:45:13,820 --> 00:45:18,610 I want this scheme to reproduce a constant solution. 815 00:45:18,610 --> 00:45:19,650 Right? 816 00:45:19,650 --> 00:45:23,060 The question is why can this scheme not 817 00:45:23,060 --> 00:45:24,585 reproduce a constant solution? 818 00:45:28,910 --> 00:45:30,370 The answer is like this. 819 00:45:30,370 --> 00:45:32,400 The answer is because this scheme 820 00:45:32,400 --> 00:45:35,890 is going to allow an exponentially diverging 821 00:45:35,890 --> 00:45:37,530 solution. 822 00:45:37,530 --> 00:45:38,200 OK. 823 00:45:38,200 --> 00:45:43,804 If I have a solution Vk is equal to V0 times 824 00:45:43,804 --> 00:45:47,020 a zeta-- let me just write out z because I can't 825 00:45:47,020 --> 00:45:51,110 write zeta-- times Z to the k. 826 00:45:51,110 --> 00:45:55,170 I'm going to see-- if I plug into this equation 827 00:45:55,170 --> 00:45:58,440 and the equation is satisfied and my V is something 828 00:45:58,440 --> 00:46:02,300 greater than 1-- then what does it mean? 829 00:46:02,300 --> 00:46:05,810 It means I have a very, very small V0, even 830 00:46:05,810 --> 00:46:09,520 if I have an extremely small V0, I 831 00:46:09,520 --> 00:46:12,980 can end up having a huge solution after sufficiently 832 00:46:12,980 --> 00:46:14,490 many time steps. 833 00:46:14,490 --> 00:46:15,370 Right? 834 00:46:15,370 --> 00:46:16,940 And when I plot this thing, I'm going 835 00:46:16,940 --> 00:46:21,710 to have a V of k plus 1 is V0 times V to the k 836 00:46:21,710 --> 00:46:28,180 plus 1 is equal to minus 4 times Vk V0 times V to the k. 837 00:46:28,180 --> 00:46:32,111 And sorry, this is k minus 1. 838 00:46:32,111 --> 00:46:32,610 Right? 839 00:46:32,610 --> 00:46:34,360 This is my scheme. 840 00:46:34,360 --> 00:46:38,830 Times V0 V to the k minus 1. 841 00:46:38,830 --> 00:46:40,140 I cancel this. 842 00:46:40,140 --> 00:46:40,870 Cancel this. 843 00:46:40,870 --> 00:46:41,950 Cancel the V0. 844 00:46:41,950 --> 00:46:45,900 Because they are the same on both sides. 845 00:46:45,900 --> 00:46:50,330 And though all of these have V to the k minus 1 at least-- 846 00:46:50,330 --> 00:46:53,670 some of them are k, some of them are k plus 1. 847 00:46:53,670 --> 00:46:55,060 So I remove that vector. 848 00:46:55,060 --> 00:47:00,630 I have V squared is equal to minus 4 Z plus 5. 849 00:47:00,630 --> 00:47:03,460 We get a quadratic equation. 850 00:47:03,460 --> 00:47:04,510 Right? 851 00:47:04,510 --> 00:47:08,400 This quadratic equation can have real or complex solutions. 852 00:47:08,400 --> 00:47:11,860 But let's see what it has. 853 00:47:11,860 --> 00:47:18,100 So a is equal to 1. b is equal to 4. c is equal to minus 5. 854 00:47:18,100 --> 00:47:24,970 So we get a minus b plus minus b squared, which is 16. 855 00:47:24,970 --> 00:47:30,068 Minus 4ac, which is plus 20. 856 00:47:30,068 --> 00:47:31,200 All right? 857 00:47:31,200 --> 00:47:34,050 So that's my two solutions. 858 00:47:34,050 --> 00:47:42,690 And I can see it is-- right? 859 00:47:42,690 --> 00:47:45,590 I'm basically factoring the two into these. 860 00:47:45,590 --> 00:47:48,280 And I minus 2, plus or minus 3. 861 00:47:48,280 --> 00:47:50,500 So that's equal to what? 862 00:47:50,500 --> 00:47:54,200 It's equal to either minus 5 or 1. 863 00:47:56,780 --> 00:47:57,820 All right. 864 00:47:57,820 --> 00:47:58,320 OK. 865 00:47:58,320 --> 00:48:01,480 If I only look at the one, that means good. 866 00:48:01,480 --> 00:48:05,520 It allows a constant solution, right? 867 00:48:05,520 --> 00:48:10,350 Z is equal to 1 basically means Vk is equal to V0 period, 868 00:48:10,350 --> 00:48:12,270 or whatever k. 869 00:48:12,270 --> 00:48:18,260 So in this sense, it does allow a numerical solution 870 00:48:18,260 --> 00:48:19,890 of this differential equation. 871 00:48:19,890 --> 00:48:21,830 du/dt is equal to zero. 872 00:48:21,830 --> 00:48:27,142 But it has another solution of Z equal to minus 5. 873 00:48:27,142 --> 00:48:32,610 And it is that solution that makes the scheme diverge. 874 00:48:32,610 --> 00:48:39,180 You only need one bad solution Z to make the scheme diverge. 875 00:48:39,180 --> 00:48:41,230 And that's what we saw over here. 876 00:48:41,230 --> 00:48:43,350 The scheme is not zero stable. 877 00:48:43,350 --> 00:48:45,860 Therefore, although it is locally accurate, 878 00:48:45,860 --> 00:48:48,590 it is not globally accurate. 879 00:48:48,590 --> 00:48:50,772 But because it's not zero stable, 880 00:48:50,772 --> 00:48:54,382 it can't solve this equation. 881 00:48:54,382 --> 00:48:56,717 All right? 882 00:48:56,717 --> 00:49:00,210 I'll see you tomorrow and continue discussing this 883 00:49:00,210 --> 00:49:03,280 and also look into Eigenvalue stability.