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,780 Commons license. 3 00:00:03,780 --> 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,477 from hundreds of MIT courses, visit MIT OpenCourseWare 7 00:00:16,477 --> 00:00:17,351 at [? ocw.mit.edu. ?] 8 00:00:25,854 --> 00:00:26,950 QIQI WANG: All right. 9 00:00:26,950 --> 00:00:31,360 So last time I did this using the video recording, 10 00:00:31,360 --> 00:00:33,540 there was some complaints about the audio quality, 11 00:00:33,540 --> 00:00:34,860 so I'm wearing this today. 12 00:00:34,860 --> 00:00:36,555 It feels a little bit weird, but I'll 13 00:00:36,555 --> 00:00:41,760 see how much better the audio quality came through. 14 00:00:41,760 --> 00:00:47,320 And you can ask your questions [INAUDIBLE]. 15 00:00:47,320 --> 00:00:48,646 All right. 16 00:00:48,646 --> 00:00:52,630 Can you pass this on, like everybody have one piece? 17 00:00:52,630 --> 00:00:56,860 So the purpose of this is for you to just write down 18 00:00:56,860 --> 00:01:01,760 anything that is not clear onto this piece of paper, 19 00:01:01,760 --> 00:01:04,864 and hand it to me after the lecture. 20 00:01:04,864 --> 00:01:06,280 What I'm going to do is, I'm going 21 00:01:06,280 --> 00:01:08,750 to draft these questions in the next section 22 00:01:08,750 --> 00:01:13,777 so that you're not going to stay confused 23 00:01:13,777 --> 00:01:16,723 on what might confuse you. 24 00:01:16,723 --> 00:01:19,110 So this important [? Monte ?] cuts. 25 00:01:19,110 --> 00:01:22,140 You can write anything that is [? Monte, ?] and I'm going 26 00:01:22,140 --> 00:01:27,532 to be lacking them and try to earmark your [? Monte ?] 27 00:01:27,532 --> 00:01:28,032 points. 28 00:01:30,930 --> 00:01:36,720 So we are continuing our discussion integration 29 00:01:36,720 --> 00:01:40,814 of ordinary differential equations. 30 00:01:40,814 --> 00:01:43,510 Let's just do a brief review of what 31 00:01:43,510 --> 00:01:46,410 we did in the last lecture. 32 00:01:46,410 --> 00:01:48,860 So we are solving differential equations 33 00:01:48,860 --> 00:01:53,472 like du over dt equal to f of u. 34 00:01:53,472 --> 00:01:55,410 All right. 35 00:01:55,410 --> 00:02:00,340 We discretized the time into t0 equal to 0, 36 00:02:00,340 --> 00:02:06,640 t1 equal to delta t, t2 equal to 2 delta t, et cetera. 37 00:02:06,640 --> 00:02:12,590 We discretized the solution, u, into u0 equally to u at t0, 38 00:02:12,590 --> 00:02:20,540 et cetera and ui denoted as the solution, u at ti. 39 00:02:20,540 --> 00:02:25,020 Forward Euler, solves the equation like this. 40 00:02:25,020 --> 00:02:30,180 It is basically taking-- I'm going to change the color. 41 00:02:30,180 --> 00:02:34,390 It is essentially taking this differential equation 42 00:02:34,390 --> 00:02:40,200 and approximate the time derivative into ui 43 00:02:40,200 --> 00:02:43,370 plus 1 minus ui divided by delta t. 44 00:02:43,370 --> 00:02:48,340 And on the right-hand side, it is the same thing. 45 00:02:48,340 --> 00:02:51,340 So this is forward Euler. 46 00:02:51,340 --> 00:02:55,340 A midpoint rule, we are going to be doing something 47 00:02:55,340 --> 00:03:00,920 very similar, but instead, I'm going 48 00:03:00,920 --> 00:03:03,430 to discretize the time derivative term 49 00:03:03,430 --> 00:03:04,320 in a different way. 50 00:03:04,320 --> 00:03:07,540 I'm going to discretize the time derivative term into ui 51 00:03:07,540 --> 00:03:11,090 plus 1 minus ui minus 1 divided by what? 52 00:03:11,090 --> 00:03:12,380 AUDIENCE: [INAUDIBLE]. 53 00:03:12,380 --> 00:03:14,000 QIQI WANG: 2 delta t. 54 00:03:14,000 --> 00:03:14,500 Exactly. 55 00:03:14,500 --> 00:03:17,200 So this is a consistent approximation 56 00:03:17,200 --> 00:03:18,966 of the time derivative. 57 00:03:18,966 --> 00:03:20,840 So that's the only difference between forward 58 00:03:20,840 --> 00:03:24,210 Euler and midpoint rule. 59 00:03:24,210 --> 00:03:28,540 And how to analyze the accuracy of these two methods? 60 00:03:28,540 --> 00:03:30,170 Taylor series expansion. 61 00:03:30,170 --> 00:03:35,880 In midpoint rule, we are expanding ui plus 1 equal to ui 62 00:03:35,880 --> 00:03:41,120 plus-- you should become very familiar with this 63 00:03:41,120 --> 00:03:43,070 as you learn more. 64 00:03:43,070 --> 00:03:46,600 So ui prime, which is a shorthand for denoting 65 00:03:46,600 --> 00:03:53,490 du dt at ti time delta t plus half of ui 66 00:03:53,490 --> 00:03:56,390 double prime, which means du d squared u, 67 00:03:56,390 --> 00:04:00,230 dt squared, the second derivative of u times dt 68 00:04:00,230 --> 00:04:03,620 squared plus all the other terms. 69 00:04:03,620 --> 00:04:06,710 The third term is ui triple prime delta 70 00:04:06,710 --> 00:04:11,710 t cubed plus et cetera. 71 00:04:11,710 --> 00:04:15,155 Is a little bit too wide, maybe. 72 00:04:15,155 --> 00:04:16,975 AUDIENCE: [INAUDIBLE]. 73 00:04:16,975 --> 00:04:18,079 QIQI WANG: You can't see. 74 00:04:18,079 --> 00:04:22,730 Let me change this to a little bit darker green. 75 00:04:22,730 --> 00:04:26,500 And ui minus 1 is equal to-- let me 76 00:04:26,500 --> 00:04:31,700 scroll this up a little bit-- is ui minus ui prime delta t. 77 00:04:31,700 --> 00:04:32,548 Why minus? 78 00:04:35,760 --> 00:04:37,510 I'm going backwards. 79 00:04:37,510 --> 00:04:42,640 I should test delta t to minus delta t, right? 80 00:04:42,640 --> 00:04:44,760 This is why there is a minus here. 81 00:04:44,760 --> 00:04:51,640 And this term should have also a minus here because minus delta 82 00:04:51,640 --> 00:04:53,050 t is squared. 83 00:04:53,050 --> 00:04:56,300 So I have, again, altered it here. 84 00:04:56,300 --> 00:04:58,710 And this term? 85 00:04:58,710 --> 00:05:01,370 Should I have a minus here? 86 00:05:01,370 --> 00:05:02,300 Yes. 87 00:05:02,300 --> 00:05:05,390 Because when I cube it, the negative sign 88 00:05:05,390 --> 00:05:08,100 is preserved [INAUDIBLE]. 89 00:05:08,100 --> 00:05:11,900 So when I take the subtraction between ui plus 1 90 00:05:11,900 --> 00:05:16,390 and ui minus 1, this got canceled this gets doubled, 91 00:05:16,390 --> 00:05:21,990 this gets canceled, and this gets doubled. 92 00:05:21,990 --> 00:05:34,590 So what I ended up having is ui prime plus 1/6 of ui 93 00:05:34,590 --> 00:05:36,900 triple prime delta t squared. 94 00:05:36,900 --> 00:05:39,902 I mean, this cube cancels with this delta t and delta 95 00:05:39,902 --> 00:05:44,090 t squared equal to f of ui. 96 00:05:44,090 --> 00:05:50,282 And this is n plus 0 delta t cubed, 97 00:05:50,282 --> 00:05:53,110 so this is the leading term of the truncation error. 98 00:05:56,270 --> 00:05:58,440 Remember my origin original differential equation. 99 00:06:02,890 --> 00:06:04,690 My origin or differential equation 100 00:06:04,690 --> 00:06:07,460 is du dt equals f of u. 101 00:06:07,460 --> 00:06:12,030 The midpoint rule gives me du dt plus something times 102 00:06:12,030 --> 00:06:17,970 delta t squared and plus an even smaller terms equals f of u. 103 00:06:17,970 --> 00:06:24,640 So these additional term that is append from the discretization 104 00:06:24,640 --> 00:06:27,820 is the truncation error of the scheme, 105 00:06:27,820 --> 00:06:31,450 is the error made by the scheme, is the approximation. 106 00:06:31,450 --> 00:06:37,160 How much accuracy do I lose by approximating 107 00:06:37,160 --> 00:06:42,900 the derivative using discretization scheme? 108 00:06:42,900 --> 00:06:43,407 All right. 109 00:06:43,407 --> 00:06:44,615 This is the truncation error. 110 00:06:44,615 --> 00:06:47,615 And if you do forward Euler, this truncation error, 111 00:06:47,615 --> 00:06:51,300 would it be delta t squared? 112 00:06:51,300 --> 00:06:51,940 No. 113 00:06:51,940 --> 00:06:54,580 It will be on the order of delta t. 114 00:06:54,580 --> 00:06:59,564 So forward Euler is less accurate than midpoint. 115 00:06:59,564 --> 00:07:00,690 All right. 116 00:07:00,690 --> 00:07:02,440 What we are talking about is this 117 00:07:02,440 --> 00:07:04,970 is called local truncation error. 118 00:07:04,970 --> 00:07:08,454 The reason why we call it local truncation error is going 119 00:07:08,454 --> 00:07:11,220 to be clear in the next section because there 120 00:07:11,220 --> 00:07:16,300 is a different way of assessing the error is the global error, 121 00:07:16,300 --> 00:07:18,384 and we're going to be talking about that next one. 122 00:07:18,384 --> 00:07:18,884 Yes. 123 00:07:18,884 --> 00:07:19,878 AUDIENCE: [INAUDIBLE]. 124 00:07:32,800 --> 00:07:33,930 QIQI WANG: Good question. 125 00:07:33,930 --> 00:07:38,750 So the question here is, why did I start expanding here. 126 00:07:38,750 --> 00:07:41,980 What if I happened to expand all the way here? 127 00:07:41,980 --> 00:07:45,090 I can write a hundred terms after this. 128 00:07:45,090 --> 00:07:49,480 There is infinitely many terms in the [INAUDIBLE] 129 00:07:49,480 --> 00:07:54,240 I stopped over here because I did this, 130 00:07:54,240 --> 00:07:56,764 and I know that this term is not attainable. 131 00:07:59,890 --> 00:08:02,650 If a term is not unattainable then 132 00:08:02,650 --> 00:08:08,462 this is going to be the leading term in your local truncation 133 00:08:08,462 --> 00:08:09,910 error. 134 00:08:09,910 --> 00:08:14,095 The way you do your own is if I do an analysis of a scheme I've 135 00:08:14,095 --> 00:08:19,900 never seen before, I would expand it somewhere maybe 136 00:08:19,900 --> 00:08:27,150 onto here, and then dug it in to see if I get 137 00:08:27,150 --> 00:08:28,970 something that is not canceled. 138 00:08:28,970 --> 00:08:32,740 What happens if I see everything cancels? 139 00:08:32,740 --> 00:08:34,512 What if I see everything cancels, 140 00:08:34,512 --> 00:08:38,499 and all I have is exactly the [? C. ?] Why I might have 141 00:08:38,499 --> 00:08:40,980 an [? o ?] delta dt? 142 00:08:40,980 --> 00:08:42,639 AUDIENCE: [INAUDIBLE]. 143 00:08:42,639 --> 00:08:44,285 QIQI WANG: I should expand more. 144 00:08:44,285 --> 00:08:45,300 Exactly. 145 00:08:45,300 --> 00:08:50,490 That means that delta t2 may be the leading 146 00:08:50,490 --> 00:08:54,050 term of the truncation error, but would it 147 00:08:54,050 --> 00:08:56,435 be the leading term of truncation error 148 00:08:56,435 --> 00:08:59,430 if it is given less than delta t2? 149 00:08:59,430 --> 00:09:00,970 I don't know. 150 00:09:00,970 --> 00:09:04,595 I only know that the scheme is at least third order accurate 151 00:09:04,595 --> 00:09:06,820 [INAUDIBLE]. 152 00:09:06,820 --> 00:09:09,043 Only when I expand it more terms, 153 00:09:09,043 --> 00:09:14,190 I can know it happened in what order of accuracy this means. 154 00:09:14,190 --> 00:09:15,110 All right. 155 00:09:15,110 --> 00:09:20,190 So when you [? file ?] a scheme, when you're not sure, 156 00:09:20,190 --> 00:09:23,560 expand it to whatever you like and see if everything cancels. 157 00:09:23,560 --> 00:09:26,030 If everything does cancel, then it 158 00:09:26,030 --> 00:09:30,990 means you need to go back and expand more terms. 159 00:09:30,990 --> 00:09:31,770 Is it clear? 160 00:09:35,520 --> 00:09:40,180 In the last lecture, I asked you to do this to find out 161 00:09:40,180 --> 00:09:42,872 who the best X schemes. 162 00:09:42,872 --> 00:09:47,800 Can I have a show of hand of who did this? 163 00:09:47,800 --> 00:09:48,380 Who did this? 164 00:09:51,260 --> 00:09:53,730 Who attempted? 165 00:09:53,730 --> 00:09:56,400 Good. 166 00:09:56,400 --> 00:10:00,552 Who attempted the first one? 167 00:10:00,552 --> 00:10:03,706 Let me see, 1, 2, 3, 4, 5. 168 00:10:06,650 --> 00:10:09,310 I don't need you to be successful, just tell me 169 00:10:09,310 --> 00:10:10,960 what is the best you got. 170 00:10:14,690 --> 00:10:16,960 AUDIENCE: The first order [INAUDIBLE]. 171 00:10:16,960 --> 00:10:17,940 QIQI WANG: No. 172 00:10:17,940 --> 00:10:19,320 What is the scheme? 173 00:10:19,320 --> 00:10:22,190 So I want the product is, what is the-- 174 00:10:22,190 --> 00:10:23,670 AUDIENCE: [INAUDIBLE] coefficients? 175 00:10:23,670 --> 00:10:24,390 QIQI WANG: Yes. 176 00:10:24,390 --> 00:10:26,322 AUDIENCE: I just got r1. 177 00:10:26,322 --> 00:10:28,580 QIQI WANG: Did you get-- 178 00:10:28,580 --> 00:10:30,644 AUDIENCE: I got u of t sub delta t. 179 00:10:30,644 --> 00:10:31,310 QIQI WANG: Good. 180 00:10:31,310 --> 00:10:35,840 You get u of [INAUDIBLE]? 181 00:10:35,840 --> 00:10:37,210 AUDIENCE: Yeah. 182 00:10:37,210 --> 00:10:38,640 QIQI WANG: --t plus delta t. 183 00:10:38,640 --> 00:10:41,229 AUDIENCE: --is equal to u prime. 184 00:10:41,229 --> 00:10:42,770 QIQI WANG: --is equal to u prime at-- 185 00:10:42,770 --> 00:10:44,354 AUDIENCE: --t plus delta t. 186 00:10:44,354 --> 00:10:45,520 QIQI WANG: --t plus delta t. 187 00:10:45,520 --> 00:10:46,900 AUDIENCE: --plus u of t. 188 00:10:46,900 --> 00:10:48,170 QIQI WANG: --plus u of t. 189 00:10:48,170 --> 00:10:49,660 AUDIENCE: --plus u prime of t. 190 00:10:49,660 --> 00:10:52,152 QIQI WANG: --plus u prime of t. 191 00:10:55,160 --> 00:10:58,440 Anybody get anything different, or everybody agrees with that? 192 00:11:02,910 --> 00:11:04,940 Five people got the same thing? 193 00:11:04,940 --> 00:11:07,652 AUDIENCE: [INAUDIBLE]. 194 00:11:07,652 --> 00:11:09,610 QIQI WANG: OK. 195 00:11:09,610 --> 00:11:12,870 Let me write this in a more concise manner. 196 00:11:12,870 --> 00:11:16,650 t plus delta t is [? as ?] in f times x. 197 00:11:16,650 --> 00:11:18,440 t is at the current time set. 198 00:11:18,440 --> 00:11:21,160 So let me denote everything, t plus [INAUDIBLE] 199 00:11:21,160 --> 00:11:25,320 i plus 1, i plus 1, and this is i and i. 200 00:11:29,136 --> 00:11:34,460 So it's just clearer not having all the tn, t plus delta t. 201 00:11:34,460 --> 00:11:36,810 If I know-- come on. 202 00:11:36,810 --> 00:11:39,450 Come on. 203 00:11:39,450 --> 00:11:42,340 And you are sure there is no delta t in front of this? 204 00:11:46,110 --> 00:11:50,480 I mean, forward Euler, and at midpoint, I have delta t's. 205 00:11:50,480 --> 00:11:50,980 Yeah? 206 00:11:50,980 --> 00:11:52,828 AUDIENCE: [INAUDIBLE]. 207 00:11:52,828 --> 00:11:53,630 QIQI WANG: OK. 208 00:11:53,630 --> 00:11:57,040 All the derivative terms should have delta t. 209 00:11:57,040 --> 00:11:58,631 AUDIENCE: [INAUDIBLE]. 210 00:11:58,631 --> 00:12:01,327 QIQI WANG: You don't have that? 211 00:12:01,327 --> 00:12:02,243 AUDIENCE: [INAUDIBLE]. 212 00:12:12,940 --> 00:12:16,510 QIQI WANG: I was listing delta t's. 213 00:12:16,510 --> 00:12:20,510 Sorry for that, but if I do the theta series 214 00:12:20,510 --> 00:12:23,310 and I just [? carefully ?], correctly, I should test that. 215 00:12:23,310 --> 00:12:25,606 And if you got something different, 216 00:12:25,606 --> 00:12:27,800 let me write it here. 217 00:12:27,800 --> 00:12:33,210 You've got something ui plus 1 equal to 1/2 of ui 218 00:12:33,210 --> 00:12:41,500 plus 1 prime delta t plus ui plus, you also get 1/2 here. 219 00:12:41,500 --> 00:12:42,000 Yes? 220 00:12:42,000 --> 00:12:43,260 You got something different? 221 00:12:43,260 --> 00:12:44,200 AUDIENCE: [INAUDIBLE]. 222 00:12:47,020 --> 00:12:48,520 QIQI WANG: u prime of delta t. 223 00:12:48,520 --> 00:12:51,750 You get negative 1/2 here. 224 00:12:51,750 --> 00:12:52,250 All right. 225 00:12:52,250 --> 00:12:55,840 So we get three answers. 226 00:12:55,840 --> 00:12:58,680 Is that all? 227 00:12:58,680 --> 00:13:00,011 Anybody get anything different? 228 00:13:02,970 --> 00:13:03,480 Good. 229 00:13:03,480 --> 00:13:05,610 We have five people got three answers, 230 00:13:05,610 --> 00:13:07,562 so at least there are some agreements. 231 00:13:11,880 --> 00:13:16,120 Who did the second one, best explicit, two-step scheme? 232 00:13:16,120 --> 00:13:18,274 Who attempted the second one? 233 00:13:18,274 --> 00:13:19,696 You? 234 00:13:19,696 --> 00:13:25,591 So we have 1, 2, 3, 4, 5, 6, 7. 235 00:13:25,591 --> 00:13:27,800 More people did this. 236 00:13:27,800 --> 00:13:30,530 What is the best scheme you guys have got? 237 00:13:33,530 --> 00:13:34,530 AUDIENCE: [INAUDIBLE]. 238 00:13:38,600 --> 00:13:41,740 QIQI WANG: The question is, can she-- what's your name? 239 00:13:41,740 --> 00:13:43,195 AUDIENCE: [INAUDIBLE]. 240 00:13:43,195 --> 00:13:44,300 QIQI WANG: OK. 241 00:13:44,300 --> 00:13:49,880 Tren's question is, can we express the derivatives 242 00:13:49,880 --> 00:13:51,960 as f as well? 243 00:13:51,960 --> 00:13:53,330 The answer is yes. 244 00:13:53,330 --> 00:13:56,870 So whenever I see a derivative, I can plug in the derivative 245 00:13:56,870 --> 00:13:58,640 into the differential equation. 246 00:13:58,640 --> 00:14:01,900 The differential equation is u prime equal to fu. 247 00:14:01,900 --> 00:14:04,720 So whenever I have f prime of i plus 1, 248 00:14:04,720 --> 00:14:09,150 I can replace this as f of ui plus 1. 249 00:14:09,150 --> 00:14:15,860 Whenever I have f prime of i, I can replace this as f of ui. 250 00:14:15,860 --> 00:14:22,590 So what is your answer of the best explicit two-step scheme? 251 00:14:22,590 --> 00:14:23,544 AUDIENCE: [INAUDIBLE]. 252 00:14:31,653 --> 00:14:33,010 QIQI WANG: Negative 4ui? 253 00:14:33,010 --> 00:14:33,510 OK. 254 00:14:33,510 --> 00:14:34,482 AUDIENCE: [INAUDIBLE]. 255 00:14:49,770 --> 00:14:52,390 QIQI WANG: f ui and-- 256 00:14:52,390 --> 00:14:59,674 AUDIENCE: --plus 2 f times u [? sub ?] i minus 1. 257 00:14:59,674 --> 00:15:00,560 QIQI WANG: OK. 258 00:15:00,560 --> 00:15:03,190 Good. 259 00:15:03,190 --> 00:15:05,700 Anybody who get anything different for best explicit 260 00:15:05,700 --> 00:15:06,470 two-step scheme? 261 00:15:10,690 --> 00:15:12,021 You all got the same thing? 262 00:15:12,021 --> 00:15:12,903 AUDIENCE: [INAUDIBLE] 263 00:15:12,903 --> 00:15:15,350 QIQI WANG: Delta t's. 264 00:15:15,350 --> 00:15:18,490 And I guess delta t here. 265 00:15:18,490 --> 00:15:19,920 You must have that too. 266 00:15:19,920 --> 00:15:21,270 All right. 267 00:15:21,270 --> 00:15:22,175 Oh, wow. 268 00:15:22,175 --> 00:15:22,800 That's amazing. 269 00:15:22,800 --> 00:15:25,190 Everybody got this? 270 00:15:25,190 --> 00:15:28,130 No disagreements? 271 00:15:28,130 --> 00:15:33,220 So the explicit two-step people are more in agreement 272 00:15:33,220 --> 00:15:36,160 than implicit one-step people. 273 00:15:36,160 --> 00:15:38,950 Who did best implicit two-step scheme? 274 00:15:38,950 --> 00:15:40,260 Who attempted? 275 00:15:40,260 --> 00:15:41,010 You also did that? 276 00:15:41,010 --> 00:15:43,290 OK, great. 277 00:15:43,290 --> 00:15:44,659 And what is the answer? 278 00:15:48,012 --> 00:15:48,970 AUDIENCE: [INAUDIBLE]. 279 00:15:48,970 --> 00:15:51,144 QIQI WANG: OK, you. 280 00:15:51,144 --> 00:15:52,060 AUDIENCE: [INAUDIBLE]. 281 00:16:02,680 --> 00:16:05,950 QIQI WANG: I'm going to append a delta to for you. 282 00:16:05,950 --> 00:16:06,930 Plus-- 283 00:16:06,930 --> 00:16:07,910 AUDIENCE: [INAUDIBLE]. 284 00:16:20,650 --> 00:16:22,120 QIQI WANG: Wonderful. 285 00:16:22,120 --> 00:16:25,264 So can you tell me how you did that? 286 00:16:25,264 --> 00:16:31,556 AUDIENCE: I did a bunch of Taylor [INAUDIBLE]. 287 00:16:31,556 --> 00:16:32,315 QIQI WANG: Nice. 288 00:16:34,980 --> 00:16:39,920 I get to have you show me how you did that later on. 289 00:16:39,920 --> 00:16:42,662 Who the best explicit three-step scheme? 290 00:16:42,662 --> 00:16:43,370 Anybody try that? 291 00:16:46,280 --> 00:16:47,360 You also tried that? 292 00:16:47,360 --> 00:16:48,590 OK, great. 293 00:16:48,590 --> 00:16:52,340 I'm going to have you as an instructor here. 294 00:16:52,340 --> 00:16:53,902 What is the answer for that? 295 00:16:53,902 --> 00:16:54,818 AUDIENCE: [INAUDIBLE]. 296 00:17:00,424 --> 00:17:01,340 QIQI WANG: Minus what? 297 00:17:01,340 --> 00:17:01,600 Sorry? 298 00:17:01,600 --> 00:17:02,556 AUDIENCE: [INAUDIBLE]. 299 00:17:08,292 --> 00:17:11,020 QIQI WANG: 18 and sorry, what? 300 00:17:11,020 --> 00:17:13,948 AUDIENCE: 9. 301 00:17:13,948 --> 00:17:14,940 QIQI WANG: 9. 302 00:17:14,940 --> 00:17:20,231 AUDIENCE: 9, 18, 10, and 3. 303 00:17:20,231 --> 00:17:24,880 QIQI WANG: 18, 10, and 3. 304 00:17:24,880 --> 00:17:26,967 So what is multiplied on this one? 305 00:17:26,967 --> 00:17:27,961 AUDIENCE: [INAUDIBLE] 306 00:17:44,859 --> 00:17:48,950 QIQI WANG: ui minus 2 and ui minus 2 prime. 307 00:17:48,950 --> 00:17:52,688 So all the primes have delta t, I guess. 308 00:17:56,600 --> 00:17:59,070 Great. 309 00:17:59,070 --> 00:18:04,230 So I'm going to do some analysis of what this is. 310 00:18:04,230 --> 00:18:06,140 I'm going to go back to this. 311 00:18:06,140 --> 00:18:08,722 So I saved the documents. 312 00:18:08,722 --> 00:18:12,470 Even if I lose it, I'm going to go back to this. 313 00:18:12,470 --> 00:18:22,110 But what I want to show right now is first of all, 314 00:18:22,110 --> 00:18:24,760 I'm sure some of you did this, so you know already 315 00:18:24,760 --> 00:18:26,300 how to do this. 316 00:18:26,300 --> 00:18:29,260 But I also want to see especially 317 00:18:29,260 --> 00:18:31,920 resolve the disagreement on the first one. 318 00:18:31,920 --> 00:18:34,570 So that is what I'm going to do first. 319 00:18:34,570 --> 00:18:37,850 What is the best implicit, one-step scheme? 320 00:18:37,850 --> 00:18:41,540 I'm going to be having a demonstration of that right 321 00:18:41,540 --> 00:18:45,330 now, and then I'm going to show you 322 00:18:45,330 --> 00:18:48,650 how to implement a scheme once you derive it. 323 00:18:48,650 --> 00:18:50,530 You derived a lot of schemes, now how 324 00:18:50,530 --> 00:18:55,600 do I put this scheme into a code, into a production, 325 00:18:55,600 --> 00:18:58,610 into something actually wrong. 326 00:18:58,610 --> 00:19:02,700 So first of all the best implicit one-step scheme. 327 00:19:02,700 --> 00:19:04,280 All right. 328 00:19:04,280 --> 00:19:11,835 I'm going to append something after this. 329 00:19:16,030 --> 00:19:21,630 The best one-step scheme. 330 00:19:21,630 --> 00:19:29,756 So the scheme has to be ui plus 1 equal to something times ui. 331 00:19:29,756 --> 00:19:30,590 All right. 332 00:19:33,180 --> 00:19:37,330 Did I say implicit, one-step scheme. 333 00:19:37,330 --> 00:19:39,800 So one-step scheme means everything 334 00:19:39,800 --> 00:19:44,220 has to be written in terms of i plus 1 and i. 335 00:19:44,220 --> 00:19:51,290 And it can depend on ui, ui prime, and ui plus 1 prime. 336 00:19:51,290 --> 00:19:53,780 This is because it can be implicit. 337 00:19:53,780 --> 00:19:57,090 If it's explicit, then it has to always depend on i. 338 00:19:57,090 --> 00:19:59,530 All right. 339 00:19:59,530 --> 00:20:08,020 So let me put in three unknowns, x plus y plus z. 340 00:20:08,020 --> 00:20:12,490 And these are unknowns. 341 00:20:12,490 --> 00:20:17,706 And what do I do to make the scheme as accurate as possible? 342 00:20:17,706 --> 00:20:19,410 AUDIENCE: Minimize the error? 343 00:20:19,410 --> 00:20:21,030 QIQI WANG: Minimize the error. 344 00:20:21,030 --> 00:20:21,820 Exactly. 345 00:20:21,820 --> 00:20:24,330 And the error can only be [? appends ?] 346 00:20:24,330 --> 00:20:27,340 from Taylor series expansion All right. 347 00:20:27,340 --> 00:20:31,980 So Taylor series expansion, again, there are two options. 348 00:20:31,980 --> 00:20:37,620 One option is to expand all the i plus 1's at i. 349 00:20:37,620 --> 00:20:43,310 So for example, ui plus 1 is equal to ui plus ui prime delta 350 00:20:43,310 --> 00:20:49,520 t plus 1/2 of ui double prime delta t square plus-- I 351 00:20:49,520 --> 00:20:52,900 can do more terms-- triple prime delta 352 00:20:52,900 --> 00:20:57,370 t cubed plus o delta t fourth. 353 00:20:57,370 --> 00:21:02,290 But actually, in this case, another equally convenient way 354 00:21:02,290 --> 00:21:06,660 is to expand ui's in terms of ui's plus 1's 355 00:21:06,660 --> 00:21:12,110 because there are as many terms that is i plus 1 356 00:21:12,110 --> 00:21:14,136 as there are many terms of ui. 357 00:21:14,136 --> 00:21:19,730 So it is equally convenient, so equally cumbersome, 358 00:21:19,730 --> 00:21:22,590 whichever way you think about it. 359 00:21:22,590 --> 00:21:29,070 u prime of i plus 1, I can also do a Taylor series expansion 360 00:21:29,070 --> 00:21:31,780 also by just taking the derivative. 361 00:21:31,780 --> 00:21:34,650 Now, I'm going to write everything [? correspondence ?] 362 00:21:34,650 --> 00:21:40,050 here, but every term has one more derivative 363 00:21:40,050 --> 00:21:41,480 than what is above. 364 00:21:45,150 --> 00:21:47,690 So when I have third derivative here, 365 00:21:47,690 --> 00:21:51,170 I actually have fourth derivative. 366 00:21:51,170 --> 00:21:53,990 And this is still to the fourth. 367 00:21:53,990 --> 00:21:55,975 So all the delta t's to the power is the same, 368 00:21:55,975 --> 00:21:59,140 but all the derivatives have one more order 369 00:21:59,140 --> 00:22:01,960 because I'm expanding the derivatives. 370 00:22:01,960 --> 00:22:04,290 And when I should take the second load of derivatives, 371 00:22:04,290 --> 00:22:06,980 I should be taking the second derivative of the derivative, 372 00:22:06,980 --> 00:22:08,700 which is the third derivative. 373 00:22:08,700 --> 00:22:09,460 All right. 374 00:22:12,262 --> 00:22:15,685 Is it clear why I'm having one more derivative here? 375 00:22:19,360 --> 00:22:22,360 Now, this is all I need to expand. 376 00:22:22,360 --> 00:22:27,750 And to figure out the error, we need to plug these extensions 377 00:22:27,750 --> 00:22:29,480 into the formula. 378 00:22:29,480 --> 00:22:35,320 We also need to plot something else in the formula. 379 00:22:35,320 --> 00:22:38,620 Oh, in this case, we don't need to plug anything else. 380 00:22:38,620 --> 00:22:45,560 If I were writing the scheme not as u prime, but as X, 381 00:22:45,560 --> 00:22:47,830 then I also need to plug in the differential equation 382 00:22:47,830 --> 00:22:48,630 into the formula. 383 00:22:48,630 --> 00:22:53,586 But I'm already writing the scheme in terms of u prime, 384 00:22:53,586 --> 00:22:55,960 then I don't need to plug anything else into the formula. 385 00:22:55,960 --> 00:22:59,910 I just need to plug in the Taylor expansions. 386 00:22:59,910 --> 00:23:05,640 When I plug in the Taylor expansions, what I get is, 387 00:23:05,640 --> 00:23:08,040 I'm going to write down the truncation error, which 388 00:23:08,040 --> 00:23:12,140 is the left-hand side. 389 00:23:12,140 --> 00:23:17,030 Well, this case, it is actually the truncation era times delta 390 00:23:17,030 --> 00:23:18,548 t is equal to this. 391 00:23:18,548 --> 00:23:21,900 I'm going to let you know why I'm writing this. 392 00:23:21,900 --> 00:23:23,780 But the truncation error times delta t 393 00:23:23,780 --> 00:23:26,742 is the left-hand side minus the right-hand side. 394 00:23:31,340 --> 00:23:35,510 And let me just write down term by term. 395 00:23:35,510 --> 00:23:40,780 The first term is ui plus, which is this. 396 00:23:40,780 --> 00:23:45,840 Let me actually start by changing the color. 397 00:23:45,840 --> 00:23:52,145 It is ui plus ui prime times delta t 398 00:23:52,145 --> 00:23:58,290 plus ui double prime times delta t squared over 2. 399 00:23:58,290 --> 00:24:11,070 plus ui triple prime delta t cubed over 6 plus our delta t 400 00:24:11,070 --> 00:24:11,700 to the fourth. 401 00:24:14,940 --> 00:24:19,110 The second term is minus x times ui-- 402 00:24:19,110 --> 00:24:24,140 let me changes the color to this-- minus x times ui. 403 00:24:26,680 --> 00:24:27,950 Now, what is the third term? 404 00:24:27,950 --> 00:24:31,120 It's minus y times ui prime. 405 00:24:31,120 --> 00:24:35,610 So it is minus y times ui prime. 406 00:24:35,610 --> 00:24:38,525 I'm writing down here so that everything responding to prime 407 00:24:38,525 --> 00:24:42,000 is on this column. 408 00:24:42,000 --> 00:24:46,490 And the third term is minus z times ui plus 1 prime, 409 00:24:46,490 --> 00:24:49,190 and we again plug in the Taylor expansion. 410 00:24:51,920 --> 00:24:54,160 But everything starts at the prime, 411 00:24:54,160 --> 00:24:59,495 so I have a minus z times ui prime here. 412 00:24:59,495 --> 00:25:07,790 I have minus ui dot double prime times z delta t. 413 00:25:07,790 --> 00:25:09,420 So this is double prime. 414 00:25:09,420 --> 00:25:13,430 This is z times this, and I also have this. 415 00:25:13,430 --> 00:25:21,650 I have minus ui triple prime times 1/2 of z delta t 416 00:25:21,650 --> 00:25:23,170 to the cube. 417 00:25:23,170 --> 00:25:26,040 And I don't even need to write this. 418 00:25:26,040 --> 00:25:37,190 I'm just going to write to this as o delta t cubed times z. 419 00:25:37,190 --> 00:25:39,660 All right. 420 00:25:39,660 --> 00:25:40,656 AUDIENCE: [INAUDIBLE]. 421 00:25:44,640 --> 00:25:46,094 QIQI WANG: Up more? 422 00:25:46,094 --> 00:25:47,260 AUDIENCE: No, the other way. 423 00:25:47,260 --> 00:25:48,680 QIQI WANG: The other way. 424 00:25:48,680 --> 00:25:50,160 Scroll up. 425 00:25:50,160 --> 00:25:51,750 Right. 426 00:25:51,750 --> 00:25:54,450 So what I'm doing is I'm expanding this term 427 00:25:54,450 --> 00:25:57,750 as the first line, this term at the second line, 428 00:25:57,750 --> 00:26:01,750 this term as the third line, and this is the fourth line. 429 00:26:01,750 --> 00:26:06,440 I'm just writing things down, and it's different lines, 430 00:26:06,440 --> 00:26:11,330 and I'm also aligning the same derivatives of u 431 00:26:11,330 --> 00:26:13,130 on the same color. 432 00:26:13,130 --> 00:26:14,154 Is it clear? 433 00:26:17,870 --> 00:26:21,390 So how do I choose the coefficients? 434 00:26:25,390 --> 00:26:28,230 I'm going to choose coefficients to cancel 435 00:26:28,230 --> 00:26:32,380 as many terms as possible. 436 00:26:32,380 --> 00:26:36,110 I'm going to choose, first of all, 437 00:26:36,110 --> 00:26:41,850 the first priority is to cancel anything that is 438 00:26:41,850 --> 00:26:44,470 ui, anything multiplied on ui. 439 00:26:44,470 --> 00:26:46,742 So what x should I choose? 440 00:26:46,742 --> 00:26:48,170 AUDIENCE: [INAUDIBLE]. 441 00:26:48,170 --> 00:26:49,330 QIQI WANG: Exactly. 442 00:26:49,330 --> 00:26:54,540 x has to be equal to 1 because I want these two terms, 443 00:26:54,540 --> 00:26:57,590 I want them to be canceled. 444 00:26:57,590 --> 00:27:00,740 The only way to cancel this is we choose x equal to 1. 445 00:27:00,740 --> 00:27:02,020 Is it clear? 446 00:27:04,610 --> 00:27:13,660 And what should I do to ensure that this term also cancels? 447 00:27:13,660 --> 00:27:14,650 AUDIENCE: [INAUDIBLE]. 448 00:27:23,460 --> 00:27:24,370 QIQI WANG: What? 449 00:27:24,370 --> 00:27:24,870 Somebody? 450 00:27:24,870 --> 00:27:26,742 AUDIENCE: [INAUDIBLE]. 451 00:27:26,742 --> 00:27:28,974 QIQI WANG: y plus z equal to 1. 452 00:27:28,974 --> 00:27:29,890 AUDIENCE: [INAUDIBLE]. 453 00:27:29,890 --> 00:27:30,880 QIQI WANG: Equal to delta t. 454 00:27:30,880 --> 00:27:31,410 Exactly. 455 00:27:31,410 --> 00:27:34,970 y plus z has to equal to delta t in order 456 00:27:34,970 --> 00:27:38,440 for this, this, and this to cancel. 457 00:27:38,440 --> 00:27:42,096 y plus z has equal to delta t. 458 00:27:42,096 --> 00:27:44,030 Do we have a question? 459 00:27:44,030 --> 00:27:46,170 No? 460 00:27:46,170 --> 00:27:48,740 All right. 461 00:27:48,740 --> 00:27:51,055 What can we do to cancel this? 462 00:27:57,826 --> 00:28:01,270 AUDIENCE: [INAUDIBLE]. 463 00:28:01,270 --> 00:28:05,160 QIQI WANG: We need to have minus z delta t 464 00:28:05,160 --> 00:28:09,460 to cancel with the last delta t squared over 2. 465 00:28:09,460 --> 00:28:15,480 So we have z delta t has to equal 466 00:28:15,480 --> 00:28:19,720 to delta t squared over 2, which means z 467 00:28:19,720 --> 00:28:21,760 has to equal to 1/2 of delta t. 468 00:28:26,960 --> 00:28:30,010 Now, the question is, can I cancel more? 469 00:28:30,010 --> 00:28:33,540 Can I make sure I can cancel more? 470 00:28:33,540 --> 00:28:34,040 No. 471 00:28:34,040 --> 00:28:38,372 I only have three terms to play around with. 472 00:28:38,372 --> 00:28:43,320 I only have three unknowns, x, y, and z. 473 00:28:43,320 --> 00:28:46,550 So I can only satisfy three equations. 474 00:28:46,550 --> 00:28:50,060 On [INAUDIBLE] that of the three equations. 475 00:28:50,060 --> 00:28:53,660 There is [INAUDIBLE] that satisfy automatically. 476 00:28:53,660 --> 00:28:56,050 So unless I'm super lucky, I can only 477 00:28:56,050 --> 00:29:02,980 tune the three unknowns to satisfy three linear equations. 478 00:29:02,980 --> 00:29:06,130 All right. 479 00:29:06,130 --> 00:29:08,610 So in this case, it is pretty clear 480 00:29:08,610 --> 00:29:13,850 that z has to equal to 1/2 of t and y 481 00:29:13,850 --> 00:29:15,610 because of the second equation also 482 00:29:15,610 --> 00:29:19,350 has to be equal to 1/2 of t. 483 00:29:19,350 --> 00:29:22,780 Or I can just take these three equations and make it a matrix. 484 00:29:25,920 --> 00:29:30,920 I can take these three equations and particularly, I 485 00:29:30,920 --> 00:29:39,050 can make, if I can-- well, let me see. 486 00:29:39,050 --> 00:29:46,290 If I can divide the second and third equations by delta t 487 00:29:46,290 --> 00:29:50,810 so that the unknowns are x and y over delta t and z 488 00:29:50,810 --> 00:29:51,990 over delta t. 489 00:29:51,990 --> 00:29:55,810 If I make them as the unknowns, then I 490 00:29:55,810 --> 00:30:09,470 can write the equation as x, y, z equal to 1, 1, and 1/2. 491 00:30:09,470 --> 00:30:10,600 How can I do that? 492 00:30:10,600 --> 00:30:12,140 How can I fill in this matrix? 493 00:30:12,140 --> 00:30:13,056 AUDIENCE: [INAUDIBLE]. 494 00:30:13,056 --> 00:30:14,990 QIQI WANG: Well, the first line is 1, 0, 0. 495 00:30:14,990 --> 00:30:16,840 Thank you. 496 00:30:16,840 --> 00:30:19,146 It's just the first equation x equal to 1. 497 00:30:19,146 --> 00:30:20,270 What about the second line? 498 00:30:20,270 --> 00:30:22,970 I mean, this is not y, z, but y over 499 00:30:22,970 --> 00:30:25,028 delta t and z over delta t. 500 00:30:25,028 --> 00:30:27,012 What is the second line? 501 00:30:27,012 --> 00:30:28,996 AUDIENCE: [INAUDIBLE]. 502 00:30:28,996 --> 00:30:30,770 QIQI WANG: 0, 1, 1. 503 00:30:30,770 --> 00:30:31,830 Exactly. 504 00:30:31,830 --> 00:30:35,260 And what is the third line? 505 00:30:35,260 --> 00:30:36,750 0, 0, 1. 506 00:30:36,750 --> 00:30:40,470 I mean, this looks obvious, but as you 507 00:30:40,470 --> 00:30:46,270 get to more complex schemes, as what's your name again? 508 00:30:46,270 --> 00:30:47,196 AUDIENCE: Jacobi. 509 00:30:47,196 --> 00:30:48,010 QIQI WANG: Jacobi. 510 00:30:48,010 --> 00:30:52,490 Jacobi has discovered in his [? one ?] contact schemes, 511 00:30:52,490 --> 00:30:55,600 you cannot just by looking at equations and get all 512 00:30:55,600 --> 00:30:56,310 the coefficients. 513 00:30:56,310 --> 00:30:58,920 You have to solve a matrix equation like that. 514 00:30:58,920 --> 00:31:01,880 And in MATLAB, let me set up this matrix. 515 00:31:01,880 --> 00:31:08,390 This matrix is 1, 0, 0, the first line. 516 00:31:08,390 --> 00:31:10,657 The second line is 0, 1, 1. 517 00:31:10,657 --> 00:31:14,020 The third line 0, 0, 1. 518 00:31:14,020 --> 00:31:15,160 This is the matrix. 519 00:31:15,160 --> 00:31:15,660 Exactly. 520 00:31:15,660 --> 00:31:17,320 The matrix. 521 00:31:17,320 --> 00:31:22,484 And the right answer is 1, 1, and 1/2. 522 00:31:22,484 --> 00:31:24,448 Yeah, that's the right-hand side. 523 00:31:24,448 --> 00:31:28,780 And by a backslash c, does everybody 524 00:31:28,780 --> 00:31:30,440 know what the backslash means? 525 00:31:33,510 --> 00:31:34,010 Yes. 526 00:31:34,010 --> 00:31:34,885 It's matrix division. 527 00:31:34,885 --> 00:31:37,600 It is solving this equation. 528 00:31:37,600 --> 00:31:40,000 It is solving ax equal to b. 529 00:31:40,000 --> 00:31:43,460 It solves for x, and I push. 530 00:31:43,460 --> 00:31:46,850 It give me the coefficients, 1, 1/2, 1/2, 531 00:31:46,850 --> 00:31:51,180 which means x equal 1, y over delta 532 00:31:51,180 --> 00:31:54,190 t equal to 1/2, which means y equal to 1/2 delta t, 533 00:31:54,190 --> 00:31:57,302 and z also equal to 1/2 delta t. 534 00:31:57,302 --> 00:31:58,410 All right. 535 00:31:58,410 --> 00:32:05,040 So some of the scheme we have is ui plus 1 536 00:32:05,040 --> 00:32:12,450 equal to ui plus 1/2 delta t ui prime plus 1/2 delta t ui 537 00:32:12,450 --> 00:32:13,450 plus 1 prime. 538 00:32:16,240 --> 00:32:17,635 All right. 539 00:32:17,635 --> 00:32:18,330 So good. 540 00:32:18,330 --> 00:32:20,040 I think we resolved the conflict. 541 00:32:22,920 --> 00:32:28,940 There is only one scheme that satisfy this matrix equation. 542 00:32:28,940 --> 00:32:31,250 And there is only one scheme that 543 00:32:31,250 --> 00:32:32,950 can make sure the first term cancels, 544 00:32:32,950 --> 00:32:36,406 the second term cancels, the third term cancels. 545 00:32:36,406 --> 00:32:37,690 All right. 546 00:32:37,690 --> 00:32:40,460 And remember z is 1/2 delta t. 547 00:32:40,460 --> 00:32:45,192 So what is the truncation error here? 548 00:32:45,192 --> 00:32:46,400 What is the truncation error? 549 00:32:46,400 --> 00:32:50,150 What is the other of the truncation error of the scheme? 550 00:32:50,150 --> 00:32:53,200 What local [? order of accuracy ?] 551 00:32:53,200 --> 00:32:54,310 does the scheme have? 552 00:32:54,310 --> 00:32:56,190 AUDIENCE: [INAUDIBLE]. 553 00:32:56,190 --> 00:32:58,360 QIQI WANG: Second order delta t square. 554 00:32:58,360 --> 00:33:06,210 Because this is square, so sorry. 555 00:33:06,210 --> 00:33:09,710 This is square, but z is equal to 1/2 of delta t, 556 00:33:09,710 --> 00:33:11,110 so this term is delta t too. 557 00:33:11,110 --> 00:33:14,720 So the truncation of the leading term of the delta t times tau 558 00:33:14,720 --> 00:33:19,240 is delta t too, which means tau is delta t square. 559 00:33:19,240 --> 00:33:21,590 All right. 560 00:33:21,590 --> 00:33:25,520 And why did I say this is delta t tau? 561 00:33:25,520 --> 00:33:34,280 This is because if you look at how accurate I am approximating 562 00:33:34,280 --> 00:33:36,995 the derivative, I mean, I'm approximating 563 00:33:36,995 --> 00:33:42,220 the average of the derivative between i and i plus 1 set. 564 00:33:42,220 --> 00:33:46,500 I'm basically approximating the average of the i-th derivative 565 00:33:46,500 --> 00:33:47,840 and i-th plus 1 derivative. 566 00:33:47,840 --> 00:33:51,560 It is equal to what? 567 00:33:51,560 --> 00:33:54,860 ui plus 1 minus ui divided by delta t. 568 00:33:57,540 --> 00:34:00,180 And if you look at the truncation 569 00:34:00,180 --> 00:34:04,780 error of this scheme, so if I do this rigorously, 570 00:34:04,780 --> 00:34:07,460 this is delta t cube here. 571 00:34:07,460 --> 00:34:12,429 And if you transform, go from this equation to this equation, 572 00:34:12,429 --> 00:34:14,634 you have to divide everything out by delta t. 573 00:34:14,634 --> 00:34:16,929 So you get delta t squared. 574 00:34:16,929 --> 00:34:21,965 So the truncation error of this scheme is delta t squared. 575 00:34:21,965 --> 00:34:22,550 All right. 576 00:34:22,550 --> 00:34:27,610 I'm basically taking these two terms, dividing by delta t, 577 00:34:27,610 --> 00:34:28,580 and getting this. 578 00:34:28,580 --> 00:34:32,744 Well, I'm moving this term to here and divide by delta t 579 00:34:32,744 --> 00:34:35,110 to get this. 580 00:34:35,110 --> 00:34:36,884 And what is that here? 581 00:34:36,884 --> 00:34:39,230 Because I'm dividing by delta t, I get delta t squared. 582 00:34:39,230 --> 00:34:42,783 So the scheme is second order accurate. 583 00:34:51,820 --> 00:34:54,850 Any questions? 584 00:34:54,850 --> 00:34:55,760 No? 585 00:34:55,760 --> 00:34:56,260 Here? 586 00:34:59,530 --> 00:35:07,540 Let's try to implement the scheme into one 587 00:35:07,540 --> 00:35:12,070 of the very simple problems. 588 00:35:12,070 --> 00:35:15,490 Let's do a pendulum problem. 589 00:35:15,490 --> 00:35:19,410 In a pendulum problem, we have a pendulum that 590 00:35:19,410 --> 00:35:22,480 is vibrating on a small angle. 591 00:35:22,480 --> 00:35:27,190 So today we restrict ourselves to small angles 592 00:35:27,190 --> 00:35:29,199 so that the equation is not going to be linear. 593 00:35:29,199 --> 00:35:31,490 And we're going to be going to [? nulling ?] equations, 594 00:35:31,490 --> 00:35:35,450 and we're going to see why [? nulling ?] equations is 595 00:35:35,450 --> 00:35:39,170 going to be more complex especially 596 00:35:39,170 --> 00:35:40,570 for implicit schemes. 597 00:35:40,570 --> 00:35:44,460 So if it's explicit scheme, actually, I'm 598 00:35:44,460 --> 00:35:46,480 going to show you later that it's 599 00:35:46,480 --> 00:35:49,240 very easy to do a non-linear scheme, 600 00:35:49,240 --> 00:35:53,230 but for implicit method, like the one 601 00:35:53,230 --> 00:35:58,190 we just derived, one that actually involves 602 00:35:58,190 --> 00:36:02,020 a u prime of i plus 1, it is going 603 00:36:02,020 --> 00:36:05,560 to be much more difficult to do a non-linear scheme. 604 00:36:05,560 --> 00:36:09,550 But here the variable is theta. 605 00:36:09,550 --> 00:36:15,590 And what determines the acceleration of theta? 606 00:36:15,590 --> 00:36:17,526 What is the second derivative of theta? 607 00:36:17,526 --> 00:36:19,504 What's that going to be determining? 608 00:36:22,456 --> 00:36:23,440 AUDIENCE: [INAUDIBLE]. 609 00:36:33,772 --> 00:36:35,380 QIQI WANG: Force balance, right. 610 00:36:35,380 --> 00:36:38,500 What is the acceleration of this ball? 611 00:36:38,500 --> 00:36:41,067 Let me just ask, what is the gradual acceleration 612 00:36:41,067 --> 00:36:42,061 of this ball? 613 00:36:42,061 --> 00:36:43,055 AUDIENCE: [INAUDIBLE]. 614 00:36:47,528 --> 00:36:48,522 QIQI WANG: [INAUDIBLE]. 615 00:36:51,520 --> 00:36:54,350 This is [INAUDIBLE], so the only way to accelerate 616 00:36:54,350 --> 00:36:59,000 is along the [? duct. ?] So it is [INAUDIBLE]. 617 00:37:02,443 --> 00:37:03,109 is acceleration? 618 00:37:10,594 --> 00:37:11,592 AUDIENCE: [INAUDIBLE]. 619 00:37:15,085 --> 00:37:16,580 QIQI WANG: It can be moving. 620 00:37:16,580 --> 00:37:18,513 Let's assume there is no friction. 621 00:37:18,513 --> 00:37:28,300 The only force it has is the gravity force and the tension 622 00:37:28,300 --> 00:37:32,362 of its spring, which is [? non-elastic. ?] 623 00:37:41,710 --> 00:37:43,678 AUDIENCE: g sine theta. 624 00:37:43,678 --> 00:37:45,450 QIQI WANG: g sine theta. 625 00:37:45,450 --> 00:37:48,540 Why is that g sine theta? 626 00:37:48,540 --> 00:37:49,500 AUDIENCE: [INAUDIBLE]. 627 00:37:57,180 --> 00:37:58,940 QIQI WANG: Yeah. 628 00:37:58,940 --> 00:38:01,410 a is equal to F over m. 629 00:38:04,300 --> 00:38:09,320 It's force over mass, and what is the force? 630 00:38:09,320 --> 00:38:11,010 The force is actually here. 631 00:38:11,010 --> 00:38:15,190 This is the force in interaction. 632 00:38:15,190 --> 00:38:22,740 So it is mg sine theta over m, which is exactly g sine theta. 633 00:38:22,740 --> 00:38:25,590 So you are right. 634 00:38:25,590 --> 00:38:26,089 Sorry. 635 00:38:26,089 --> 00:38:27,478 What's your name again? 636 00:38:27,478 --> 00:38:28,857 AUDIENCE: Libby. 637 00:38:28,857 --> 00:38:29,565 QIQI WANG: Libby. 638 00:38:29,565 --> 00:38:30,065 OK. 639 00:38:30,065 --> 00:38:32,060 Libby is right. 640 00:38:32,060 --> 00:38:33,640 So this is acceleration. 641 00:38:33,640 --> 00:38:38,330 And what is the acceleration in theta? 642 00:38:38,330 --> 00:38:42,280 What is the second derivative of theta? 643 00:38:42,280 --> 00:38:45,190 Just a geometric relation between this and deceleration. 644 00:38:53,126 --> 00:38:56,078 AUDIENCE: [INAUDIBLE]. 645 00:38:56,078 --> 00:38:58,640 QIQI WANG: Divide by L. Exactly. 646 00:38:58,640 --> 00:39:02,730 So this is a over L. So it is going 647 00:39:02,730 --> 00:39:08,100 to be g sine theta over L. All right. 648 00:39:10,630 --> 00:39:14,500 And we are going to be linearizing this equation 649 00:39:14,500 --> 00:39:18,530 because sine theta, when theta is very small, is what? 650 00:39:18,530 --> 00:39:21,860 Again, Taylor series expansion sine theta 651 00:39:21,860 --> 00:39:27,620 is equal to sine theta at 0 plus theta 652 00:39:27,620 --> 00:39:30,030 times the derivative of sine theta 653 00:39:30,030 --> 00:39:33,880 at 0, which is equal to 1, which is 654 00:39:33,880 --> 00:39:36,970 cosine theta at 0, et cetera. 655 00:39:36,970 --> 00:39:38,360 So this term is 0. 656 00:39:38,360 --> 00:39:43,739 This term is theta, so it is approximately equal to theta. 657 00:39:47,050 --> 00:39:51,190 I mean, this is probably sine theta 658 00:39:51,190 --> 00:39:53,573 can be approximated by theta a lot of times, 659 00:39:53,573 --> 00:39:57,660 but in order to know this is, again, a result of Taylor's 660 00:39:57,660 --> 00:39:58,330 rule expansion. 661 00:40:01,800 --> 00:40:04,610 So sine theta can be an approximate of theta, 662 00:40:04,610 --> 00:40:08,060 and now I'm going to write the question 663 00:40:08,060 --> 00:40:10,700 theta double dot is equal to g theta over L 664 00:40:10,700 --> 00:40:12,720 in a pretty strange way. 665 00:40:12,720 --> 00:40:25,926 I'm going to write it as theta doubled dot equal to-- no, 666 00:40:25,926 --> 00:40:27,426 I'm not going to write it like this. 667 00:40:30,620 --> 00:40:33,750 I'm going to be writing this as d theta 668 00:40:33,750 --> 00:40:40,500 dt equal to theta dot, which is just a definition of theta dot. 669 00:40:40,500 --> 00:40:43,310 It's equivalent. 670 00:40:43,310 --> 00:40:48,740 d theta dot dt is actually equal to, again, 671 00:40:48,740 --> 00:40:55,810 the second derivative, which is equal to g theta 672 00:40:55,810 --> 00:41:02,600 provided by L. You may be asking, why am I doing this? 673 00:41:02,600 --> 00:41:03,570 Again, I'm sorry. 674 00:41:03,570 --> 00:41:05,635 So there is a negative sign on this. 675 00:41:05,635 --> 00:41:09,991 I have forgotten because-- so if theta is positive, 676 00:41:09,991 --> 00:41:12,580 the acceleration is in the negative direction, 677 00:41:12,580 --> 00:41:13,900 so I'm using a negative sign. 678 00:41:16,880 --> 00:41:20,655 Yeah, if you discover something like this, please shout. 679 00:41:24,250 --> 00:41:28,320 I'm writing it this way because remember 680 00:41:28,320 --> 00:41:31,100 the schemes we have been deriving 681 00:41:31,100 --> 00:41:34,700 are also first-order ODEs. 682 00:41:34,700 --> 00:41:38,550 What are first-order ODE mean? 683 00:41:38,550 --> 00:41:40,840 It means I can only have first-order derivatives 684 00:41:40,840 --> 00:41:42,906 through time. 685 00:41:42,906 --> 00:41:44,730 All right. 686 00:41:44,730 --> 00:41:46,960 So if I have a second-order ODE here, 687 00:41:46,960 --> 00:41:50,970 I need to convert it into a first-order ODE. 688 00:41:50,970 --> 00:41:54,480 Actually, I cannot convert it into first-order ODE. 689 00:41:54,480 --> 00:42:00,440 I have to convert it into two coupled first-order ODEs. 690 00:42:00,440 --> 00:42:02,110 Whenever you have a second-order ODE, 691 00:42:02,110 --> 00:42:05,530 you can always convert it into two first-order ODEs. 692 00:42:05,530 --> 00:42:07,260 Whenever you have a third-order ODE, 693 00:42:07,260 --> 00:42:10,460 you can always convert it into three first-order ODEs 694 00:42:10,460 --> 00:42:16,690 by introducing basically a separate variable, theta dot. 695 00:42:16,690 --> 00:42:21,630 So here I am having theta and theta dot as two variables. 696 00:42:21,630 --> 00:42:24,950 The derivative of both theta and theta dot 697 00:42:24,950 --> 00:42:26,974 is a function of theta and theta dot. 698 00:42:31,510 --> 00:42:33,310 Is it clear? 699 00:42:33,310 --> 00:42:35,130 Let me do this. 700 00:42:35,130 --> 00:42:42,160 Let me go to MATLAB and start to write this. 701 00:42:42,160 --> 00:42:43,080 Let me do this. 702 00:42:43,080 --> 00:42:45,840 I'm going to create an empty file here. 703 00:42:45,840 --> 00:42:50,840 I'm going to find out why MATLAB doesn't listen to me. 704 00:42:50,840 --> 00:42:51,846 It's called pendulum.m. 705 00:42:55,010 --> 00:42:56,838 Open Pendulum. 706 00:42:59,960 --> 00:43:00,636 That's good. 707 00:43:03,590 --> 00:43:09,020 So when Windows doesn't work, you can always rely on MATLAB. 708 00:43:09,020 --> 00:43:10,020 So what do we do? 709 00:43:12,540 --> 00:43:15,070 We want to code up this differential equation first. 710 00:43:15,070 --> 00:43:18,960 Whatever scheme we use, we need to have a function that 711 00:43:18,960 --> 00:43:21,730 computes f of u. 712 00:43:21,730 --> 00:43:24,180 The function should do like this. 713 00:43:24,180 --> 00:43:32,070 You give me u, the function returns f of u or du dt. 714 00:43:32,070 --> 00:43:39,740 I just call this function du dt equal to f of u. 715 00:43:39,740 --> 00:43:41,590 And what is u? 716 00:43:41,590 --> 00:43:50,270 u actually contains theta and theta dot. 717 00:43:50,270 --> 00:43:51,650 Let me do this. 718 00:43:51,650 --> 00:43:53,920 Theta equal to this first element of u. 719 00:43:56,600 --> 00:44:00,290 Theta dot is equal to the second element of u. 720 00:44:03,010 --> 00:44:06,720 I'm going to be computing theta double dot-- 721 00:44:06,720 --> 00:44:16,410 I'm just going to call this d dot-- is equal to minus g times 722 00:44:16,410 --> 00:44:25,250 theta divided by L. L is equal to minus g times theta divided 723 00:44:25,250 --> 00:44:25,846 by 5. 724 00:44:29,040 --> 00:44:32,180 So now what I need to return is du dt. 725 00:44:32,180 --> 00:44:34,920 And what is du dt now? 726 00:44:34,920 --> 00:44:38,430 Remember u is theta and theta dot. 727 00:44:38,430 --> 00:44:43,529 So du dt is actually theta dot and theta double dot. 728 00:44:49,330 --> 00:44:50,190 That's all. 729 00:44:50,190 --> 00:44:53,620 That function encodes the differential equation. 730 00:44:53,620 --> 00:44:56,880 If it will give me a u, I'm going to return a du dt. 731 00:44:56,880 --> 00:44:58,485 Here u is a [? matter ?] of 2. du 732 00:44:58,485 --> 00:45:02,140 dt is going to be a [? matter ?] of 2. 733 00:45:02,140 --> 00:45:04,183 Clear what this function does? 734 00:45:06,960 --> 00:45:16,530 Now, let's see what forward Euler does. 735 00:45:16,530 --> 00:45:19,420 Let's define dt equals 0.1, and we 736 00:45:19,420 --> 00:45:22,850 want to simulate this thing for how long? 737 00:45:22,850 --> 00:45:25,280 AUDIENCE: [INAUDIBLE]. 738 00:45:25,280 --> 00:45:26,373 QIQI WANG: 60 seconds. 739 00:45:26,373 --> 00:45:29,660 All right. 740 00:45:29,660 --> 00:45:33,990 I have to initialize a u0. 741 00:45:33,990 --> 00:45:36,280 So what u0 you want? 742 00:45:36,280 --> 00:45:40,731 I mean, it has to be a theta and d theta dt. 743 00:45:40,731 --> 00:45:41,230 Anybody? 744 00:45:41,230 --> 00:45:42,146 AUDIENCE: [INAUDIBLE]. 745 00:45:45,316 --> 00:45:49,980 QIQI WANG: If you look at this, when 746 00:45:49,980 --> 00:45:54,920 I am approximating sine theta as theta, 747 00:45:54,920 --> 00:45:58,060 is it in degrees or radius? 748 00:45:58,060 --> 00:45:59,930 It can be either, right? 749 00:45:59,930 --> 00:46:02,860 It can be either, so it doesn't matter. 750 00:46:02,860 --> 00:46:05,780 Oh, yes, it doesn't matter because it's a linear equation. 751 00:46:05,780 --> 00:46:10,060 It's a linear equation, which means if I change a unit, 752 00:46:10,060 --> 00:46:14,110 I'm just scaling everything by the unit conversion factor. 753 00:46:14,110 --> 00:46:18,000 And the linear equation still should be satisfied. 754 00:46:18,000 --> 00:46:20,670 That's like the definition of the linear equation. 755 00:46:20,670 --> 00:46:24,330 If you scale everything by some factor, it's still satisfied. 756 00:46:24,330 --> 00:46:29,850 So it can be either degrees or radius. 757 00:46:29,850 --> 00:46:31,185 AUDIENCE: 0. 758 00:46:31,185 --> 00:46:31,905 QIQI WANG: 0. 759 00:46:34,900 --> 00:46:35,440 And 1? 760 00:46:35,440 --> 00:46:36,326 AUDIENCE: [? It's ?] really boring. 761 00:46:36,326 --> 00:46:37,460 QIQI WANG: It's really boring. 762 00:46:37,460 --> 00:46:37,960 OK. 763 00:46:37,960 --> 00:46:39,740 AUDIENCE: [INAUDIBLE]. 764 00:46:39,740 --> 00:46:42,375 QIQI WANG: So this is my initial condition. 765 00:46:42,375 --> 00:46:45,100 And here's the forward Euler. 766 00:46:45,100 --> 00:46:47,220 Did you see this? 767 00:46:47,220 --> 00:46:49,270 Actually, I think I should do this. 768 00:46:49,270 --> 00:46:53,880 I should make another file so that you don't 769 00:46:53,880 --> 00:46:56,410 have to look at the bottom. 770 00:47:05,690 --> 00:47:06,315 forwardeuler.m. 771 00:47:10,090 --> 00:47:11,422 I have to open forwardeuler.m. 772 00:47:14,100 --> 00:47:21,060 So I'm going to be setting u0 1, 0. 773 00:47:21,060 --> 00:47:26,930 I'm going to set dt equal to 0.1, t equal to 60, 774 00:47:26,930 --> 00:47:36,200 and I'm going to forward t equal to 0, 0, dt t. 775 00:47:36,200 --> 00:47:39,190 So this is a loop. 776 00:47:39,190 --> 00:47:48,630 I'm going to compute du dt is equal to actually, 777 00:47:48,630 --> 00:47:51,215 I think pendulum of u0. 778 00:47:54,340 --> 00:47:57,900 And what is forward Euler? 779 00:47:57,900 --> 00:48:00,660 Forward Euler is my u at the next step 780 00:48:00,660 --> 00:48:06,334 is equal u0 plus u dt times dt. 781 00:48:11,270 --> 00:48:17,540 And when I go to the next step, u0 should be set to u. 782 00:48:17,540 --> 00:48:19,210 All right. 783 00:48:19,210 --> 00:48:25,650 And now what I want to do is I want to record the history. 784 00:48:25,650 --> 00:48:26,900 All right. 785 00:48:26,900 --> 00:48:29,060 I want to record the history. 786 00:48:29,060 --> 00:48:34,810 What I want to do is history is equal to an empty array, 787 00:48:34,810 --> 00:48:43,550 and the history is equal to history u0. 788 00:48:43,550 --> 00:48:45,716 I'm going to run this. 789 00:48:45,716 --> 00:48:48,220 It doesn't run. 790 00:48:48,220 --> 00:48:52,190 I'm going to forward Euler. 791 00:48:52,190 --> 00:48:55,610 It's done, and I have an array history, 792 00:48:55,610 --> 00:49:00,010 and I can plot history. 793 00:49:00,010 --> 00:49:03,860 It is done to give me a history of the numerical solutions. 794 00:49:07,480 --> 00:49:10,830 Does it look good? 795 00:49:10,830 --> 00:49:12,570 No. 796 00:49:12,570 --> 00:49:14,890 So this is forward Euler, and I'm 797 00:49:14,890 --> 00:49:19,230 using a pretty big time step. 798 00:49:19,230 --> 00:49:24,420 Let's see what can I do better if I use a smaller time step. 799 00:49:29,870 --> 00:49:33,550 You see, this is better. 800 00:49:33,550 --> 00:49:38,210 And a difficult solution should be an oscillating pendulum 801 00:49:38,210 --> 00:49:44,480 with no increase in the magnitude of the oscillation. 802 00:49:44,480 --> 00:49:51,300 And now can I do even better by having an even smaller delta t? 803 00:49:57,660 --> 00:49:58,310 All right. 804 00:49:58,310 --> 00:50:00,252 AUDIENCE: It's still working. 805 00:50:00,252 --> 00:50:01,293 QIQI WANG: Still working? 806 00:50:04,131 --> 00:50:06,701 Am I having a reasonable delta t? 807 00:50:06,701 --> 00:50:07,200 Yeah. 808 00:50:07,200 --> 00:50:20,480 It should be-- It actually takes a pretty absurd value 809 00:50:20,480 --> 00:50:25,990 of delta t to make this like not even super accurate. 810 00:50:25,990 --> 00:50:31,570 And you can visually see the empty array still growing. 811 00:50:31,570 --> 00:50:35,020 It is visually different from the analytical solution 812 00:50:35,020 --> 00:50:35,780 of the ODE. 813 00:50:40,180 --> 00:50:46,120 So let's go back and try to implement 814 00:50:46,120 --> 00:50:49,150 something that is better. 815 00:50:49,150 --> 00:50:52,119 Just name a scheme that is better than forward Euler. 816 00:50:52,119 --> 00:50:52,910 AUDIENCE: Midpoint. 817 00:50:52,910 --> 00:50:55,390 QIQI WANG: Midpoint. 818 00:50:55,390 --> 00:50:57,610 Sure, let's do midpoint. 819 00:50:57,610 --> 00:51:03,990 To avoid repetitive work, let me just rename the Forward 820 00:51:03,990 --> 00:51:07,300 Euler using Midpoint. 821 00:51:07,300 --> 00:51:09,700 And just copying the Forward Euler 822 00:51:09,700 --> 00:51:11,580 and to rename it as Midpoint. 823 00:51:11,580 --> 00:51:14,800 I'm going to load this guy. 824 00:51:14,800 --> 00:51:18,010 I'm opening midpoint.m. 825 00:51:18,010 --> 00:51:20,860 So this is Forward Euler, but I'm 826 00:51:20,860 --> 00:51:24,950 going to change it to Midpoint. 827 00:51:30,290 --> 00:51:32,227 How should midpoint be different. 828 00:51:35,640 --> 00:51:39,208 What is the difference between midpoint and forward Euler? 829 00:51:39,208 --> 00:51:40,144 AUDIENCE: [INAUDIBLE]. 830 00:51:44,830 --> 00:51:54,330 QIQI WANG: If you look at the formula, midpoint, 831 00:51:54,330 --> 00:52:00,770 instead of ui plus 1 equal to ui, plus delta t times fu, 832 00:52:00,770 --> 00:52:05,530 I have ui plus 1 equal to ui minus 1 plus 2 delta t times 833 00:52:05,530 --> 00:52:07,790 fu. 834 00:52:07,790 --> 00:52:12,630 So two things, instead of ui, I need ui minus 1. 835 00:52:12,630 --> 00:52:17,280 Instead of just f of ui, I need 2f of ui. 836 00:52:17,280 --> 00:52:20,340 The second one is easy the change. 837 00:52:20,340 --> 00:52:22,810 I can just multiply by 2. 838 00:52:22,810 --> 00:52:25,195 Now, how do I do about this? 839 00:52:27,810 --> 00:52:32,925 How can I get ui minus 1? 840 00:52:32,925 --> 00:52:33,923 AUDIENCE: [INAUDIBLE]. 841 00:52:41,408 --> 00:52:42,550 QIQI WANG: Good point. 842 00:52:42,550 --> 00:52:46,460 I need to do one step of forward Euler or something 843 00:52:46,460 --> 00:52:49,070 to get me one actually at time step first. 844 00:52:49,070 --> 00:52:50,220 So let me do this. 845 00:52:50,220 --> 00:52:52,350 Let me do one step of forward Euler. 846 00:52:52,350 --> 00:52:57,190 u dt equal to pendulum of u0. 847 00:52:57,190 --> 00:53:03,090 I'll just put a comment here, one step of forward Euler. 848 00:53:03,090 --> 00:53:10,840 And the u is equal to u0 plus du dt times dt. 849 00:53:10,840 --> 00:53:15,380 u0 is equal to u. 850 00:53:15,380 --> 00:53:16,197 Is that enough? 851 00:53:19,119 --> 00:53:20,100 AUDIENCE: [INAUDIBLE]. 852 00:53:23,306 --> 00:53:24,590 QIQI WANG: Exactly. 853 00:53:24,590 --> 00:53:27,920 I also need to store something else. 854 00:53:27,920 --> 00:53:29,540 I need to do this. 855 00:53:29,540 --> 00:53:33,340 See, I need to store an extra step. 856 00:53:33,340 --> 00:53:36,990 That is why I think now we get to why 857 00:53:36,990 --> 00:53:40,330 it's called a two-step scheme because you 858 00:53:40,330 --> 00:53:44,630 need to store two steps in your computer's memory. 859 00:53:44,630 --> 00:53:46,090 While for the forward Euler scheme, 860 00:53:46,090 --> 00:53:48,140 you only need to store u0. 861 00:53:48,140 --> 00:53:51,310 You only need to store one step in your memory. 862 00:53:51,310 --> 00:53:55,170 So implement of midpoint scheme, the two-step scheme, 863 00:53:55,170 --> 00:54:00,980 you need to store two time steps in your memory. 864 00:54:00,980 --> 00:54:03,100 You have to store u0 and u00. 865 00:54:03,100 --> 00:54:07,704 So here I'm just going to use u00, and u00 equal to u0, 866 00:54:07,704 --> 00:54:12,530 and u0 equal to u. 867 00:54:12,530 --> 00:54:16,218 You keep two previous time steps in memory. 868 00:54:19,730 --> 00:54:23,680 So this is going to work, and let me go back 869 00:54:23,680 --> 00:54:30,300 to 0.1 L of t, for which forward Euler did a pretty crazy job. 870 00:54:30,300 --> 00:54:35,030 It went over to something like 1,000 or something. 871 00:54:35,030 --> 00:54:36,640 Let's try midpoint. 872 00:54:40,140 --> 00:54:45,084 Let me Close All. 873 00:54:45,084 --> 00:54:47,056 Did I Close All? 874 00:54:51,900 --> 00:54:55,840 Midpoint, but I didn't plot it. 875 00:54:55,840 --> 00:54:58,070 Plot History. 876 00:55:03,006 --> 00:55:05,191 How about this? 877 00:55:05,191 --> 00:55:07,930 Is this much better? 878 00:55:07,930 --> 00:55:08,430 It's pretty. 879 00:55:11,040 --> 00:55:14,300 So let me actually plot it not just the history. 880 00:55:14,300 --> 00:55:17,890 Let me actually make the primes to be accurate. 881 00:55:17,890 --> 00:55:19,870 Let me make 0 dt t. 882 00:55:27,010 --> 00:55:31,510 Let me actually plot the x-axis as the real time instead 883 00:55:31,510 --> 00:55:32,960 of time steps. 884 00:55:32,960 --> 00:55:35,160 If I plot x-axis as unset, it would 885 00:55:35,160 --> 00:55:37,767 be like 1, 2, 3, 4, et cetera. [? It's ?] the number 886 00:55:37,767 --> 00:55:38,350 of time steps. 887 00:55:38,350 --> 00:55:45,480 Now, this is plotting [INAUDIBLE] solution 888 00:55:45,480 --> 00:55:46,940 of the ODE. 889 00:55:46,940 --> 00:55:52,050 So here we get visibly the same answer 890 00:55:52,050 --> 00:55:56,520 as an original solution, which is sinusoidal oscillations. 891 00:56:01,920 --> 00:56:06,256 We got half an hour left, but I want to take a short break now 892 00:56:06,256 --> 00:56:09,420 to get our minds back. 893 00:56:09,420 --> 00:56:12,316 And next what we're going to be doing is [INAUDIBLE]. 894 00:56:15,310 --> 00:56:18,715 I'll see how it goes. 895 00:56:18,715 --> 00:56:19,340 Let's continue. 896 00:56:19,340 --> 00:56:23,250 Is everybody back? 897 00:56:23,250 --> 00:56:25,810 Before I try your schemes, let me just 898 00:56:25,810 --> 00:56:30,600 show how easy it is to change what 899 00:56:30,600 --> 00:56:34,550 we did with the linear differential equation, 900 00:56:34,550 --> 00:56:35,810 change it to a non-linear one. 901 00:56:38,246 --> 00:56:39,620 To change it to a non-linear one, 902 00:56:39,620 --> 00:56:43,570 we just keep the sine theta here. 903 00:56:43,570 --> 00:56:48,320 Just replace this with sine of theta. 904 00:56:48,320 --> 00:56:51,460 As soon as we change this to sine of theta, 905 00:56:51,460 --> 00:56:53,140 we lose the analytical solution. 906 00:56:56,040 --> 00:56:58,000 I don't know how to solve this anymore. 907 00:56:58,000 --> 00:57:00,010 Do you? 908 00:57:00,010 --> 00:57:05,075 But once we have a computer, what I need to change 909 00:57:05,075 --> 00:57:06,875 is let me open my pendulum.m. 910 00:57:10,740 --> 00:57:11,680 What do I need to do? 911 00:57:16,290 --> 00:57:20,890 I just need to replace theta with sine of theta. 912 00:57:20,890 --> 00:57:23,700 And let me do this again. 913 00:57:23,700 --> 00:57:25,720 Midpoint. 914 00:57:25,720 --> 00:57:30,730 I get a solution as easy as I'm solving the linear equations. 915 00:57:30,730 --> 00:57:39,975 And when I plot it, I want to plot this with history. 916 00:57:43,002 --> 00:57:44,960 Here I get the solution to non-linear equation. 917 00:57:44,960 --> 00:57:49,570 I mean, it may look similar, but let me zoom in a little bit 918 00:57:49,570 --> 00:57:51,200 to see. 919 00:57:51,200 --> 00:57:54,240 Oh, man, MATLAB doesn't listen to me. 920 00:57:54,240 --> 00:57:56,950 I'm going to the x lim, which is kind of a scale 921 00:57:56,950 --> 00:58:01,800 to x limit to 0 and 5 here. 922 00:58:01,800 --> 00:58:05,610 And well, it doesn't really show a lot of difference, 923 00:58:05,610 --> 00:58:09,505 but let me change the initial condition. 924 00:58:12,740 --> 00:58:17,270 Open midpoint to m. 925 00:58:17,270 --> 00:58:21,140 Let me change the initial condition to how much? 926 00:58:21,140 --> 00:58:24,960 1 is like 60 degrees. 927 00:58:24,960 --> 00:58:27,640 AUDIENCE: [INAUDIBLE]. 928 00:58:27,640 --> 00:58:30,740 QIQI WANG: [? Reads ?] very close to pi, so you get 929 00:58:30,740 --> 00:58:32,798 like-- so for a let's try it. 930 00:58:36,144 --> 00:58:39,970 Let me plot it again. 931 00:58:39,970 --> 00:58:41,010 Oh, yeah. 932 00:58:41,010 --> 00:58:53,590 So we get [INAUDIBLE] So it will stay close to the top 933 00:58:53,590 --> 00:58:57,646 for a while and then drop down and stay at [? 0 ?] for a while 934 00:58:57,646 --> 00:59:00,840 and then [? conditions ?] the other part of the chart. 935 00:59:00,840 --> 00:59:03,526 To get over the middle of the area [INAUDIBLE]. 936 00:59:06,560 --> 00:59:09,920 So you get some very nonlinear behaviors 937 00:59:09,920 --> 00:59:12,600 where you have a big initial [? areas ?] as easy 938 00:59:12,600 --> 00:59:14,930 as you're solving a linear equation. 939 00:59:14,930 --> 00:59:18,740 And that is why we use numerical methods. 940 00:59:18,740 --> 00:59:21,610 This is when we don't have an easy way 941 00:59:21,610 --> 00:59:24,110 to get analytical solutions, and it is certainly 942 00:59:24,110 --> 00:59:27,130 true for most of the engineering problems we are dealing with. 943 00:59:29,680 --> 00:59:32,170 And it's quite rare you can find any of the solutions. 944 00:59:32,170 --> 00:59:35,485 And most of the time you have to use numerical solutions. 945 00:59:38,370 --> 00:59:42,510 So let's go back and start to try-- we did forward Euler. 946 00:59:42,510 --> 00:59:44,290 We did midpoint rule. 947 00:59:44,290 --> 00:59:48,290 Let's try best implicit one-step scheme. 948 00:59:48,290 --> 00:59:53,480 And I think I remember what it is. 949 00:59:53,480 --> 00:59:55,190 Please correct me if I'm wrong. 950 00:59:55,190 --> 01:00:02,637 ui plus 1 equal to ui plus 1/2 of delta t ui prime plus 1/2 951 01:00:02,637 --> 01:00:05,811 of delta t ui plus 1 prime. 952 01:00:09,210 --> 01:00:15,040 This is the one we just derived from writing down all the rows, 953 01:00:15,040 --> 01:00:19,400 having all the columns aligned and canceled three 954 01:00:19,400 --> 01:00:22,630 of the most important terms. 955 01:00:22,630 --> 01:00:23,410 And we get this. 956 01:00:23,410 --> 01:00:30,420 We get x being 1, y being 1/2 L of t, z being 1/2 L of t. 957 01:00:30,420 --> 01:00:34,330 Now, let's try this, and for this, I 958 01:00:34,330 --> 01:00:38,140 have to say for now I am going to go back 959 01:00:38,140 --> 01:00:42,150 cowardly to a linear equation. 960 01:00:42,150 --> 01:00:46,670 And we are going to be talking about how 961 01:00:46,670 --> 01:00:53,526 to do the same with nonlinear equations. 962 01:00:53,526 --> 01:00:58,574 I think three lectures from now or something like that, 963 01:00:58,574 --> 01:00:59,990 we are going to be using something 964 01:00:59,990 --> 01:01:01,710 called the Newton-Raphson. 965 01:01:01,710 --> 01:01:06,210 But for now, let's stick to linear equations. 966 01:01:09,350 --> 01:01:23,240 What we have is we have u prime is equal to something times ui. 967 01:01:23,240 --> 01:01:24,870 ui is 2 by 1 vector. 968 01:01:27,930 --> 01:01:40,940 ui prime is going to be 1, 0, 0, minus g over L times ui. 969 01:01:43,910 --> 01:01:45,940 Why is that true? 970 01:01:45,940 --> 01:01:49,000 You just need to look at the code. 971 01:01:49,000 --> 01:01:54,130 du dt equal to theta dot, which is u2. 972 01:01:54,130 --> 01:02:04,880 And theta d dot, which is u1 times minus g over L. 973 01:02:04,880 --> 01:02:16,530 This is the same as saying du dt equal to u2, u1 times minus g 974 01:02:16,530 --> 01:02:22,400 over L. I can basically rephrase all of this by this. 975 01:02:28,450 --> 01:02:30,710 Or all the stuff I just commented out, 976 01:02:30,710 --> 01:02:32,690 made green is just this. 977 01:02:35,870 --> 01:02:38,120 Agreed? 978 01:02:38,120 --> 01:02:43,990 And this is just saying that du dt is 979 01:02:43,990 --> 01:02:47,960 equal to this matrix times u. 980 01:02:47,960 --> 01:02:50,540 And this is going back and forth between the formula 981 01:02:50,540 --> 01:02:55,970 and the code, the formula and the code. 982 01:02:55,970 --> 01:02:58,713 Doing this is the same as a matrix multiplication. 983 01:03:01,680 --> 01:03:07,450 Now, the same thing happens for ui plus 1 prime. 984 01:03:15,762 --> 01:03:18,570 Agreed? 985 01:03:18,570 --> 01:03:23,120 So now if you plot this into the equation, 986 01:03:23,120 --> 01:03:27,280 we removed all the primes and replaced all the primes 987 01:03:27,280 --> 01:03:30,770 with ui and ui plus 1. 988 01:03:30,770 --> 01:03:33,360 We can actually collect the terms. 989 01:03:33,360 --> 01:03:41,890 We can say ui plus 1 equal to-- I'm 990 01:03:41,890 --> 01:03:46,740 going to use ui plus 1/2 delta t times ui prime, which 991 01:03:46,740 --> 01:03:56,820 I am going to write like this plus-- and switching to red, 992 01:03:56,820 --> 01:04:05,350 and you're going to see why I am doing this-- ui plus 1. 993 01:04:07,819 --> 01:04:08,860 This is really important. 994 01:04:08,860 --> 01:04:11,240 It is ui plus 1. 995 01:04:11,240 --> 01:04:13,101 Why am I using different colors? 996 01:04:16,050 --> 01:04:17,960 What is the difference between the blue terms 997 01:04:17,960 --> 01:04:18,800 and the red terms? 998 01:04:18,800 --> 01:04:20,210 AUDIENCE: [INAUDIBLE]. 999 01:04:20,210 --> 01:04:21,922 QIQI WANG: Yes. 1000 01:04:21,922 --> 01:04:23,846 AUDIENCE: [INAUDIBLE]. 1001 01:04:23,846 --> 01:04:26,630 QIQI WANG: The current step versus next step. 1002 01:04:26,630 --> 01:04:28,060 The known versus unknown. 1003 01:04:30,930 --> 01:04:34,080 When I have the i-th step when I want to go to the next step, 1004 01:04:34,080 --> 01:04:36,300 ui plus 1 is unknown. 1005 01:04:36,300 --> 01:04:38,670 ui is known. 1006 01:04:38,670 --> 01:04:39,170 All right. 1007 01:04:39,170 --> 01:04:54,080 So what I can do is I can rearrange and say 1008 01:04:54,080 --> 01:05:02,801 that something times ui plus 1 is equal to something times ui. 1009 01:05:02,801 --> 01:05:05,940 Can somebody tell me what these two matrices are? 1010 01:05:05,940 --> 01:05:08,510 If I move all the red terms to the left 1011 01:05:08,510 --> 01:05:14,580 and let me actually write it as blue because [? there-- ?] 1012 01:05:14,580 --> 01:05:19,160 and all the blue terms to the right, what do I get? 1013 01:05:19,160 --> 01:05:20,390 What is the red matrix? 1014 01:05:20,390 --> 01:05:21,445 What is the blue matrix? 1015 01:05:35,025 --> 01:05:36,480 The red matrix. 1016 01:05:36,480 --> 01:05:38,088 AUDIENCE: [INAUDIBLE] 1017 01:05:38,088 --> 01:05:39,884 QIQI WANG: 1. 1018 01:05:39,884 --> 01:05:40,878 AUDIENCE: [INAUDIBLE]. 1019 01:06:04,237 --> 01:06:05,820 QIQI WANG: Great. 1020 01:06:05,820 --> 01:06:10,950 The two 1's are because ui can be represented 1021 01:06:10,950 --> 01:06:14,060 as identity matrix times ui. 1022 01:06:14,060 --> 01:06:17,890 Identity matrix is 1, 0, 0, 1. 1023 01:06:17,890 --> 01:06:21,030 Now, these two terms are basically this matrix times 1024 01:06:21,030 --> 01:06:21,880 delta t over 2. 1025 01:06:26,060 --> 01:06:27,580 Any questions about this matrix? 1026 01:06:27,580 --> 01:06:29,896 I mean, this is just combining these two terms. 1027 01:06:33,680 --> 01:06:34,620 How about the red one? 1028 01:06:37,560 --> 01:06:40,928 AUDIENCE: [INAUDIBLE]. 1029 01:06:40,928 --> 01:06:42,290 QIQI WANG: 1. 1030 01:06:42,290 --> 01:06:43,652 AUDIENCE: [INAUDIBLE] 1031 01:06:43,652 --> 01:06:46,580 QIQI WANG: Negative delta t over 2, yeah. 1032 01:06:46,580 --> 01:06:47,340 Exactly. 1033 01:06:47,340 --> 01:06:50,690 It is basically taking this matrix takes the negative sign. 1034 01:06:50,690 --> 01:06:51,340 And here? 1035 01:06:51,340 --> 01:06:52,280 AUDIENCE: [INAUDIBLE]. 1036 01:06:55,570 --> 01:06:56,520 QIQI WANG: 1. 1037 01:06:56,520 --> 01:06:57,740 Great. 1038 01:06:57,740 --> 01:07:02,100 So to implement the scheme, we need this. 1039 01:07:02,100 --> 01:07:06,870 We need to basically multiply this matrix with ui 1040 01:07:06,870 --> 01:07:09,430 and invert this matrix. 1041 01:07:09,430 --> 01:07:11,450 All right. 1042 01:07:11,450 --> 01:07:13,680 Clear? 1043 01:07:13,680 --> 01:07:17,990 I'm going to make another file. 1044 01:07:17,990 --> 01:07:26,330 I'm going to call it-- this scheme by the way, 1045 01:07:26,330 --> 01:07:29,960 the best one step implicit schemes you guys invented 1046 01:07:29,960 --> 01:07:31,440 is called the trapezoidal rule. 1047 01:07:36,310 --> 01:07:37,460 Not tx, t. 1048 01:07:37,460 --> 01:07:43,665 Sorry-- dot m. 1049 01:07:43,665 --> 01:07:44,165 Yes. 1050 01:07:53,260 --> 01:07:54,954 Oh, I should be copying. 1051 01:08:02,560 --> 01:08:04,410 this. 1052 01:08:04,410 --> 01:08:06,865 Trapezoidal rule is a one-step scheme, 1053 01:08:06,865 --> 01:08:08,790 is a single-step scheme, so I don't 1054 01:08:08,790 --> 01:08:12,050 need to do any forward Euler or anything like that. 1055 01:08:12,050 --> 01:08:13,080 All right. 1056 01:08:13,080 --> 01:08:16,250 And I'm going to be doing t equal to 0 1057 01:08:16,250 --> 01:08:22,680 to dt equal to t. [? I'm ?] going to type m first. 1058 01:08:22,680 --> 01:08:24,790 What I'm going to do is I'm going to talk 1059 01:08:24,790 --> 01:08:27,460 to construct two matrices. 1060 01:08:27,460 --> 01:08:33,990 One matrix, I call it A. Let me just 1061 01:08:33,990 --> 01:08:39,368 call it red matrix and the blue matrix. 1062 01:08:44,290 --> 01:08:51,479 The red matrix is 1 minus delta t, delta t over 2-- 1063 01:08:51,479 --> 01:08:53,779 I'm just going to type it. 1064 01:08:53,779 --> 01:08:58,600 1 delta t over 2, delta t over 2. 1065 01:08:58,600 --> 01:09:00,390 I think that I get a minus here. 1066 01:09:00,390 --> 01:09:09,260 That will be over 2 times g over L. and 1. 1067 01:09:09,260 --> 01:09:12,439 The blue matrix is roughly the same, 1068 01:09:12,439 --> 01:09:17,810 but there is a minus sign here, and there is a plus sign here. 1069 01:09:17,810 --> 01:09:18,723 All right. 1070 01:09:24,569 --> 01:09:26,281 Then what do I do next? 1071 01:09:31,000 --> 01:09:33,790 u is equal to what? 1072 01:09:37,290 --> 01:09:40,526 AUDIENCE: [INAUDIBLE]. 1073 01:09:40,526 --> 01:09:51,720 QIQI WANG: Inverse of the red matrix times blue matrix times 1074 01:09:51,720 --> 01:09:53,149 u0. 1075 01:09:53,149 --> 01:09:55,890 So actually here I need to do u0 transpose 1076 01:09:55,890 --> 01:09:57,900 because u0 is a row matrix. 1077 01:09:57,900 --> 01:10:00,280 I need to make it a column matrix, 1078 01:10:00,280 --> 01:10:05,495 and the whole thing has to be transposed again. 1079 01:10:08,650 --> 01:10:09,520 All right. 1080 01:10:09,520 --> 01:10:12,300 And I'm going to be doing history 1081 01:10:12,300 --> 01:10:19,172 is equal to history is 0 and u0 equal to [? u. ?] 1082 01:10:23,620 --> 01:10:24,860 Is it clear what I did? 1083 01:10:28,130 --> 01:10:30,240 Anybody have any comments on this? 1084 01:10:30,240 --> 01:10:34,073 Is it a good way of doing matrix inverse? 1085 01:10:34,073 --> 01:10:36,378 AUDIENCE: [INAUDIBLE]. 1086 01:10:36,378 --> 01:10:39,444 QIQI WANG: Yes. u is ui plus 1. 1087 01:10:39,444 --> 01:10:40,398 AUDIENCE: [INAUDIBLE]. 1088 01:10:44,900 --> 01:10:46,654 QIQI WANG: Yes. 1089 01:10:46,654 --> 01:10:47,598 AUDIENCE: [INAUDIBLE]. 1090 01:10:52,790 --> 01:10:54,080 QIQI WANG: Good point. 1091 01:10:54,080 --> 01:10:55,080 AUDIENCE: [INAUDIBLE]. 1092 01:10:59,890 --> 01:11:06,580 QIQI WANG: Inverse red matrix is equal to inverse red matrix. 1093 01:11:06,580 --> 01:11:07,580 AUDIENCE: [INAUDIBLE]. 1094 01:11:11,014 --> 01:11:11,680 QIQI WANG: Sure. 1095 01:11:11,680 --> 01:11:13,480 I can also combine blue matrix. 1096 01:11:13,480 --> 01:11:14,950 AUDIENCE: [INAUDIBLE]. 1097 01:11:14,950 --> 01:11:15,920 QIQI WANG: Sure. 1098 01:11:15,920 --> 01:11:16,920 A good suggestion. 1099 01:11:22,370 --> 01:11:28,630 I'm just going to be u equal to a times-- all right. 1100 01:11:28,630 --> 01:11:32,140 And let me actually make u0 a column matrix 1101 01:11:32,140 --> 01:11:34,190 by replacing this with this. 1102 01:11:34,190 --> 01:11:35,510 All right. 1103 01:11:35,510 --> 01:11:39,790 So I don't need to do all of this. 1104 01:11:39,790 --> 01:11:44,870 Here, instead of when I recall the history, I transpose it. 1105 01:11:44,870 --> 01:11:45,370 Great. 1106 01:11:45,370 --> 01:11:48,378 Let's see how it works. 1107 01:11:48,378 --> 01:11:48,878 Trapezoidal. 1108 01:11:56,091 --> 01:11:56,590 Yeah. 1109 01:11:59,100 --> 01:12:02,524 Let's plot this again. 1110 01:12:02,524 --> 01:12:04,666 It works perfectly. 1111 01:12:04,666 --> 01:12:06,300 All right. 1112 01:12:06,300 --> 01:12:08,210 So good job, guys. 1113 01:12:08,210 --> 01:12:11,790 The scheme you came up with works. 1114 01:12:11,790 --> 01:12:15,260 That's it for today. 1115 01:12:15,260 --> 01:12:19,030 And so next class, what I want to do 1116 01:12:19,030 --> 01:12:21,080 is bring your own computer because we're 1117 01:12:21,080 --> 01:12:24,730 going to be doing something-- next class, 1118 01:12:24,730 --> 01:12:29,430 not I'm going to code up things like this. 1119 01:12:29,430 --> 01:12:35,050 Instead I'm going to ask you to code up something like this, 1120 01:12:35,050 --> 01:12:37,320 and you are going to be working class 1121 01:12:37,320 --> 01:12:42,500 and hopefully, show things on this screen.