1 00:00:00,090 --> 00:00:01,260 Welcome to my office. 2 00:00:01,260 --> 00:00:04,410 I'll now give you a look into 3 00:00:04,410 --> 00:00:07,380 how I usually work when I'm interested in 4 00:00:07,380 --> 00:00:10,290 following the behavior of a dynamical system 5 00:00:10,290 --> 00:00:12,030 and I need to rely on 6 00:00:12,030 --> 00:00:15,243 numerical integration or numerical methods to do so. 7 00:00:17,130 --> 00:00:21,780 I will start by introducing 8 00:00:21,780 --> 00:00:26,760 a slightly more general Lotka-Volterra model 9 00:00:26,760 --> 00:00:29,283 than what we've seen so far. 10 00:00:32,880 --> 00:00:36,450 The reason being that the two-species Lotka-Volterra model, 11 00:00:36,450 --> 00:00:38,490 we've been able to study its behavior, 12 00:00:38,490 --> 00:00:41,280 look at its fixed points, flow diagrams 13 00:00:41,280 --> 00:00:43,800 with simple pen and paper, using null clients 14 00:00:43,800 --> 00:00:47,400 and solving for our derivatives equal to zero. 15 00:00:47,400 --> 00:00:50,487 In the case of three-species Lotka-Volterra, 16 00:00:50,487 --> 00:00:51,320 things get a little crazier. 17 00:00:51,320 --> 00:00:54,750 The algebra gets more complex, it gets harder to visualize. 18 00:00:54,750 --> 00:00:56,730 Four-species, even more so. 19 00:00:56,730 --> 00:00:59,883 And so numerical methods are going to be very useful. 20 00:01:01,770 --> 00:01:04,050 We're gonna use a slightly different form. 21 00:01:04,050 --> 00:01:07,500 So this is what's often called a generalized Lotka-Volterra. 22 00:01:07,500 --> 00:01:11,160 It has those logistic-looking terms 23 00:01:11,160 --> 00:01:13,890 like X I times one minus X I, 24 00:01:13,890 --> 00:01:17,640 but it's also defined over a set N of species, 25 00:01:17,640 --> 00:01:21,990 so we would have X one, x two, up to X N for X N species, 26 00:01:21,990 --> 00:01:26,400 and all interactions would be encoded in this A matrix. 27 00:01:26,400 --> 00:01:29,160 So the interaction A I J would represent 28 00:01:29,160 --> 00:01:32,703 the impact of species J on species I. 29 00:01:33,660 --> 00:01:36,210 So for example, in the case of a three-species model, 30 00:01:36,210 --> 00:01:39,390 we would have nine parameters for interaction. 31 00:01:39,390 --> 00:01:40,860 It's important to realize that 32 00:01:40,860 --> 00:01:44,550 because of the negative term in the first equation, 33 00:01:44,550 --> 00:01:49,140 so the fact that the X I times X J term Is a negative term, 34 00:01:49,140 --> 00:01:53,580 it means that by default a positive parameter, 35 00:01:53,580 --> 00:01:56,130 a positive element in the matrix A 36 00:01:56,130 --> 00:01:58,080 would result in a decrease, 37 00:01:58,080 --> 00:02:01,740 so a negative impact on the derivative of X I. 38 00:02:01,740 --> 00:02:05,310 And so positive parameter means something like predation, 39 00:02:05,310 --> 00:02:09,600 So A one two means the species two is hurting species one 40 00:02:09,600 --> 00:02:11,070 if A one two is positive, 41 00:02:11,070 --> 00:02:13,413 and a negative parameter means the opposite. 42 00:02:15,600 --> 00:02:17,880 We're gonna rely on numerical integration. 43 00:02:17,880 --> 00:02:20,700 If I go back, you'll see that in this matrix 44 00:02:20,700 --> 00:02:22,800 throughout this video, I'll keep everything fixed, 45 00:02:22,800 --> 00:02:25,980 but only vary this A three one element 46 00:02:25,980 --> 00:02:27,993 which I'm parameterizing with alpha. 47 00:02:28,920 --> 00:02:31,980 And so you can think of the middle row, species two, 48 00:02:31,980 --> 00:02:36,980 as being a species that predates on the other two 49 00:02:38,460 --> 00:02:40,740 because both species one and two have a positive 50 00:02:40,740 --> 00:02:42,540 because of the negative element here. 51 00:02:42,540 --> 00:02:44,520 It's hard to wrap your arm around that. 52 00:02:44,520 --> 00:02:46,470 Positive impact on species two, 53 00:02:46,470 --> 00:02:47,760 meaning it's sort of a predator. 54 00:02:47,760 --> 00:02:49,740 And by tuning species three, 55 00:02:49,740 --> 00:02:51,570 we're deciding whether it's above. 56 00:02:51,570 --> 00:02:53,190 By tuning alpha, I mean, 57 00:02:53,190 --> 00:02:54,720 we're deciding whether species three 58 00:02:54,720 --> 00:02:57,750 is above or below species one in the trophic level. 59 00:02:57,750 --> 00:02:59,340 That's one way to think about this. 60 00:02:59,340 --> 00:03:00,960 But we don't have to bother too much 61 00:03:00,960 --> 00:03:03,180 about the reality behind this model. 62 00:03:03,180 --> 00:03:04,710 We're really interested in 63 00:03:04,710 --> 00:03:06,510 following its behavior through time, 64 00:03:06,510 --> 00:03:08,970 and the only way to do it well 65 00:03:08,970 --> 00:03:11,073 is to rely on numerical integration. 66 00:03:11,940 --> 00:03:14,370 So from a previous video, you might remember 67 00:03:14,370 --> 00:03:16,650 that the goal of numerical integration 68 00:03:16,650 --> 00:03:20,280 is to capture how X of T is gonna change 69 00:03:20,280 --> 00:03:22,410 and if I make small steps in time, 70 00:03:22,410 --> 00:03:26,013 so if I go to S T plus H forward in time. 71 00:03:26,880 --> 00:03:29,130 And from Taylor series, we know that 72 00:03:29,130 --> 00:03:32,970 value X at time T plus H is gonna depend on the derivative. 73 00:03:32,970 --> 00:03:34,560 Here I'm using primes to denote 74 00:03:34,560 --> 00:03:39,560 the derivative of the solution X of T, 75 00:03:39,750 --> 00:03:42,690 but it's also gonna depend on the curvature X prime prime, 76 00:03:42,690 --> 00:03:45,060 on the third derivative, X prime prime prime, 77 00:03:45,060 --> 00:03:47,040 and so on and so forth. 78 00:03:47,040 --> 00:03:51,210 But all we have is the current value at that time T, X of T, 79 00:03:51,210 --> 00:03:55,320 and a function that gives me my time derivative. 80 00:03:55,320 --> 00:03:56,850 So remember that whether it is, 81 00:03:56,850 --> 00:04:00,660 D X over D T, X that or X prime, that doesn't matter. 82 00:04:00,660 --> 00:04:03,780 That derivative is equal to some function 83 00:04:03,780 --> 00:04:06,243 of X of T and maybe even of time. 84 00:04:07,320 --> 00:04:10,020 So we know how to calculate X prime 85 00:04:10,020 --> 00:04:12,462 because we have this differential equation, 86 00:04:12,462 --> 00:04:14,340 we have that function that gives us X prime. 87 00:04:14,340 --> 00:04:16,607 So Euler's method is going to say 88 00:04:16,607 --> 00:04:18,720 I'm just gonna take what I have, 89 00:04:18,720 --> 00:04:20,940 I'm gonna assume that I'm at point X of T, 90 00:04:20,940 --> 00:04:24,810 I'm following a slope for a length of time H 91 00:04:24,810 --> 00:04:27,360 and I'm gonna know where I am after that time step H, 92 00:04:27,360 --> 00:04:30,300 and I'm gonna ignore the higher order derivatives. 93 00:04:30,300 --> 00:04:33,750 So the local error, as we've said before, is of order two 94 00:04:33,750 --> 00:04:36,180 because we're ignoring this H squared term 95 00:04:36,180 --> 00:04:38,100 which is the most significant 96 00:04:38,100 --> 00:04:40,800 since we're assuming H to be small. 97 00:04:40,800 --> 00:04:43,620 So the local error is gonna be an H square. 98 00:04:43,620 --> 00:04:45,820 The global error over, you know, 99 00:04:45,820 --> 00:04:50,760 a capital T times series that we're interested in 100 00:04:50,760 --> 00:04:53,820 is gonna take capital T over H steps 101 00:04:53,820 --> 00:04:58,230 and therefore capital T over H time H square error, 102 00:04:58,230 --> 00:05:01,740 which gives us the capital T times H global error. 103 00:05:01,740 --> 00:05:03,240 The global error is always gonna be 104 00:05:03,240 --> 00:05:07,350 one order smaller than the local error 105 00:05:07,350 --> 00:05:10,173 since the smaller H is, the more steps we have to do. 106 00:05:11,790 --> 00:05:14,100 Heun's method we've also discussed, 107 00:05:14,100 --> 00:05:17,390 as you'll see in the assignment, 108 00:05:17,390 --> 00:05:19,920 it combines multiple derivatives 109 00:05:19,920 --> 00:05:22,740 to evaluation of the derivative in a smart way 110 00:05:22,740 --> 00:05:24,360 such that it can approximate 111 00:05:24,360 --> 00:05:26,460 or actually capture, I should say, 112 00:05:26,460 --> 00:05:27,840 the second derivative 113 00:05:27,840 --> 00:05:31,500 without having an explicit equation for it. 114 00:05:31,500 --> 00:05:34,650 So that's a really interesting method. 115 00:05:34,650 --> 00:05:36,540 So we're gonna use some of those. 116 00:05:36,540 --> 00:05:39,030 I'll rely on some more advanced integrator 117 00:05:39,030 --> 00:05:41,970 to give us an idea of a ground truth. 118 00:05:41,970 --> 00:05:43,680 There is no, like, ground truth here, 119 00:05:43,680 --> 00:05:45,840 no perfect solution for X of T, 120 00:05:45,840 --> 00:05:48,810 but we can rely on more advanced algorithms 121 00:05:48,810 --> 00:05:52,803 just to compare the error done by those simpler approaches. 122 00:05:53,910 --> 00:05:55,350 So we're gonna look at 123 00:05:55,350 --> 00:06:00,350 this generalized Lotka-Volterra model in greater detail. 124 00:06:08,229 --> 00:06:13,229 To do so, I will take you to my usual workspace. 125 00:06:14,280 --> 00:06:17,550 So you should now be seeing my codes in the top left. 126 00:06:17,550 --> 00:06:19,510 So this is in Octave, 127 00:06:19,510 --> 00:06:23,981 which is an open-source free alternative to Matlab. 128 00:06:23,981 --> 00:06:25,950 If you want to use this code, 129 00:06:25,950 --> 00:06:28,800 I'll be sharing it on the Blackboard site. 130 00:06:28,800 --> 00:06:30,930 There are some special functions, 131 00:06:30,930 --> 00:06:33,900 like calls to this O D E four five integrator, 132 00:06:33,900 --> 00:06:35,820 which you might have to play with a little bit, 133 00:06:35,820 --> 00:06:38,460 but it should almost work out of the box 134 00:06:38,460 --> 00:06:42,090 for your Matlab software. 135 00:06:42,090 --> 00:06:43,890 In the bottom left, what you have is 136 00:06:43,890 --> 00:06:46,530 this command line version of Octave 137 00:06:46,530 --> 00:06:49,413 where I'll be able to simply call my different functions. 138 00:06:51,600 --> 00:06:54,750 So very briefly here, how these codes are organized is 139 00:06:54,750 --> 00:06:57,840 you have a main code, which is "compare integrators" 140 00:06:57,840 --> 00:06:59,520 because that's what I want to do. 141 00:06:59,520 --> 00:07:03,630 I'm gonna be taking this time series, capital T, 142 00:07:03,630 --> 00:07:06,660 which is gonna go from time zero to time max T 143 00:07:06,660 --> 00:07:11,100 which is an argument of my function call, 144 00:07:11,100 --> 00:07:13,530 and I'm gonna do so in time steps of H. 145 00:07:13,530 --> 00:07:16,527 So this, you know, in Python you might use NumPy 146 00:07:16,527 --> 00:07:18,270 and this would be something like 147 00:07:18,270 --> 00:07:20,650 a NumPy arrange function call 148 00:07:21,602 --> 00:07:24,630 to create my different time steps. 149 00:07:24,630 --> 00:07:27,720 And I also have initial conditions X zero, 150 00:07:27,720 --> 00:07:30,990 so that's the starting point for my time series. 151 00:07:30,990 --> 00:07:33,690 And with that I can call O D E four five 152 00:07:33,690 --> 00:07:36,120 which is an advanced Runge Kutta 153 00:07:36,120 --> 00:07:38,340 with adaptive time step integrator 154 00:07:38,340 --> 00:07:39,990 that we're gonna assume is 155 00:07:39,990 --> 00:07:42,570 close enough to the the real trajectory, 156 00:07:42,570 --> 00:07:43,530 and we're gonna compare 157 00:07:43,530 --> 00:07:46,233 what the orders method does against it. 158 00:07:47,370 --> 00:07:50,670 So this line here is creating my time series of T 159 00:07:50,670 --> 00:07:54,150 and Y, the state of my system over time. 160 00:07:54,150 --> 00:07:55,320 And I'm doing the same thing 161 00:07:55,320 --> 00:07:57,120 with my own Euler integrator 162 00:07:57,120 --> 00:07:59,673 using time step H to compare them. 163 00:08:01,980 --> 00:08:03,483 We won't need this for now. 164 00:08:05,340 --> 00:08:07,200 And we'll be comparing, 165 00:08:07,200 --> 00:08:10,000 actually, let me remove all of this for now. 166 00:08:11,487 --> 00:08:14,403 Remove this, keep it as simple as possible. 167 00:08:23,610 --> 00:08:27,180 Okay, so let's look at what these other functions look like. 168 00:08:27,180 --> 00:08:31,279 So O D E four five, I have to pass it my set of equation. 169 00:08:31,279 --> 00:08:32,520 It's always good to define your equation 170 00:08:32,520 --> 00:08:34,440 in a separate function. 171 00:08:34,440 --> 00:08:37,680 The reason for that is if I wanted to use a different model, 172 00:08:37,680 --> 00:08:38,910 I wouldn't have to change 173 00:08:38,910 --> 00:08:41,010 my comparing integrator function here. 174 00:08:41,010 --> 00:08:43,590 I would just have to give O D E four five 175 00:08:43,590 --> 00:08:46,020 a different function to integrate. 176 00:08:46,020 --> 00:08:49,350 So if I go in my predator-prey system, 177 00:08:49,350 --> 00:08:52,230 I've hard-coded my matrix A 178 00:08:52,230 --> 00:08:55,110 because I said the only parameter that I'll be changing 179 00:08:55,110 --> 00:08:58,173 is this alpha, or here lower case A, 180 00:09:00,630 --> 00:09:02,940 and this is in matricial form 181 00:09:02,940 --> 00:09:07,480 calculating the three equations for D X one over D T, 182 00:09:07,480 --> 00:09:10,020 D X two over D T, and D X three over D T 183 00:09:10,020 --> 00:09:12,395 that you saw in the previous slide. 184 00:09:12,395 --> 00:09:15,690 It's doing the exact same form in matricial operations. 185 00:09:15,690 --> 00:09:17,750 The reason for that being that 186 00:09:17,750 --> 00:09:20,070 whenever you use Matlab, remember that the "Mat" in Matlab 187 00:09:20,070 --> 00:09:22,350 is for matrix not mathematics 188 00:09:22,350 --> 00:09:26,190 and therefore it is optimized for matricial operations. 189 00:09:26,190 --> 00:09:27,840 So this makes it a little faster than 190 00:09:27,840 --> 00:09:30,693 if I have three lines for the three different equations. 191 00:09:31,980 --> 00:09:35,010 The other function that I have is for my Euler integrator, 192 00:09:35,010 --> 00:09:36,990 so this Euler pred-prey function, 193 00:09:36,990 --> 00:09:39,270 and if we go in there, we'll see that 194 00:09:39,270 --> 00:09:42,360 I have the same hard-coding of my A matrix 195 00:09:42,360 --> 00:09:44,283 with my interaction parameter A, 196 00:09:45,990 --> 00:09:49,110 and then I wanna go from time T equals zero 197 00:09:49,110 --> 00:09:51,420 to time T F in steps of H. 198 00:09:51,420 --> 00:09:53,850 So again, this is similar to a N P R range. 199 00:09:53,850 --> 00:09:57,780 So this four loop is just doing as many time steps of H 200 00:09:57,780 --> 00:10:00,660 that it needs until it reaches time T F. 201 00:10:00,660 --> 00:10:03,870 At every time step, it's gonna look at the current state. 202 00:10:03,870 --> 00:10:04,703 So this is just asking, 203 00:10:04,703 --> 00:10:07,686 "What is the current state of the system?" 204 00:10:07,686 --> 00:10:08,570 It's gonna calculate 205 00:10:08,570 --> 00:10:11,400 in the same way I was doing in the previous code, 206 00:10:11,400 --> 00:10:14,700 calculate the derivatives at this point, 207 00:10:14,700 --> 00:10:17,673 and then simply add to the current state Y, 208 00:10:19,770 --> 00:10:22,050 add this first term of our Taylor series. 209 00:10:22,050 --> 00:10:27,050 So follow this derivative D X times D T for time H, right? 210 00:10:27,360 --> 00:10:31,500 So D X times D T, or D X D T, I should say, 211 00:10:31,500 --> 00:10:33,330 is simply the speed of our system. 212 00:10:33,330 --> 00:10:35,790 So if you're going somewhere at 10 miles per hour, 213 00:10:35,790 --> 00:10:38,370 and you're doing it for a time step of 10 hours, 214 00:10:38,370 --> 00:10:40,560 you get 10 times 10 miles per hour, 215 00:10:40,560 --> 00:10:42,750 a hundred miles done in 10 hours. 216 00:10:42,750 --> 00:10:47,750 So this is the amount of change done over that time steps, 217 00:10:47,850 --> 00:10:49,800 and we simply keep track of all time 218 00:10:49,800 --> 00:10:51,780 with this little trick here. 219 00:10:51,780 --> 00:10:55,113 Going back to that compare integrators code, 220 00:10:55,973 --> 00:10:58,920 you'll see another good programming practice here, which is 221 00:10:58,920 --> 00:11:01,715 always include a sample call. 222 00:11:01,715 --> 00:11:02,970 So when I first opened this file, 223 00:11:02,970 --> 00:11:04,170 there were some modifications 224 00:11:04,170 --> 00:11:06,390 that were done in discussion groups. 225 00:11:06,390 --> 00:11:10,200 Sometimes you'll modify things and forget about them. 226 00:11:10,200 --> 00:11:12,060 It might take a year until you go back to your code. 227 00:11:12,060 --> 00:11:13,830 It's always useful to be able 228 00:11:13,830 --> 00:11:16,773 and copy a function call to get it running. 229 00:11:18,810 --> 00:11:20,640 So here, if I run this command, 230 00:11:20,640 --> 00:11:25,573 we expect to run our model with alpha equals 0.25. 231 00:11:26,550 --> 00:11:28,920 I have my initial values here set randomly 232 00:11:28,920 --> 00:11:30,120 or set to 0.2, 0.3, 0.4, 233 00:11:31,080 --> 00:11:34,020 but those were not chosen in any meaningful way. 234 00:11:34,020 --> 00:11:39,020 I expect two figures to appear. 235 00:11:40,350 --> 00:11:43,885 Figure one will be a plot comparing the time series 236 00:11:43,885 --> 00:11:48,030 from the O D E four five, our really advanced integrator, 237 00:11:48,030 --> 00:11:52,350 to, in black, the time series obtained 238 00:11:52,350 --> 00:11:54,461 with our Euler integrator. 239 00:11:54,461 --> 00:11:57,720 So red will be the advanced integrator and black the Euler. 240 00:11:57,720 --> 00:12:00,630 And in the second plot, I expect a three D plot 241 00:12:00,630 --> 00:12:03,280 of the trajectory obtained with the Euler integrator. 242 00:12:09,420 --> 00:12:11,310 So here's our first figure. 243 00:12:11,310 --> 00:12:15,480 So again, in red you have the really advanced integrator 244 00:12:15,480 --> 00:12:19,413 and in black, which we obtained with Euler, 245 00:12:20,610 --> 00:12:23,790 so notice that they don't match perfectly, right? 246 00:12:23,790 --> 00:12:28,320 And we know that Euler does better when H is very small. 247 00:12:28,320 --> 00:12:31,530 So if H was five times as big, for example, 248 00:12:31,530 --> 00:12:33,960 the difference would be even bigger. 249 00:12:33,960 --> 00:12:37,143 And if H was instead five times smaller, 250 00:12:39,090 --> 00:12:40,650 it takes a little longer to run 251 00:12:40,650 --> 00:12:43,563 but now my two curves are almost indistinguishable. 252 00:12:45,090 --> 00:12:47,073 So going back to that original one, 253 00:12:49,290 --> 00:12:54,290 and then we get the three D plot of that Euler curve. 254 00:12:54,990 --> 00:12:57,930 So we start in the bottom corner here, 255 00:12:57,930 --> 00:12:59,580 and as our system evolves, 256 00:12:59,580 --> 00:13:02,400 it eventually settles around a steady state, 257 00:13:02,400 --> 00:13:05,220 and we know this by looking at this curve. 258 00:13:05,220 --> 00:13:07,713 Now, I can play with different values of alpha, 259 00:13:10,950 --> 00:13:14,644 so if I make alpha bigger, we'll notice that 260 00:13:14,644 --> 00:13:18,480 an overall time, capital T, equal to 50 261 00:13:18,480 --> 00:13:21,120 is not enough to settle to a steady state value, 262 00:13:21,120 --> 00:13:22,530 and we see the same thing here. 263 00:13:22,530 --> 00:13:25,590 We see an interesting behavior, the start of a swirl here. 264 00:13:25,590 --> 00:13:30,540 So we can integrate for longer, to 200, 265 00:13:30,540 --> 00:13:32,841 to see what that actually looks like. 266 00:13:32,841 --> 00:13:35,130 And we see that we're spiraling toward a fixed point, 267 00:13:35,130 --> 00:13:36,543 eventually reaching it. 268 00:13:38,280 --> 00:13:41,403 And our integrator is still doing really well. 269 00:13:43,770 --> 00:13:46,710 But again, if our time steps were too big, 270 00:13:46,710 --> 00:13:50,370 we might get a little, or a big difference in this case, 271 00:13:50,370 --> 00:13:53,610 in time to convergence to that fixed point. 272 00:13:53,610 --> 00:13:56,526 So one important thing to keep in mind is that 273 00:13:56,526 --> 00:13:58,890 it's really important to run your integrator 274 00:13:58,890 --> 00:14:02,130 with maybe a larger time step at first 275 00:14:02,130 --> 00:14:04,892 just to make sure that things are running. 276 00:14:04,892 --> 00:14:06,600 But if you want to validate a trajectory, 277 00:14:06,600 --> 00:14:09,600 then you really have to decrease your time step 278 00:14:09,600 --> 00:14:12,300 and make sure that you get the same time series. 279 00:14:12,300 --> 00:14:14,370 And once your time series doesn't change much, 280 00:14:14,370 --> 00:14:16,356 with your times that you know 281 00:14:16,356 --> 00:14:18,506 you've converged onto something meaningful. 282 00:14:20,370 --> 00:14:22,863 Let's bring alpha up a little more, 283 00:14:27,120 --> 00:14:30,723 and now we get some really nice oscillations going. 284 00:14:32,880 --> 00:14:34,710 We can integrate further to see 285 00:14:34,710 --> 00:14:37,180 if it's just the beginning of a spiral 286 00:14:38,100 --> 00:14:39,750 like we had with the previous value. 287 00:14:39,750 --> 00:14:42,120 And in this case we see that we're really capturing 288 00:14:42,120 --> 00:14:45,150 sort of what we had in the numerical simulation 289 00:14:45,150 --> 00:14:46,770 of the Lotka-Volterra model. 290 00:14:46,770 --> 00:14:50,670 We have resilient, robust oscillation and a periodic system 291 00:14:50,670 --> 00:14:54,270 where we're now orbiting, our system is orbiting 292 00:14:54,270 --> 00:14:57,000 in this three D space, three different species. 293 00:14:57,000 --> 00:15:01,503 I'm only plotting species one on this top figure right now. 294 00:15:02,801 --> 00:15:07,801 I can increase my alpha value a little bit more. 295 00:15:08,190 --> 00:15:10,110 If I do so, what I notice is that 296 00:15:10,110 --> 00:15:15,110 I'm slightly changing the length and duration of my orbit. 297 00:15:15,120 --> 00:15:16,500 So in this case here, 298 00:15:16,500 --> 00:15:20,190 the period is now roughly twice as much, 299 00:15:20,190 --> 00:15:22,530 and we see that we have two different peaks 300 00:15:22,530 --> 00:15:24,270 and we see the same thing in the trajectories 301 00:15:24,270 --> 00:15:28,230 where we're now tracing a slightly more complicated orbit. 302 00:15:28,230 --> 00:15:29,550 So as we increase alpha, 303 00:15:29,550 --> 00:15:31,833 we're getting richer and richer behavior. 304 00:15:34,230 --> 00:15:36,900 If I increase it to 1.5, 305 00:15:36,900 --> 00:15:40,530 well, now, even though I'm still using a times step of 0.1, 306 00:15:40,530 --> 00:15:43,590 I start seeing a really big difference between the behavior 307 00:15:43,590 --> 00:15:46,890 predicted by my good integrator in red 308 00:15:46,890 --> 00:15:49,680 and by Euler in black. 309 00:15:49,680 --> 00:15:52,440 So let me decrease my time steps, 310 00:15:52,440 --> 00:15:53,940 because I know something's up. 311 00:15:56,250 --> 00:15:58,923 That's starting to look better, but not quite. 312 00:15:59,810 --> 00:16:00,643 Smaller time step. 313 00:16:02,105 --> 00:16:04,211 And remember, every time I decrease my time step, 314 00:16:04,211 --> 00:16:07,233 I have to do more integration steps so it runs slower. 315 00:16:07,233 --> 00:16:08,983 So that's why I'm not overdoing it. 316 00:16:11,730 --> 00:16:14,010 All right, I'm starting to do pretty well 317 00:16:14,010 --> 00:16:17,613 and I can go see what the future looks like for my system. 318 00:16:19,080 --> 00:16:21,300 This is gonna take three times as long 319 00:16:21,300 --> 00:16:25,443 because I'm now gonna integrate to 150 instead of 50. 320 00:16:30,150 --> 00:16:34,860 And once that's done running we'll be able to see 321 00:16:34,860 --> 00:16:39,480 whether or not our time step of 0.005 was really enough. 322 00:16:39,480 --> 00:16:42,600 Eventually we see that we're not predicting 323 00:16:42,600 --> 00:16:44,310 the state of the system accurately at all. 324 00:16:44,310 --> 00:16:46,239 We have a huge error. 325 00:16:46,239 --> 00:16:47,550 So even though we were doing great at first, 326 00:16:47,550 --> 00:16:50,790 eventually our two trajectories have sort of diverged. 327 00:16:50,790 --> 00:16:55,410 And we could keep playing with our time step 328 00:16:55,410 --> 00:16:59,400 to try and see if there is a time step small enough 329 00:16:59,400 --> 00:17:02,250 to get the same trajectories over a long period of time, 330 00:17:02,250 --> 00:17:05,400 but I can tell you that's not gonna be the case. 331 00:17:05,400 --> 00:17:08,580 The reason for that is that the system has now entered 332 00:17:08,580 --> 00:17:10,980 what we call a chaotic regime 333 00:17:10,980 --> 00:17:13,920 where it is sensitive to initial conditions. 334 00:17:13,920 --> 00:17:16,290 So there is a buildup of error 335 00:17:16,290 --> 00:17:18,060 in that first part of the integration 336 00:17:18,060 --> 00:17:21,082 and eventually the two trajectories, 337 00:17:21,082 --> 00:17:23,430 even though they're very similar, diverge. 338 00:17:23,430 --> 00:17:26,230 To visualize that, I'll introduce 339 00:17:28,110 --> 00:17:30,220 a second set of initial conditions 340 00:17:35,520 --> 00:17:37,650 in my code in the top left. 341 00:17:37,650 --> 00:17:40,170 So I'm gonna use random numbers 342 00:17:40,170 --> 00:17:43,410 to just randomly perturb my initial condition 343 00:17:43,410 --> 00:17:46,290 between zero and 0.001. 344 00:17:46,290 --> 00:17:50,460 So very small perturbation on my initial conditions. 345 00:17:50,460 --> 00:17:53,070 I'm gonna integrate from these conditions 346 00:17:53,070 --> 00:17:56,010 using the Euler algorithm. 347 00:17:56,010 --> 00:18:00,420 So I'm gonna have a series of trajectories, 348 00:18:00,420 --> 00:18:02,700 a trajectory used by Euler, 349 00:18:02,700 --> 00:18:04,780 which I'm now gonna call Y E two here 350 00:18:05,635 --> 00:18:06,893 for a second set of Euler. 351 00:18:07,920 --> 00:18:11,150 Instead of comparing the Euler integrator 352 00:18:14,010 --> 00:18:17,766 to the good integrator, the O D E four five approach, 353 00:18:17,766 --> 00:18:20,670 we're gonna compare our old set of initial condition in red 354 00:18:20,670 --> 00:18:25,129 and the new randomly perturbed set of initial condition 355 00:18:25,129 --> 00:18:25,962 in black. 356 00:18:26,865 --> 00:18:31,230 We're gonna keep our three D figure here for figure two, 357 00:18:31,230 --> 00:18:34,110 but I'm gonna introduce a third figure 358 00:18:34,110 --> 00:18:35,610 which is gonna look at the difference 359 00:18:35,610 --> 00:18:38,640 between the black and the red curve in the first plot. 360 00:18:41,430 --> 00:18:43,553 To edit some of this out, 361 00:18:43,553 --> 00:18:45,210 you've probably just had a weird time jump. 362 00:18:45,210 --> 00:18:48,120 I was getting an error because in my code I forgotten 363 00:18:48,120 --> 00:18:51,120 to transpose my second solution. 364 00:18:51,120 --> 00:18:54,180 I just need, because of my comment here, 365 00:18:54,180 --> 00:18:56,040 I wanted to make sure I was using 366 00:18:56,040 --> 00:18:57,870 vector columns for all of them. 367 00:18:57,870 --> 00:19:02,870 Okay, so we can run and now compare our solutions. 368 00:19:03,690 --> 00:19:05,463 You see one already on the screen. 369 00:19:06,540 --> 00:19:11,540 So now our two Euler solutions are very similar. 370 00:19:11,740 --> 00:19:14,883 Again, I'm gonna go a little further in time, 371 00:19:24,914 --> 00:19:28,710 and eventually they start diverging, right? 372 00:19:28,710 --> 00:19:31,260 So even though they start very closely 373 00:19:31,260 --> 00:19:34,110 and you might expect them to stay close forever, 374 00:19:34,110 --> 00:19:37,650 because there is this sensitivity to initial conditions 375 00:19:37,650 --> 00:19:39,183 in the chaotic system, 376 00:19:40,410 --> 00:19:45,410 we expect eventually the two time series to diverge. 377 00:19:45,780 --> 00:19:47,460 And this is what I'm looking at 378 00:19:47,460 --> 00:19:51,236 with my third Figure here 379 00:19:51,236 --> 00:19:55,080 is the distance between these two time series 380 00:19:55,080 --> 00:19:56,313 as time goes on. 381 00:19:57,780 --> 00:20:00,660 And I'm using a semi-log scale, 382 00:20:00,660 --> 00:20:02,970 so a log scale in the Y axis, 383 00:20:02,970 --> 00:20:06,810 and as you can see, it goes up roughly exponentially. 384 00:20:06,810 --> 00:20:09,210 I say roughly because you're in a limited space 385 00:20:09,210 --> 00:20:11,148 and they're moving around, 386 00:20:11,148 --> 00:20:13,023 so you do expect some up and down 387 00:20:13,023 --> 00:20:15,423 just due to the constrained nature of the space. 388 00:20:16,260 --> 00:20:20,730 But you have this general exponential increase in distance 389 00:20:20,730 --> 00:20:22,290 between the two time series. 390 00:20:22,290 --> 00:20:24,510 The exponent, this exponential rate 391 00:20:24,510 --> 00:20:27,210 of separation after initial condition 392 00:20:27,210 --> 00:20:29,370 is what's called the Lyapunov exponent 393 00:20:29,370 --> 00:20:33,750 and it is a measure of how chaotic the system might be, 394 00:20:33,750 --> 00:20:37,713 how sensitive to initial condition the system might be. 395 00:20:39,000 --> 00:20:43,293 So I can now run my system for much longer, 396 00:20:46,620 --> 00:20:49,560 say, up time 400 instead of 200. 397 00:20:49,560 --> 00:20:51,090 I expect the two time series 398 00:20:51,090 --> 00:20:53,910 to eventually diverge completely. 399 00:20:53,910 --> 00:20:55,500 I've used a larger time step 400 00:20:55,500 --> 00:20:59,523 to be a little more accurate in my integration. 401 00:21:09,300 --> 00:21:11,875 And I'm gonna go even further. 402 00:21:11,875 --> 00:21:12,977 Eventually, as you'll see, 403 00:21:12,977 --> 00:21:14,160 the separation between the two methods, 404 00:21:14,160 --> 00:21:16,830 we see this right now in figure three, saturates 405 00:21:16,830 --> 00:21:20,070 because there is just so much separation they can take 406 00:21:20,070 --> 00:21:23,070 because they're constrained to a final space 407 00:21:23,070 --> 00:21:25,593 just because of the way this model is specified. 408 00:21:28,650 --> 00:21:31,408 And as we integrate more and more, 409 00:21:31,408 --> 00:21:32,400 we keep exploring more of the space. 410 00:21:32,400 --> 00:21:36,420 If you're looking at this bottom figure here, my figure two, 411 00:21:36,420 --> 00:21:39,783 I keep exploring a larger fraction of space. 412 00:21:41,010 --> 00:21:46,010 So this nice shape here is what's called a strange attractor 413 00:21:46,770 --> 00:21:49,500 because no matter where I start in this system, 414 00:21:49,500 --> 00:21:52,720 I will end up on this surface or volume 415 00:21:54,270 --> 00:21:56,550 and once I'm on it, my chaotic system 416 00:21:56,550 --> 00:21:59,310 will keep exploring it, exploring every point. 417 00:21:59,310 --> 00:22:01,230 This has come up in discussion before. 418 00:22:01,230 --> 00:22:05,040 It's what we call ergodicity, and chaos is exactly this. 419 00:22:05,040 --> 00:22:07,469 It's a deterministic system. 420 00:22:07,469 --> 00:22:09,450 It's not about being unpredictable or random. 421 00:22:09,450 --> 00:22:11,100 It's a deterministic system 422 00:22:11,100 --> 00:22:13,560 that is sensitive to initial condition. 423 00:22:13,560 --> 00:22:18,180 It has this exponential separation of trajectories 424 00:22:18,180 --> 00:22:23,180 and those trajectories are aperiodic and ergodic, 425 00:22:23,250 --> 00:22:24,990 meaning they don't repeat themselves 426 00:22:24,990 --> 00:22:27,930 and they eventually explore all of space. 427 00:22:27,930 --> 00:22:30,750 So if I kept integrating for longer and longer, 428 00:22:30,750 --> 00:22:34,773 I would fill up this entire surface completely. 429 00:22:36,750 --> 00:22:39,990 So that's it for now, that's our first taste of chaos. 430 00:22:39,990 --> 00:22:42,740 The textbook has an entire chapter of it. 431 00:22:42,740 --> 00:22:44,790 There is a great book by Steven Strogatz 432 00:22:44,790 --> 00:22:47,940 on non-linear dynamics if you wanna study this topic more. 433 00:22:47,940 --> 00:22:51,745 We also have a course on it here at UVM. 434 00:22:51,745 --> 00:22:53,760 Otherwise, I'll see you in discussion 435 00:22:53,760 --> 00:22:56,324 where we can go a little deeply into 436 00:22:56,324 --> 00:22:58,620 the theory and implications of chaos. 437 00:22:58,620 --> 00:23:01,500 And I'll see you in the next video to wrap up 438 00:23:01,500 --> 00:23:03,810 our first module on dynamical systems. 439 00:23:03,810 --> 00:23:04,643 Thank you.