1 00:00:03,680 --> 00:00:09,570 all right so I'm gonna be talking about 2 00:00:07,500 --> 00:00:11,910 numerical methods I'm Erin kado 3 00:00:09,570 --> 00:00:15,840 I work at Blizzard Entertainment I do a 4 00:00:11,910 --> 00:00:17,099 lot of physics work there you've 5 00:00:15,840 --> 00:00:18,660 probably seen something you'll see a 6 00:00:17,099 --> 00:00:21,080 video later I'll give you a little hint 7 00:00:18,660 --> 00:00:25,980 about that 8 00:00:21,080 --> 00:00:29,359 so numerical methods in game physics 9 00:00:25,980 --> 00:00:32,010 were often dealing with really complex 10 00:00:29,359 --> 00:00:36,559 problems we have lots of moving objects 11 00:00:32,010 --> 00:00:38,790 friction contact constraints and these 12 00:00:36,559 --> 00:00:41,100 situations are very complex and we have 13 00:00:38,790 --> 00:00:45,329 to use numerical methods to solve these 14 00:00:41,100 --> 00:00:47,489 situations in particular we need to use 15 00:00:45,329 --> 00:00:49,710 numerical methods to get approximate 16 00:00:47,489 --> 00:00:52,020 solutions to all the differential 17 00:00:49,710 --> 00:00:55,920 equations that pop up as a result of 18 00:00:52,020 --> 00:00:59,010 that so there are many numerical methods 19 00:00:55,920 --> 00:01:01,109 available but it isn't always clear what 20 00:00:59,010 --> 00:01:04,020 is the best numerical method to use what 21 00:01:01,109 --> 00:01:09,000 is the right accuracy performance 22 00:01:04,020 --> 00:01:11,070 trade-off and also once you've chosen 23 00:01:09,000 --> 00:01:13,080 method you know what are the stability 24 00:01:11,070 --> 00:01:15,000 limits of that method you know is this 25 00:01:13,080 --> 00:01:18,360 method gonna blow up you know you don't 26 00:01:15,000 --> 00:01:20,220 want to be about to ship your game and 27 00:01:18,360 --> 00:01:23,189 all some things are blowing up and you 28 00:01:20,220 --> 00:01:25,140 didn't know why so it's it's good to 29 00:01:23,189 --> 00:01:27,270 have a good understanding of numerical 30 00:01:25,140 --> 00:01:30,000 methods and that's kind of the goal of 31 00:01:27,270 --> 00:01:34,020 this presentation is to give you a 32 00:01:30,000 --> 00:01:36,299 better kind of background and way to 33 00:01:34,020 --> 00:01:40,290 think about these things to help you 34 00:01:36,299 --> 00:01:42,329 make better games all right I'm gonna 35 00:01:40,290 --> 00:01:47,250 start out first with a little bit of 36 00:01:42,329 --> 00:01:49,610 review of differential calculus so 37 00:01:47,250 --> 00:01:53,070 differential calculus starts with the 38 00:01:49,610 --> 00:01:55,229 derivative so here's the the definition 39 00:01:53,070 --> 00:01:59,759 of the derivative so say I have some 40 00:01:55,229 --> 00:02:02,430 function f that is a function of the x 41 00:01:59,759 --> 00:02:06,240 coordinate and so that's this curve here 42 00:02:02,430 --> 00:02:09,179 and I want to get the slope at any point 43 00:02:06,240 --> 00:02:11,009 along that curve well you can imagine I 44 00:02:09,179 --> 00:02:12,780 can approximate that by drawing the 45 00:02:11,009 --> 00:02:15,110 secant through there sample of two 46 00:02:12,780 --> 00:02:18,920 points and then I have a step with 47 00:02:15,110 --> 00:02:21,010 eh and I draw a line through there and 48 00:02:18,920 --> 00:02:25,040 that's kind of the approximate slope 49 00:02:21,010 --> 00:02:26,960 well if we take that step size H and 50 00:02:25,040 --> 00:02:29,830 make that smaller and smaller and the 51 00:02:26,960 --> 00:02:33,050 limit as that goes to zero our 52 00:02:29,830 --> 00:02:38,030 approximation of the slope at that point 53 00:02:33,050 --> 00:02:40,970 X is improved and so that's kind of the 54 00:02:38,030 --> 00:02:43,340 definition of the derivative and lots of 55 00:02:40,970 --> 00:02:45,260 formulas for computing derivatives and 56 00:02:43,340 --> 00:02:48,320 giving you the exact slope at any point 57 00:02:45,260 --> 00:02:50,720 along a curve if those have been 58 00:02:48,320 --> 00:02:55,250 developed and so I'll go through a few 59 00:02:50,720 --> 00:02:58,730 examples of that so the derivative of a 60 00:02:55,250 --> 00:03:01,700 line is well that's that's the slope of 61 00:02:58,730 --> 00:03:04,190 the line so the slope of the line is the 62 00:03:01,700 --> 00:03:07,100 same everywhere along the line so if I 63 00:03:04,190 --> 00:03:10,160 have my my function for the line you 64 00:03:07,100 --> 00:03:14,750 know y equals MX plus B the derivative 65 00:03:10,160 --> 00:03:25,459 of that is just M the slope or you know 66 00:03:14,750 --> 00:03:28,940 the rise over the over the run so what 67 00:03:25,459 --> 00:03:31,130 if we have a like a curve like a 68 00:03:28,940 --> 00:03:34,970 parabola well there that the slope is 69 00:03:31,130 --> 00:03:37,220 not the same at every point so so we 70 00:03:34,970 --> 00:03:39,170 need to use you know the ideas from 71 00:03:37,220 --> 00:03:42,340 differential calculus to compute that 72 00:03:39,170 --> 00:03:45,800 that slope so if I have this parabola 73 00:03:42,340 --> 00:03:49,459 that's y equals ax squared where a is 74 00:03:45,800 --> 00:03:52,580 just an arbitrary constant I take the 75 00:03:49,459 --> 00:03:56,750 derivative of that and that's to ax so 76 00:03:52,580 --> 00:04:04,190 the the slope varies linearly along the 77 00:03:56,750 --> 00:04:06,080 parabola when it gets when we get two 78 00:04:04,190 --> 00:04:08,360 differential equations the the 79 00:04:06,080 --> 00:04:12,410 exponential function is very interesting 80 00:04:08,360 --> 00:04:14,780 so the exponential function has this 81 00:04:12,410 --> 00:04:17,419 really interesting property that if I 82 00:04:14,780 --> 00:04:20,630 take its derivative I get back another 83 00:04:17,419 --> 00:04:23,780 accident exponential so here I have y 84 00:04:20,630 --> 00:04:26,810 equals e a X where a is yet again an 85 00:04:23,780 --> 00:04:27,980 arbitrary constant and then I compute 86 00:04:26,810 --> 00:04:33,200 the derivative of that 87 00:04:27,980 --> 00:04:38,660 get a e to the ax so friend of mine 88 00:04:33,200 --> 00:04:41,060 found this cute little cartoon so a wild 89 00:04:38,660 --> 00:04:43,430 exponential function appeared you use 90 00:04:41,060 --> 00:04:47,950 differentiate it is not very effective 91 00:04:43,430 --> 00:04:47,950 because you get back another exponential 92 00:04:49,510 --> 00:04:56,330 when we look at you know simulation the 93 00:04:54,080 --> 00:04:59,780 sine wave is another function that's 94 00:04:56,330 --> 00:05:01,880 very important if I compute the 95 00:04:59,780 --> 00:05:04,700 derivative of a sine wave I get back 96 00:05:01,880 --> 00:05:07,130 basically another sine wave shifted by 97 00:05:04,700 --> 00:05:09,620 90 degrees so when I compute the 98 00:05:07,130 --> 00:05:14,200 derivative of sine I get cosine which is 99 00:05:09,620 --> 00:05:14,200 basically signed as shifted 90 degrees 100 00:05:15,880 --> 00:05:21,860 all right so in videogame physics were 101 00:05:18,950 --> 00:05:24,620 often dealing with functions of time so 102 00:05:21,860 --> 00:05:27,680 I'm going to move now to having time is 103 00:05:24,620 --> 00:05:29,660 my independent variable and then I'm 104 00:05:27,680 --> 00:05:31,850 going to use X a lot as the dependent 105 00:05:29,660 --> 00:05:35,420 variable so X you can think of it as the 106 00:05:31,850 --> 00:05:37,070 position and then if I compute the time 107 00:05:35,420 --> 00:05:39,830 derivative of the position I get the 108 00:05:37,070 --> 00:05:40,490 velocity and then the time derivative of 109 00:05:39,830 --> 00:05:43,490 velocity 110 00:05:40,490 --> 00:05:46,370 I get the acceleration well since we're 111 00:05:43,490 --> 00:05:49,760 doing this a lot writing all this dx/dt 112 00:05:46,370 --> 00:05:52,730 stuff gets a little busy so it's a 113 00:05:49,760 --> 00:05:57,320 standard notation to use the dot to 114 00:05:52,730 --> 00:06:01,250 represent the derivative so X dot is DX 115 00:05:57,320 --> 00:06:04,070 DT and I can have for the acceleration X 116 00:06:01,250 --> 00:06:08,000 double dot so I'm gonna be using that 117 00:06:04,070 --> 00:06:10,000 notation going forward all right a 118 00:06:08,000 --> 00:06:13,100 little bit more review let's talk about 119 00:06:10,000 --> 00:06:16,250 differential equations so the main 120 00:06:13,100 --> 00:06:18,860 premise here is okay what if I had the 121 00:06:16,250 --> 00:06:21,680 slope and I want to get back to the 122 00:06:18,860 --> 00:06:23,990 original function why is that 123 00:06:21,680 --> 00:06:26,780 interesting because you know Newton's 124 00:06:23,990 --> 00:06:28,610 law as you'll see that kind of that 125 00:06:26,780 --> 00:06:32,150 gives us the slope and when I get back 126 00:06:28,610 --> 00:06:34,160 to the original function so in 127 00:06:32,150 --> 00:06:36,170 particular we're going to be looking at 128 00:06:34,160 --> 00:06:38,450 what's called a first-order differential 129 00:06:36,170 --> 00:06:40,169 equation so here's the standard form for 130 00:06:38,450 --> 00:06:43,770 that 131 00:06:40,169 --> 00:06:48,990 to X dot so that's DX DT is equal to 132 00:06:43,770 --> 00:06:51,150 some function of time and X itself so 133 00:06:48,990 --> 00:06:54,060 what's this function well there's all 134 00:06:51,150 --> 00:06:56,610 kinds of functions that we might come 135 00:06:54,060 --> 00:07:00,870 upon but we can kind of talk about this 136 00:06:56,610 --> 00:07:03,150 in a generic way but important point 137 00:07:00,870 --> 00:07:04,800 here is that it is a first order 138 00:07:03,150 --> 00:07:14,219 differential equation because it's just 139 00:07:04,800 --> 00:07:16,169 involving X not X double dot so this is 140 00:07:14,219 --> 00:07:20,430 kind of what we want to do we want to 141 00:07:16,169 --> 00:07:25,229 have this formula for the derivative of 142 00:07:20,430 --> 00:07:27,029 X and then get back X so that's that is 143 00:07:25,229 --> 00:07:33,629 what it means to solve a differential 144 00:07:27,029 --> 00:07:36,210 equation if that function is just a 145 00:07:33,629 --> 00:07:38,810 function of time it's kind of 146 00:07:36,210 --> 00:07:42,270 straightforward it's just like reverse 147 00:07:38,810 --> 00:07:46,020 differentiation or integration if that 148 00:07:42,270 --> 00:07:48,180 function depends on X then things get a 149 00:07:46,020 --> 00:07:52,199 little bit more complex I'll show some 150 00:07:48,180 --> 00:07:54,330 examples all right so I want to create a 151 00:07:52,199 --> 00:07:55,860 bunch of examples and I was thinking 152 00:07:54,330 --> 00:07:58,529 well what's the simplest differential 153 00:07:55,860 --> 00:08:01,759 equation I could have well how about X 154 00:07:58,529 --> 00:08:04,199 dot equals 0 so that means that the 155 00:08:01,759 --> 00:08:12,360 derivative is zero everywhere that means 156 00:08:04,199 --> 00:08:16,740 the slope is constant and it's zero and 157 00:08:12,360 --> 00:08:18,749 so then the solution is just that the X 158 00:08:16,740 --> 00:08:23,419 is just constant it's equal to the 159 00:08:18,749 --> 00:08:28,050 initial value of X x0 pretty exciting 160 00:08:23,419 --> 00:08:31,050 all right so let's look at something a 161 00:08:28,050 --> 00:08:35,279 little bit more complex so X dot is now 162 00:08:31,050 --> 00:08:37,289 equal to a nonzero constant a so the 163 00:08:35,279 --> 00:08:40,199 solution to that is well that's a line 164 00:08:37,289 --> 00:08:42,630 because before we saw well if you have a 165 00:08:40,199 --> 00:08:46,920 constant slope and that unless that's 166 00:08:42,630 --> 00:08:49,709 what a line has so x equals 80 plus x 0 167 00:08:46,920 --> 00:08:53,600 where X 0 again is the initial value of 168 00:08:49,709 --> 00:08:53,600 X at T equals 0 169 00:08:54,470 --> 00:08:59,430 okay but as I said things see a little 170 00:08:57,060 --> 00:09:03,180 bit more complex if X is on the 171 00:08:59,430 --> 00:09:07,019 right-hand side so here X dot equals a X 172 00:09:03,180 --> 00:09:11,940 so that means that the slope of X 173 00:09:07,019 --> 00:09:13,500 depends on the value of X well it turns 174 00:09:11,940 --> 00:09:17,610 out the solution to this is the 175 00:09:13,500 --> 00:09:21,329 exponential so x equals x 0 e to the a T 176 00:09:17,610 --> 00:09:21,899 so X 0 is the initial value that makes 177 00:09:21,329 --> 00:09:25,589 sense 178 00:09:21,899 --> 00:09:31,160 if I put in T equals 0 exponential e to 179 00:09:25,589 --> 00:09:33,779 the 0 is 1 so X at T equals 0 is X 0 and 180 00:09:31,160 --> 00:09:37,440 this is where for me it gets a little 181 00:09:33,779 --> 00:09:39,509 bit weird that you don't really derive 182 00:09:37,440 --> 00:09:42,420 the solution you kind of guess the 183 00:09:39,509 --> 00:09:46,800 solution and prove that it's true this 184 00:09:42,420 --> 00:09:48,509 is how calculus works so you can see 185 00:09:46,800 --> 00:09:53,459 down here on the the bottom is the proof 186 00:09:48,509 --> 00:09:55,079 I'll substitute in my solution X 0 into 187 00:09:53,459 --> 00:09:57,569 that top equation take its derivative 188 00:09:55,079 --> 00:10:01,709 and then kind of regroup things and say 189 00:09:57,569 --> 00:10:06,510 oh that is that is equal to a times X so 190 00:10:01,709 --> 00:10:08,670 that's the proof all right of course in 191 00:10:06,510 --> 00:10:10,050 game physics were we're interested in 192 00:10:08,670 --> 00:10:15,269 Newton's second law so mass times 193 00:10:10,050 --> 00:10:18,180 acceleration equals force and in 194 00:10:15,269 --> 00:10:20,160 particular if we look in one dimension x 195 00:10:18,180 --> 00:10:23,339 is just positioned along one the 196 00:10:20,160 --> 00:10:26,490 coordinate axis and X double dot equals 197 00:10:23,339 --> 00:10:28,410 force so that's a scalar equation and 198 00:10:26,490 --> 00:10:29,639 this is called a second order 199 00:10:28,410 --> 00:10:35,970 differential equation because it 200 00:10:29,639 --> 00:10:37,860 involves the second derivative of X all 201 00:10:35,970 --> 00:10:43,010 right so let's look at some examples of 202 00:10:37,860 --> 00:10:48,660 that so gravity this is a constant force 203 00:10:43,010 --> 00:10:54,089 so the force is negative mg and pulling 204 00:10:48,660 --> 00:10:58,199 this point mass down alright let's let's 205 00:10:54,089 --> 00:10:59,670 try to solve that well since the right 206 00:10:58,199 --> 00:11:01,920 hand side is constant I can just 207 00:10:59,670 --> 00:11:04,240 basically get the solution through 208 00:11:01,920 --> 00:11:07,339 reverse differentiation 209 00:11:04,240 --> 00:11:09,440 so I divide through by the mass so that 210 00:11:07,339 --> 00:11:11,600 shows you that the acceleration doesn't 211 00:11:09,440 --> 00:11:15,800 depend on the mass so y double dot 212 00:11:11,600 --> 00:11:19,430 equals negative mg and then I integrate 213 00:11:15,800 --> 00:11:23,560 that once and I get y dot equals v-0 214 00:11:19,430 --> 00:11:27,230 minus GT V zero is the initial velocity 215 00:11:23,560 --> 00:11:31,870 and then I integrate that one more time 216 00:11:27,230 --> 00:11:36,800 and I get this quadratic function for 217 00:11:31,870 --> 00:11:38,839 the vertical position Y and that 218 00:11:36,800 --> 00:11:41,029 involves the the initial position there 219 00:11:38,839 --> 00:11:43,550 as well and so that's how you get 220 00:11:41,029 --> 00:11:48,160 parabolic motion because it's this 221 00:11:43,550 --> 00:11:48,160 parabolic quadratic function 222 00:11:50,769 --> 00:11:55,700 alright so things get a little bit more 223 00:11:53,089 --> 00:11:57,529 interesting with Newton's law when we 224 00:11:55,700 --> 00:12:00,320 start to involve forces that depend on 225 00:11:57,529 --> 00:12:01,790 position and I'm gonna be using this 226 00:12:00,320 --> 00:12:03,529 example throughout the rest of the 227 00:12:01,790 --> 00:12:05,930 presentation this mass spring system 228 00:12:03,529 --> 00:12:09,470 so I'll gonna go over it and in a bit of 229 00:12:05,930 --> 00:12:11,089 detail so I have this point mass m and 230 00:12:09,470 --> 00:12:14,720 it can only move along the horizontal 231 00:12:11,089 --> 00:12:18,140 axis and it's connected to ground by 232 00:12:14,720 --> 00:12:20,329 this spring K and the way it works is 233 00:12:18,140 --> 00:12:23,180 well the more I stretch that spring the 234 00:12:20,329 --> 00:12:26,810 harder the force pulls back so the force 235 00:12:23,180 --> 00:12:29,120 is actually negative KX convention if 236 00:12:26,810 --> 00:12:32,420 you look on Wikipedia you're gonna see 237 00:12:29,120 --> 00:12:35,019 it written this form so the the force 238 00:12:32,420 --> 00:12:39,519 has been move over to the left hand side 239 00:12:35,019 --> 00:12:39,519 so there's your KX there 240 00:12:43,410 --> 00:12:54,629 all right so here's how you might expect 241 00:12:45,509 --> 00:13:02,239 that system to move so you would expect 242 00:12:54,629 --> 00:13:05,329 some nice oscillation all right so let's 243 00:13:02,239 --> 00:13:08,009 try to solve this thing so first thing 244 00:13:05,329 --> 00:13:10,859 I'm going to do is reduce the number of 245 00:13:08,009 --> 00:13:14,429 parameters so if I divide through by the 246 00:13:10,859 --> 00:13:15,629 mass then I can reduce down to one 247 00:13:14,429 --> 00:13:17,989 parameter and I'm going to use this 248 00:13:15,629 --> 00:13:19,169 parameter Omega which is the angular 249 00:13:17,989 --> 00:13:22,139 frequency 250 00:13:19,169 --> 00:13:25,199 so that's has units of radians per 251 00:13:22,139 --> 00:13:30,089 second and that is defined as the square 252 00:13:25,199 --> 00:13:33,299 root of K over m so notice that as the 253 00:13:30,089 --> 00:13:36,359 spring stiffness K increases the angular 254 00:13:33,299 --> 00:13:38,939 frequency also increases but as the mass 255 00:13:36,359 --> 00:13:47,449 increases the angular frequency goes 256 00:13:38,939 --> 00:13:50,729 down so here is the solution to that 257 00:13:47,449 --> 00:13:54,569 it's a combination of a sine and cosine 258 00:13:50,729 --> 00:13:58,350 and we have the initial position and 259 00:13:54,569 --> 00:14:01,229 velocity in there I noticed that the the 260 00:13:58,350 --> 00:14:04,199 angular frequency shows up inside the 261 00:14:01,229 --> 00:14:05,939 the sine and cosine functions so that 262 00:14:04,199 --> 00:14:12,319 kind of scales the number of 263 00:14:05,939 --> 00:14:15,659 oscillations per second that you get so 264 00:14:12,319 --> 00:14:17,669 here's a plot of the solution for the 265 00:14:15,659 --> 00:14:20,189 position and for the velocity versus 266 00:14:17,669 --> 00:14:24,529 time for some specific initial 267 00:14:20,189 --> 00:14:26,970 conditions and angular frequency of one 268 00:14:24,529 --> 00:14:32,339 so you can see it's just you know it's a 269 00:14:26,970 --> 00:14:34,829 cosine and a sine so what's interesting 270 00:14:32,339 --> 00:14:39,239 is like well if I'm going to be solving 271 00:14:34,829 --> 00:14:42,509 lots of these things maybe I could have 272 00:14:39,239 --> 00:14:46,169 another way to look at this data and I 273 00:14:42,509 --> 00:14:49,229 it an interesting way to look at the 274 00:14:46,169 --> 00:14:53,209 data is actually plotting the velocity 275 00:14:49,229 --> 00:14:57,040 versus the position and in this case 276 00:14:53,209 --> 00:15:00,740 time becomes a parameter on that curve 277 00:14:57,040 --> 00:15:03,400 so in this case the solution starts at 278 00:15:00,740 --> 00:15:10,070 that green circle and then goes around 279 00:15:03,400 --> 00:15:13,400 in a clockwise fashion and the angular 280 00:15:10,070 --> 00:15:17,200 frequency determines how quickly it goes 281 00:15:13,400 --> 00:15:17,200 around that circle as time progresses 282 00:15:18,760 --> 00:15:24,980 all right so I want to make this a 283 00:15:22,790 --> 00:15:27,980 little bit clearer so I I created an 284 00:15:24,980 --> 00:15:31,730 animation here and so what you're gonna 285 00:15:27,980 --> 00:15:34,490 see is on the upper left there that's 286 00:15:31,730 --> 00:15:37,850 going to be position versus time the 287 00:15:34,490 --> 00:15:40,250 bottom left is velocity versus time and 288 00:15:37,850 --> 00:15:43,970 then over on the right that's the phase 289 00:15:40,250 --> 00:15:51,980 portrait which is velocity versus 290 00:15:43,970 --> 00:15:54,280 position I'm sorry it's going it's going 291 00:15:51,980 --> 00:15:54,280 counterclockwise 292 00:15:59,780 --> 00:16:05,000 all right we're gonna be seeing more of 293 00:16:01,430 --> 00:16:08,390 these all right so like I said in the 294 00:16:05,000 --> 00:16:11,500 beginning and game physics we're dealing 295 00:16:08,390 --> 00:16:13,820 with lots of moving objects lots of 296 00:16:11,500 --> 00:16:16,610 complicated effects like friction and 297 00:16:13,820 --> 00:16:22,310 constraints and we can't come up with 298 00:16:16,610 --> 00:16:27,400 these nice pretty exact solutions so we 299 00:16:22,310 --> 00:16:29,540 have to turn to numerical methods so 300 00:16:27,400 --> 00:16:31,160 right off the bat when you when you 301 00:16:29,540 --> 00:16:34,940 switch to numerical methods there's a 302 00:16:31,160 --> 00:16:37,010 few things that you have to to deal with 303 00:16:34,940 --> 00:16:38,450 first is you're gonna be sacrificing 304 00:16:37,010 --> 00:16:39,980 accuracy you're not going to get the 305 00:16:38,450 --> 00:16:42,230 exact sine wave of more 306 00:16:39,980 --> 00:16:48,530 you may not get the exact parabolic 307 00:16:42,230 --> 00:16:51,320 motion anymore for video games this is 308 00:16:48,530 --> 00:16:55,850 almost always okay it's okay to 309 00:16:51,320 --> 00:16:59,240 sacrifice accuracy for video games we 310 00:16:55,850 --> 00:17:01,010 can sacrifice a lot of accuracy but a 311 00:16:59,240 --> 00:17:03,170 nice thing that you pick up is with 312 00:17:01,010 --> 00:17:04,760 numerical methods that they work very 313 00:17:03,170 --> 00:17:07,490 well with time stepping so if you have a 314 00:17:04,760 --> 00:17:10,459 game loop that is you know moving the 315 00:17:07,490 --> 00:17:12,740 game forward by the time step delta-t 316 00:17:10,459 --> 00:17:14,980 every frame well that's how are the 317 00:17:12,740 --> 00:17:17,689 numerical methods work they just like to 318 00:17:14,980 --> 00:17:20,449 move from from one solution to the next 319 00:17:17,689 --> 00:17:22,459 solution by a small amount of time so 320 00:17:20,449 --> 00:17:24,470 that works out well another thing that's 321 00:17:22,459 --> 00:17:28,459 nice about numerical methods is they 322 00:17:24,470 --> 00:17:31,970 they let you treat everything in a kind 323 00:17:28,459 --> 00:17:33,320 of generic way so if I have a whole 324 00:17:31,970 --> 00:17:35,510 bunch of differential equations I can 325 00:17:33,320 --> 00:17:37,880 just slam them into an array and then 326 00:17:35,510 --> 00:17:40,640 feed that off to my numerical method and 327 00:17:37,880 --> 00:17:43,670 then it fits back a result so it can 328 00:17:40,640 --> 00:17:51,010 kind of view your system that's kind of 329 00:17:43,670 --> 00:17:53,300 a black box but it does require your 330 00:17:51,010 --> 00:17:57,350 differential equations to be set up in a 331 00:17:53,300 --> 00:17:59,810 specific way so numerical methods like 332 00:17:57,350 --> 00:18:02,900 the differential equations to be set up 333 00:17:59,810 --> 00:18:04,640 in a first-order form so you can have 334 00:18:02,900 --> 00:18:06,230 just an array of first-order 335 00:18:04,640 --> 00:18:08,330 differential equations they can be 336 00:18:06,230 --> 00:18:10,550 mutually dependent they can be super 337 00:18:08,330 --> 00:18:11,850 complex in terms of the the forces or 338 00:18:10,550 --> 00:18:15,770 whatever there 339 00:18:11,850 --> 00:18:18,360 but it should be in this format so 340 00:18:15,770 --> 00:18:20,960 that's usually not a problem at all we 341 00:18:18,360 --> 00:18:23,730 can't we can assemble this and this is 342 00:18:20,960 --> 00:18:28,500 you know computers are good at this sort 343 00:18:23,730 --> 00:18:31,260 of thing alright but there is a problem 344 00:18:28,500 --> 00:18:34,470 so we're usually dealing with Newton's 345 00:18:31,260 --> 00:18:37,560 law so mass times X double dot equals 346 00:18:34,470 --> 00:18:40,110 have that's at least for one dimension 347 00:18:37,560 --> 00:18:43,230 you could view that maybe as X is a 348 00:18:40,110 --> 00:18:46,370 three dimensional quantity and forces 349 00:18:43,230 --> 00:18:49,320 three dimensional that would work too 350 00:18:46,370 --> 00:18:51,000 but the problem is that's a second-order 351 00:18:49,320 --> 00:18:52,950 differential equation and we want to get 352 00:18:51,000 --> 00:18:54,420 that into first-order form so that we 353 00:18:52,950 --> 00:18:59,310 can and hand it off to a numerical 354 00:18:54,420 --> 00:19:01,380 method so it's actually not that bad you 355 00:18:59,310 --> 00:19:03,090 can you can do that and and and the way 356 00:19:01,380 --> 00:19:05,250 you do it is by introducing another 357 00:19:03,090 --> 00:19:06,870 state so instead of just having X my 358 00:19:05,250 --> 00:19:09,210 position is my state I'm also gonna have 359 00:19:06,870 --> 00:19:14,820 velocity as my state and so then I can 360 00:19:09,210 --> 00:19:16,110 say okay the V dot the velocity the 361 00:19:14,820 --> 00:19:19,050 derivative of the velocity which is 362 00:19:16,110 --> 00:19:21,780 acceleration that's equal to F divided 363 00:19:19,050 --> 00:19:24,600 by M and then I can introduce this kind 364 00:19:21,780 --> 00:19:27,990 of trivial differential equation X dot 365 00:19:24,600 --> 00:19:29,580 equals V it is trivial but is it 366 00:19:27,990 --> 00:19:31,970 actually a valid differential equation 367 00:19:29,580 --> 00:19:34,110 and then I can just stack those two 368 00:19:31,970 --> 00:19:36,870 first order differential equations 369 00:19:34,110 --> 00:19:39,360 together and so what I am able to do is 370 00:19:36,870 --> 00:19:41,160 convert that one second order 371 00:19:39,360 --> 00:19:47,040 differential equation into two first 372 00:19:41,160 --> 00:19:50,730 order differential equations so it's 373 00:19:47,040 --> 00:19:54,240 pretty simple to make this happen all 374 00:19:50,730 --> 00:19:57,600 right so when we're getting into 375 00:19:54,240 --> 00:20:02,180 numerical methods a good place to start 376 00:19:57,600 --> 00:20:05,220 is looking at the finite difference so 377 00:20:02,180 --> 00:20:10,220 what I want to do here is show you a way 378 00:20:05,220 --> 00:20:13,320 to approximate the derivative with a 379 00:20:10,220 --> 00:20:15,720 basically a secant like we saw before so 380 00:20:13,320 --> 00:20:20,100 this is pretty much going back to that 381 00:20:15,720 --> 00:20:24,090 first slide about calculus so I can 382 00:20:20,100 --> 00:20:27,780 approximate the derivative by 383 00:20:24,090 --> 00:20:31,440 say I have my current state x1 and then 384 00:20:27,780 --> 00:20:32,940 I take a future state x2 compute the 385 00:20:31,440 --> 00:20:35,730 difference and divide that by the time 386 00:20:32,940 --> 00:20:39,750 step and that's an approximation of the 387 00:20:35,730 --> 00:20:41,760 the slope at X 1 now you can see there 388 00:20:39,750 --> 00:20:44,340 I've drawn with the dashed line the the 389 00:20:41,760 --> 00:20:46,800 actual slope at X 1 and then the green 390 00:20:44,340 --> 00:20:49,800 line that's my approximation so you can 391 00:20:46,800 --> 00:20:52,500 see it's kind of crude it gets better as 392 00:20:49,800 --> 00:20:56,810 the time step gets smaller obviously and 393 00:20:52,500 --> 00:20:56,810 this is called the forward difference 394 00:20:57,560 --> 00:21:03,210 you can also have the backwards 395 00:20:59,880 --> 00:21:07,500 difference so say I'm gonna store my old 396 00:21:03,210 --> 00:21:10,020 state x0 and then I'll take my current 397 00:21:07,500 --> 00:21:11,730 state - my old state and divide that by 398 00:21:10,020 --> 00:21:14,670 the time step and that's called the 399 00:21:11,730 --> 00:21:17,340 backward difference you're gonna see 400 00:21:14,670 --> 00:21:23,250 later how that how these impact our 401 00:21:17,340 --> 00:21:24,480 solution ok but let's let's go ahead and 402 00:21:23,250 --> 00:21:27,630 apply that forward difference 403 00:21:24,480 --> 00:21:29,910 approximation so what I'm doing here is 404 00:21:27,630 --> 00:21:32,670 I have my generic first order 405 00:21:29,910 --> 00:21:35,670 differential equation X dot equals F of 406 00:21:32,670 --> 00:21:39,090 T and X and I'm just gonna get rid of X 407 00:21:35,670 --> 00:21:42,720 dot by by inserting my approximation for 408 00:21:39,090 --> 00:21:46,590 X dot using the forward difference x2 409 00:21:42,720 --> 00:21:48,930 minus x1 over H and what that does is 410 00:21:46,590 --> 00:21:52,410 that turns that differential equation 411 00:21:48,930 --> 00:21:54,660 into an algebraic equation so down on 412 00:21:52,410 --> 00:21:58,890 the bottom there I just rearranged the 413 00:21:54,660 --> 00:22:02,540 terms so now I can compute X to my new 414 00:21:58,890 --> 00:22:05,460 state in terms of my current state x1 415 00:22:02,540 --> 00:22:08,390 plus the function 416 00:22:05,460 --> 00:22:12,150 you know multiplied by the time step H 417 00:22:08,390 --> 00:22:16,470 and that function is evaluated at the 418 00:22:12,150 --> 00:22:18,870 current state t1 and x1 and this method 419 00:22:16,470 --> 00:22:21,810 is called explicit Euler it's called 420 00:22:18,870 --> 00:22:24,630 explicit Euler because one boiler I 421 00:22:21,810 --> 00:22:27,450 think he invented it but it's also 422 00:22:24,630 --> 00:22:30,810 called it's it called explicit because I 423 00:22:27,450 --> 00:22:32,760 can explicitly just plug in values I 424 00:22:30,810 --> 00:22:34,530 already have and get the result back 425 00:22:32,760 --> 00:22:37,760 there's nothing enough to solve at this 426 00:22:34,530 --> 00:22:37,760 point it's just plug and chug 427 00:22:39,039 --> 00:22:46,610 so here's explicit Euler applied to 428 00:22:42,410 --> 00:22:48,049 Newton's law so remember I converted 429 00:22:46,610 --> 00:22:50,270 Newton's law into two first order 430 00:22:48,049 --> 00:22:53,570 differential equations so I have a first 431 00:22:50,270 --> 00:22:55,370 type of velocity update so I compute the 432 00:22:53,570 --> 00:22:59,750 new velocity in terms of the current 433 00:22:55,370 --> 00:23:02,120 velocity and the current force and then 434 00:22:59,750 --> 00:23:04,070 I also update the position in terms of 435 00:23:02,120 --> 00:23:08,809 the current position and the current 436 00:23:04,070 --> 00:23:11,120 velocity and if you look at this you can 437 00:23:08,809 --> 00:23:15,830 see this is kind of like extrapolation 438 00:23:11,120 --> 00:23:18,460 so that's gonna be important but anyhow 439 00:23:15,830 --> 00:23:21,980 here here's how that looks in pseudocode 440 00:23:18,460 --> 00:23:26,450 so I just formed some state structure it 441 00:23:21,980 --> 00:23:30,140 has time position velocity and then I 442 00:23:26,450 --> 00:23:32,059 have my current state I called the 443 00:23:30,140 --> 00:23:35,150 explicit order function I pass in the 444 00:23:32,059 --> 00:23:36,860 time step H as well you know then I use 445 00:23:35,150 --> 00:23:38,419 that state to compute the force whatever 446 00:23:36,860 --> 00:23:41,210 that force may be it may be springs 447 00:23:38,419 --> 00:23:45,020 dampers friction who-knows-what wind 448 00:23:41,210 --> 00:23:46,970 whatever and it just shows that all into 449 00:23:45,020 --> 00:23:49,220 that flow here I'm just doing one 450 00:23:46,970 --> 00:23:51,140 dimensional motion but you could have 451 00:23:49,220 --> 00:23:54,200 three dimensional motion it's two 452 00:23:51,140 --> 00:23:55,490 dimensional motion it would all work you 453 00:23:54,200 --> 00:24:01,880 just add vectors instead of adding 454 00:23:55,490 --> 00:24:05,960 scalars and then I apply the explicit 455 00:24:01,880 --> 00:24:08,780 euler so compute the new time compute 456 00:24:05,960 --> 00:24:11,890 the new velocity and then the new 457 00:24:08,780 --> 00:24:14,600 position and then return the new state 458 00:24:11,890 --> 00:24:19,580 so it's pretty straightforward very easy 459 00:24:14,600 --> 00:24:20,809 to implement in practice all right so we 460 00:24:19,580 --> 00:24:23,210 had that mass spring system 461 00:24:20,809 --> 00:24:26,990 so let's look explicit order for that 462 00:24:23,210 --> 00:24:32,360 mass spring system so here the force is 463 00:24:26,990 --> 00:24:34,690 coming in as a mega squared x1 divided 464 00:24:32,360 --> 00:24:34,690 by the mass 465 00:24:34,990 --> 00:24:39,220 all right so that's the only change I 466 00:24:36,940 --> 00:24:43,210 made so I'm just now I have an explicit 467 00:24:39,220 --> 00:24:45,760 force it's not generic anymore so here 468 00:24:43,210 --> 00:24:49,559 is the phase portrait for that so the 469 00:24:45,760 --> 00:24:51,909 exact solution is that red dashed circle 470 00:24:49,559 --> 00:24:55,960 that's that's that phase portrait I 471 00:24:51,909 --> 00:24:57,460 showed earlier and then if I did a 472 00:24:55,960 --> 00:24:59,710 numerical simulation with the explicit 473 00:24:57,460 --> 00:25:03,190 Euler and you can see the solution is 474 00:24:59,710 --> 00:25:05,590 that blue line spiraling outwards all 475 00:25:03,190 --> 00:25:07,720 right this is not stable this thing is 476 00:25:05,590 --> 00:25:11,440 blowing up actually it's just getting 477 00:25:07,720 --> 00:25:16,899 moving further and further as time goes 478 00:25:11,440 --> 00:25:19,270 on why is it unstable 479 00:25:16,899 --> 00:25:20,980 well I hinted this before you know 480 00:25:19,270 --> 00:25:25,419 explicit order is kind of like 481 00:25:20,980 --> 00:25:26,890 extrapolation if I if I'm at equal zero 482 00:25:25,419 --> 00:25:28,990 there you can see I have that dashed 483 00:25:26,890 --> 00:25:30,309 line that's my extrapolation but then 484 00:25:28,990 --> 00:25:35,289 the curve is shooting off in another 485 00:25:30,309 --> 00:25:37,480 direction well the force is well the 486 00:25:35,289 --> 00:25:39,789 sine wave is kind of a curve and if I'm 487 00:25:37,480 --> 00:25:41,730 trying to approximate the sine wave that 488 00:25:39,789 --> 00:25:43,929 the mass spring system would generate by 489 00:25:41,730 --> 00:25:47,830 extrapolating I'm gonna I'm gonna mean 490 00:25:43,929 --> 00:25:51,510 that a curve and so the actually the 491 00:25:47,830 --> 00:25:51,510 bigger the time step the worse it gets 492 00:25:53,669 --> 00:25:59,620 but we can we can actually analyze the 493 00:25:56,260 --> 00:26:03,580 stability of explicit or there a bit 494 00:25:59,620 --> 00:26:04,299 more precisely so here's how I'm gonna 495 00:26:03,580 --> 00:26:07,450 do that 496 00:26:04,299 --> 00:26:09,130 so I've assembled the the two first 497 00:26:07,450 --> 00:26:11,500 order differential equations kind of in 498 00:26:09,130 --> 00:26:13,510 an array here and if you look at it for 499 00:26:11,500 --> 00:26:17,380 a second on top on top there you can see 500 00:26:13,510 --> 00:26:19,210 whoa that's actually just like taking 501 00:26:17,380 --> 00:26:21,610 the current state and then multiplying 502 00:26:19,210 --> 00:26:23,950 by a matrix and then I get the new state 503 00:26:21,610 --> 00:26:26,409 so there and that in the bottom left I 504 00:26:23,950 --> 00:26:30,490 rearranged that into kind of a matrix 505 00:26:26,409 --> 00:26:32,679 vector form so I have this this matrix 506 00:26:30,490 --> 00:26:38,890 that just depends on the time step and 507 00:26:32,679 --> 00:26:41,080 Omega and then I just wrote that in 508 00:26:38,890 --> 00:26:44,020 vector notation X to the new state 509 00:26:41,080 --> 00:26:47,300 equals the matrix a times the current 510 00:26:44,020 --> 00:26:49,340 state x1 all right 511 00:26:47,300 --> 00:26:51,260 and then if you think about it well it's 512 00:26:49,340 --> 00:26:52,070 time goes on I'm doing time step after 513 00:26:51,260 --> 00:26:54,260 time step 514 00:26:52,070 --> 00:26:58,310 that's just like multiplying by that 515 00:26:54,260 --> 00:27:02,420 matrix a over and over again so I can 516 00:26:58,310 --> 00:27:04,550 get X 2 is a times X 1 X 3 is a squared 517 00:27:02,420 --> 00:27:08,870 times X 1 and I just keep multiplying 518 00:27:04,550 --> 00:27:11,210 mother that a matrix alright now if you 519 00:27:08,870 --> 00:27:13,970 think about a matrixes well that matrix 520 00:27:11,210 --> 00:27:18,320 can be have a kind of a magnifying 521 00:27:13,970 --> 00:27:19,970 influence or scaling down influence so 522 00:27:18,320 --> 00:27:21,560 if you think about the magnitude of the 523 00:27:19,970 --> 00:27:25,400 matrix well if the magnitude of that 524 00:27:21,560 --> 00:27:27,950 matrix is greater than 1 that means I'm 525 00:27:25,400 --> 00:27:31,130 kind of like scaling up that vector 526 00:27:27,950 --> 00:27:34,790 every time I multiply and so that would 527 00:27:31,130 --> 00:27:37,460 be unstable but if the magnitude is less 528 00:27:34,790 --> 00:27:43,520 than or equal to 1 then this solution 529 00:27:37,460 --> 00:27:48,380 would be stable so we can actually go 530 00:27:43,520 --> 00:27:50,810 further so the you can express the the 531 00:27:48,380 --> 00:27:55,940 magnitude of a matrix in terms of eigen 532 00:27:50,810 --> 00:27:57,350 values so there's you know there's lots 533 00:27:55,940 --> 00:27:58,730 of references for eigenvalues and 534 00:27:57,350 --> 00:28:01,310 eigenvectors I don't have time to go 535 00:27:58,730 --> 00:28:05,450 through this in detail but the 536 00:28:01,310 --> 00:28:08,270 eigenvalue and eigenvector have a very 537 00:28:05,450 --> 00:28:10,520 special property so if I if I find a 538 00:28:08,270 --> 00:28:13,550 specific direction so if I'm 2d I find 539 00:28:10,520 --> 00:28:17,080 some specific direction and then I 540 00:28:13,550 --> 00:28:17,080 multiply the matrix times that vector 541 00:28:17,200 --> 00:28:23,210 that is equivalent to multiplying that 542 00:28:19,970 --> 00:28:25,520 vector by a scalar the eigen value 543 00:28:23,210 --> 00:28:27,620 lambda so U is the eigen vector and 544 00:28:25,520 --> 00:28:29,390 lambda there that's nice I don't really 545 00:28:27,620 --> 00:28:31,160 care about the eigenvectors going 546 00:28:29,390 --> 00:28:35,570 forward I just care about the eigenvalue 547 00:28:31,160 --> 00:28:37,970 lambda so what I want to do is find this 548 00:28:35,570 --> 00:28:43,240 eigen value and then I know and I like 549 00:28:37,970 --> 00:28:46,130 what the magnitude of a matrix is so 550 00:28:43,240 --> 00:28:50,630 using the definition of the eigen value 551 00:28:46,130 --> 00:28:55,700 I can form this a linear algebra problem 552 00:28:50,630 --> 00:28:57,590 and then I know that well if so so let 553 00:28:55,700 --> 00:28:58,510 me say that that landed there that's a 554 00:28:57,590 --> 00:29:00,970 scalar 555 00:28:58,510 --> 00:29:04,570 and then I that's the identity matrix so 556 00:29:00,970 --> 00:29:07,690 in for backward I showed you that be the 557 00:29:04,570 --> 00:29:10,360 two by two I did any matrix and then 558 00:29:07,690 --> 00:29:13,299 there's my a matrix so I know that thing 559 00:29:10,360 --> 00:29:16,480 times U is zero so that means that the 560 00:29:13,299 --> 00:29:18,280 the rows must be linearly dependent so 561 00:29:16,480 --> 00:29:21,429 that must mean that the determinant is 562 00:29:18,280 --> 00:29:24,030 zero so that's the determinant down 563 00:29:21,429 --> 00:29:27,640 there on the bottom 564 00:29:24,030 --> 00:29:31,480 alright so bringing back in the mass 565 00:29:27,640 --> 00:29:33,100 spring system matrix there I can 566 00:29:31,480 --> 00:29:35,700 actually go ahead and compute that 567 00:29:33,100 --> 00:29:38,679 determinate and that results in a 568 00:29:35,700 --> 00:29:39,669 quadratic polynomial and that that 569 00:29:38,679 --> 00:29:43,510 polynomial is called the characteristic 570 00:29:39,669 --> 00:29:46,840 polynomial so I have a quadratic 571 00:29:43,510 --> 00:29:49,120 polynomial in lambda now and it depends 572 00:29:46,840 --> 00:29:53,669 on these constants I have the time step 573 00:29:49,120 --> 00:29:57,730 and Omega so I can actually compute 574 00:29:53,669 --> 00:29:59,770 lambda in terms of the time step in 575 00:29:57,730 --> 00:30:02,260 Omega using you know the quadratic 576 00:29:59,770 --> 00:30:04,330 formula and so I went ahead and do that 577 00:30:02,260 --> 00:30:06,340 and then I found that the magnitude of 578 00:30:04,330 --> 00:30:10,150 lambda is equal to square root of 1 plus 579 00:30:06,340 --> 00:30:12,250 h squared Omega squared well if you look 580 00:30:10,150 --> 00:30:15,160 at that for a second you can see well if 581 00:30:12,250 --> 00:30:17,140 H and Omega are greater than 1 then it's 582 00:30:15,160 --> 00:30:19,780 like we're sorry greater than zero then 583 00:30:17,140 --> 00:30:21,370 it's like 1 plus something square root 584 00:30:19,780 --> 00:30:24,610 of that well I know that's always bigger 585 00:30:21,370 --> 00:30:27,370 than 1 so that shows you that explicit 586 00:30:24,610 --> 00:30:32,260 Euler is always unstable for a mass 587 00:30:27,370 --> 00:30:34,179 spring system you can add some damping 588 00:30:32,260 --> 00:30:37,480 in your simulation try to get explicit 589 00:30:34,179 --> 00:30:39,100 order working but I'm going to show you 590 00:30:37,480 --> 00:30:42,750 a better method you don't even have to 591 00:30:39,100 --> 00:30:46,090 bother doing that all right 592 00:30:42,750 --> 00:30:47,950 let's return back to the original 593 00:30:46,090 --> 00:30:50,290 problem of solving differential 594 00:30:47,950 --> 00:30:52,360 equations numerically so instead of 595 00:30:50,290 --> 00:30:54,190 using the forward difference what if I 596 00:30:52,360 --> 00:30:57,340 use that backwards difference I showed 597 00:30:54,190 --> 00:31:01,179 you so there I put in that backwards 598 00:30:57,340 --> 00:31:06,299 difference x1 minus x0 divided by H and 599 00:31:01,179 --> 00:31:06,299 the functions Y evaluated at x1 600 00:31:07,500 --> 00:31:12,570 all right so that that succeeds in 601 00:31:10,700 --> 00:31:14,820 converting that differential equation 602 00:31:12,570 --> 00:31:17,669 into an algebraic equation but let's 603 00:31:14,820 --> 00:31:19,559 just say I want to I want to get to X - 604 00:31:17,669 --> 00:31:22,200 I don't want to I already have X 0 I'm 605 00:31:19,559 --> 00:31:24,780 at X 1 I want to get to X 2 from X 1 606 00:31:22,200 --> 00:31:27,090 well I can do that just by shifting the 607 00:31:24,780 --> 00:31:29,309 indices forward 1 so I add 1 to each of 608 00:31:27,090 --> 00:31:32,309 those indices and then I get this update 609 00:31:29,309 --> 00:31:36,240 formula in the bottom excuse X 1 plus 610 00:31:32,309 --> 00:31:36,990 the function evaluated at X 2 so that's 611 00:31:36,240 --> 00:31:44,520 different 612 00:31:36,990 --> 00:31:47,730 now the the function depends on the new 613 00:31:44,520 --> 00:31:49,500 state not the current state so I can't 614 00:31:47,730 --> 00:31:53,640 just plug and chug here I can't just 615 00:31:49,500 --> 00:31:55,470 compute X 2 and and move on because the 616 00:31:53,640 --> 00:31:57,900 X 2 appears on both sides of the 617 00:31:55,470 --> 00:31:59,580 equation so it it's actually I have to 618 00:31:57,900 --> 00:32:01,320 solve something now I have to solve an 619 00:31:59,580 --> 00:32:03,900 algebraic equation it could be linear it 620 00:32:01,320 --> 00:32:09,210 could be nonlinear we'll see examples of 621 00:32:03,900 --> 00:32:12,890 both so implicit order for the 622 00:32:09,210 --> 00:32:17,760 mass-spring system here's how it looks 623 00:32:12,890 --> 00:32:21,419 so I have V 2 depends on V 1 plus the 624 00:32:17,760 --> 00:32:24,289 influence of the spring at X 2 and then 625 00:32:21,419 --> 00:32:29,309 I have an update for X 2 that depends on 626 00:32:24,289 --> 00:32:32,309 the new V 2 so this is not too bad 627 00:32:29,309 --> 00:32:35,070 this is actually just a couple linear 628 00:32:32,309 --> 00:32:40,230 equations I can solve this I can get the 629 00:32:35,070 --> 00:32:41,970 new X 2 and V 2 it's not too bad but if 630 00:32:40,230 --> 00:32:46,860 you had a huge system this would be 631 00:32:41,970 --> 00:32:48,960 really expensive to do in a game so we 632 00:32:46,860 --> 00:32:50,760 generally don't apply this to large 633 00:32:48,960 --> 00:32:53,480 systems but anyhow here's how the 634 00:32:50,760 --> 00:32:56,640 solution looks with the phase portrait 635 00:32:53,480 --> 00:33:00,150 so it starts at the green circle and the 636 00:32:56,640 --> 00:33:02,909 solution de spirals inwards so it's 637 00:33:00,150 --> 00:33:04,860 stable this is great you know explicitly 638 00:33:02,909 --> 00:33:11,010 where's unstable now we have something 639 00:33:04,860 --> 00:33:13,799 it's stable and also there's some like 640 00:33:11,010 --> 00:33:15,179 there's some damping there's there's 641 00:33:13,799 --> 00:33:17,340 actually no damping in the mass spring 642 00:33:15,179 --> 00:33:19,530 system it should just go on that red - 643 00:33:17,340 --> 00:33:20,440 circle forever but here this thing is 644 00:33:19,530 --> 00:33:23,290 actually 645 00:33:20,440 --> 00:33:25,660 it's damned so this is called numerical 646 00:33:23,290 --> 00:33:28,750 damping so implicit order it kind of 647 00:33:25,660 --> 00:33:31,360 brings along its own damping and this is 648 00:33:28,750 --> 00:33:33,640 generally not a problem for games it's 649 00:33:31,360 --> 00:33:36,550 okay if things dampen games off most of 650 00:33:33,640 --> 00:33:39,630 the time all right so I went ahead and 651 00:33:36,550 --> 00:33:44,370 and did the eigen value computation and 652 00:33:39,630 --> 00:33:48,700 here's what I found for implicit order 653 00:33:44,370 --> 00:33:50,790 the eigen value is 1 over square root of 654 00:33:48,700 --> 00:33:54,610 1 plus h squared Omega squared so 655 00:33:50,790 --> 00:33:57,220 positive H and Omega that means that the 656 00:33:54,610 --> 00:34:01,000 denominator is greater than 1 so that 657 00:33:57,220 --> 00:34:05,020 means the magnitude that lambda there is 658 00:34:01,000 --> 00:34:09,310 is always less than 1 so that shows the 659 00:34:05,020 --> 00:34:13,000 implicit where they're stable all right 660 00:34:09,310 --> 00:34:15,669 so I had this really nice cheap method 661 00:34:13,000 --> 00:34:20,080 explicit Euler completely unstable it's 662 00:34:15,669 --> 00:34:23,379 not usable then I had an implicit order 663 00:34:20,080 --> 00:34:25,210 it's like super stable but also it's 664 00:34:23,379 --> 00:34:27,399 super expensive to compute I don't put 665 00:34:25,210 --> 00:34:33,730 that he's like my general solver in my 666 00:34:27,399 --> 00:34:36,700 game so what do we do it'd be nice to 667 00:34:33,730 --> 00:34:38,980 have kind of a balance this can I have 668 00:34:36,700 --> 00:34:41,169 like the cheapness of explicit order but 669 00:34:38,980 --> 00:34:43,240 some of the stability maybe not all the 670 00:34:41,169 --> 00:34:45,250 stability of implicit were there and 671 00:34:43,240 --> 00:34:48,490 there it turns out there's a way to do 672 00:34:45,250 --> 00:34:51,940 that so here I'm gonna make a small 673 00:34:48,490 --> 00:34:53,889 modification to explicit Euler so if you 674 00:34:51,940 --> 00:34:56,649 look at explicit order there first you 675 00:34:53,889 --> 00:34:58,090 see there's the velocity update I get 676 00:34:56,649 --> 00:35:00,640 that new velocity in terms of their 677 00:34:58,090 --> 00:35:03,460 current velocity and current Forest so 678 00:35:00,640 --> 00:35:05,560 that's this plug in charge and then when 679 00:35:03,460 --> 00:35:08,170 I go do the position update well I could 680 00:35:05,560 --> 00:35:10,570 use the current velocity or I could use 681 00:35:08,170 --> 00:35:13,390 I could use that new velocity and I just 682 00:35:10,570 --> 00:35:15,370 computed well that's over and what I'm 683 00:35:13,390 --> 00:35:17,410 doing there over on the right is I'm 684 00:35:15,370 --> 00:35:19,600 gonna use that new velocity in my 685 00:35:17,410 --> 00:35:22,150 position update and that actually makes 686 00:35:19,600 --> 00:35:23,860 a world of difference it seems it seems 687 00:35:22,150 --> 00:35:25,960 kludgy it seems hacky but there's 688 00:35:23,860 --> 00:35:27,790 actually a lot of good theory and 689 00:35:25,960 --> 00:35:30,600 developed around this and you can find 690 00:35:27,790 --> 00:35:30,600 that in the references 691 00:35:32,530 --> 00:35:38,420 so I went ahead and I did the the eigen 692 00:35:35,420 --> 00:35:41,420 value computation it's kind of tricky 693 00:35:38,420 --> 00:35:45,849 but I worked it out and it turns out 694 00:35:41,420 --> 00:35:50,000 that semi implicit Euler is 695 00:35:45,849 --> 00:35:53,750 conditionally stable so the the eigen 696 00:35:50,000 --> 00:35:57,680 value is equal to one the magnitude of 697 00:35:53,750 --> 00:36:01,340 it is equal to 1 as long as H times 698 00:35:57,680 --> 00:36:03,170 Omega is less than equal to 2 that might 699 00:36:01,340 --> 00:36:07,190 seem really restrictive but keep in mind 700 00:36:03,170 --> 00:36:10,250 games H is often a thirtieth of a second 701 00:36:07,190 --> 00:36:12,560 or a sixtieth of a second so Omega can 702 00:36:10,250 --> 00:36:15,580 actually get to be fairly significant so 703 00:36:12,560 --> 00:36:18,290 you can't handle Springs with this and 704 00:36:15,580 --> 00:36:23,440 so here are some examples of some 705 00:36:18,290 --> 00:36:27,710 implicit order so here's here's the 706 00:36:23,440 --> 00:36:30,050 phase portrait you can see it's tracking 707 00:36:27,710 --> 00:36:33,500 that exact solution very well it's a 708 00:36:30,050 --> 00:36:36,320 little bit off of it but what is really 709 00:36:33,500 --> 00:36:38,480 cool is that it is actually coming back 710 00:36:36,320 --> 00:36:41,480 on itself it's not damped and it's not 711 00:36:38,480 --> 00:36:44,990 blowing up it's actually conserving 712 00:36:41,480 --> 00:36:47,930 energy on average and here's an example 713 00:36:44,990 --> 00:36:50,599 where I increased Omega from 1 to 4 you 714 00:36:47,930 --> 00:36:53,570 can see it's still stable the the 715 00:36:50,599 --> 00:36:59,300 accuracy is diminishing so here's an 716 00:36:53,570 --> 00:37:01,700 animation of that with omega equals 4 so 717 00:36:59,300 --> 00:37:07,580 it is you know kind of a nice looking 718 00:37:01,700 --> 00:37:09,500 sine wave increase the angular frequency 719 00:37:07,580 --> 00:37:12,500 to 10 which puts me at half the 720 00:37:09,500 --> 00:37:14,990 stability limit so H Omega equals 1 and 721 00:37:12,500 --> 00:37:19,099 I got this really strange result it's 722 00:37:14,990 --> 00:37:25,640 like this very geometric shape there so 723 00:37:19,099 --> 00:37:33,980 I did the animation for that is pretty 724 00:37:25,640 --> 00:37:39,260 wild you know increasing I make it even 725 00:37:33,980 --> 00:37:41,839 more no no I'm a goes 16 I get it's kind 726 00:37:39,260 --> 00:37:43,369 of getting kind of crazy it's kind of 727 00:37:41,839 --> 00:37:44,630 like some sort of Spira graph but it's 728 00:37:43,369 --> 00:37:47,000 still stable 729 00:37:44,630 --> 00:37:49,730 it may not be the solution you want your 730 00:37:47,000 --> 00:37:53,210 games maybe it's a little wobbly or 731 00:37:49,730 --> 00:37:55,850 something okay but to show you that same 732 00:37:53,210 --> 00:37:59,930 implicit warrior can blow up I increased 733 00:37:55,850 --> 00:38:03,220 a mega to 24 so now HM it's 2.4 yeah it 734 00:37:59,930 --> 00:38:06,080 blows up it does have its limits 735 00:38:03,220 --> 00:38:10,640 alright so that's all the numerical 736 00:38:06,080 --> 00:38:12,680 methods I'm going to show there's lots 737 00:38:10,640 --> 00:38:15,260 of numerical methods out there that you 738 00:38:12,680 --> 00:38:17,510 can use in games and you may be tempted 739 00:38:15,260 --> 00:38:20,990 like oh I can use Rangga CUDA because 740 00:38:17,510 --> 00:38:22,910 it's more accurate most of the time 99% 741 00:38:20,990 --> 00:38:25,970 of the time you don't need the accuracy 742 00:38:22,910 --> 00:38:28,370 that those other methods might offer and 743 00:38:25,970 --> 00:38:32,120 they're more to be potentially much more 744 00:38:28,370 --> 00:38:36,710 expensive to implement in terms of 745 00:38:32,120 --> 00:38:38,600 processor and memory so and then of 746 00:38:36,710 --> 00:38:41,480 course you should just completely avoid 747 00:38:38,600 --> 00:38:43,190 explicit order these are their methods 748 00:38:41,480 --> 00:38:44,690 they're you know it's there's obviously 749 00:38:43,190 --> 00:38:46,580 no harm in knowing about them and maybe 750 00:38:44,690 --> 00:38:48,680 they'll be useful at some point but I I 751 00:38:46,580 --> 00:38:51,440 recommend your first choice should be 752 00:38:48,680 --> 00:38:53,870 some implicit order it's fast and it's 753 00:38:51,440 --> 00:38:55,640 pretty stable and it you know has some 754 00:38:53,870 --> 00:38:59,570 nice properties like that it kind of 755 00:38:55,640 --> 00:39:01,160 conserves energy all right so now I'm 756 00:38:59,570 --> 00:39:04,570 going to look at a more advanced topic 757 00:39:01,160 --> 00:39:07,850 and I'm gonna show you how I applied 758 00:39:04,570 --> 00:39:11,120 these methods to a kind of difficult 759 00:39:07,850 --> 00:39:17,660 problem that I ran into in in games and 760 00:39:11,120 --> 00:39:19,700 that's under a gyroscopic motion alright 761 00:39:17,660 --> 00:39:22,460 so getting a little bit advanced here so 762 00:39:19,700 --> 00:39:25,130 rigid body rotation I don't have time to 763 00:39:22,460 --> 00:39:27,320 go into all the details of this but I'm 764 00:39:25,130 --> 00:39:30,290 gonna make a short attempt at that 765 00:39:27,320 --> 00:39:33,910 so this equation up here on the left 766 00:39:30,290 --> 00:39:38,090 that's basically like Newton's law for 767 00:39:33,910 --> 00:39:39,650 rotation so there's also some notational 768 00:39:38,090 --> 00:39:42,730 challenges here because rigid body 769 00:39:39,650 --> 00:39:45,650 dynamics has a really bad notation so I 770 00:39:42,730 --> 00:39:47,750 here is not the identity matrix it's the 771 00:39:45,650 --> 00:39:50,090 inertia tensor the inertia tensor is a 772 00:39:47,750 --> 00:39:51,540 three by three matrix so this equation 773 00:39:50,090 --> 00:39:56,250 up here is a 774 00:39:51,540 --> 00:39:59,180 vector equation so I the 3x3 matrix it's 775 00:39:56,250 --> 00:40:03,000 symmetric and it represents the 776 00:39:59,180 --> 00:40:06,960 rotational inertia of an object of a 777 00:40:03,000 --> 00:40:08,370 rigid body then Omega here is no longer 778 00:40:06,960 --> 00:40:11,850 than any other frequency it's the 779 00:40:08,370 --> 00:40:18,180 angular velocity and it is a three 780 00:40:11,850 --> 00:40:22,230 vector so you can and then there's this 781 00:40:18,180 --> 00:40:23,970 T this applied torque okay so you could 782 00:40:22,230 --> 00:40:26,730 kind of see that like as mass times 783 00:40:23,970 --> 00:40:28,380 acceleration equals force except there's 784 00:40:26,730 --> 00:40:32,520 this weird term in the middle 785 00:40:28,380 --> 00:40:36,090 Omega across I Omega and that is the 786 00:40:32,520 --> 00:40:38,430 gyroscopic torque that's something that 787 00:40:36,090 --> 00:40:40,460 only shows up and rotational dynamics 788 00:40:38,430 --> 00:40:44,460 you don't get this weird kind of term in 789 00:40:40,460 --> 00:40:47,040 linear motion and so I'm gonna try try 790 00:40:44,460 --> 00:40:50,940 to explain briefly where that gyroscopic 791 00:40:47,040 --> 00:40:53,100 torque comes from so in linear motion 792 00:40:50,940 --> 00:40:56,330 you have linear momentum which is mass 793 00:40:53,100 --> 00:40:59,460 times velocity okay 794 00:40:56,330 --> 00:41:02,520 likewise in angular motion you have 795 00:40:59,460 --> 00:41:06,030 angular momentum which is the inertia 796 00:41:02,520 --> 00:41:08,160 tensor times the angular velocity so 797 00:41:06,030 --> 00:41:12,060 that's L down there on the bottom right 798 00:41:08,160 --> 00:41:13,820 that's the angular momentum so one thing 799 00:41:12,060 --> 00:41:18,230 that's interesting that shows up in 800 00:41:13,820 --> 00:41:20,970 rigid body rotation is that my angular 801 00:41:18,230 --> 00:41:23,430 velocity my point in a different 802 00:41:20,970 --> 00:41:29,870 direction than my angular momentum and 803 00:41:23,430 --> 00:41:32,880 this creates discharge scopic torque 804 00:41:29,870 --> 00:41:34,740 because the object is moving in a 805 00:41:32,880 --> 00:41:39,600 different way than it wants to be kind 806 00:41:34,740 --> 00:41:44,520 of like that so let's look at an example 807 00:41:39,600 --> 00:41:47,610 of that so once I had this box and every 808 00:41:44,520 --> 00:41:49,350 side of that box is it different with so 809 00:41:47,610 --> 00:41:52,230 if I compute the rotational inertia 810 00:41:49,350 --> 00:41:55,290 about each of those three axes I'm gonna 811 00:41:52,230 --> 00:41:58,920 get a different value and those so I 812 00:41:55,290 --> 00:42:00,330 have ixx iyy is easy and when I build 813 00:41:58,920 --> 00:42:00,990 the inertia tensor that's a diagonal 814 00:42:00,330 --> 00:42:03,540 matrix 815 00:42:00,990 --> 00:42:05,880 but each element on the diagonal is 816 00:42:03,540 --> 00:42:08,320 different 817 00:42:05,880 --> 00:42:13,210 and then when I compute the angular 818 00:42:08,320 --> 00:42:16,150 momentum that's I times Omega that is 819 00:42:13,210 --> 00:42:21,130 basically like applying non-uniform 820 00:42:16,150 --> 00:42:22,390 scale to Omega so then that sends that 821 00:42:21,130 --> 00:42:25,210 creates a vector that points in a 822 00:42:22,390 --> 00:42:30,310 different direction than Omega so when I 823 00:42:25,210 --> 00:42:33,940 go compute Omega cross Omega which is 824 00:42:30,310 --> 00:42:37,270 basically the same as Omega cross L well 825 00:42:33,940 --> 00:42:38,830 that cross product is nonzero so the 826 00:42:37,270 --> 00:42:41,610 only case where that cross product turns 827 00:42:38,830 --> 00:42:45,340 out to be 0 is if you have a sphere and 828 00:42:41,610 --> 00:42:54,720 the rotational inertia is the same about 829 00:42:45,340 --> 00:42:57,460 any axis all right so let's go ahead and 830 00:42:54,720 --> 00:43:01,900 solve that rigid body rotation using 831 00:42:57,460 --> 00:43:03,730 semi implicit order so the first first 832 00:43:01,900 --> 00:43:05,650 equation there is the velocity update 833 00:43:03,730 --> 00:43:08,920 and it's it's basically what you've seen 834 00:43:05,650 --> 00:43:12,250 before I have the current angular 835 00:43:08,920 --> 00:43:14,020 velocity the time step the kind of the 836 00:43:12,250 --> 00:43:16,110 inverse mass which is the inverse 837 00:43:14,020 --> 00:43:20,800 inertia tensor and then I have my forces 838 00:43:16,110 --> 00:43:24,610 the second equation there that is when 839 00:43:20,800 --> 00:43:26,800 we get into rotation we don't have such 840 00:43:24,610 --> 00:43:28,810 a simple way of representing rotation we 841 00:43:26,800 --> 00:43:29,730 have we're using the quarter knee in 842 00:43:28,810 --> 00:43:32,530 here 843 00:43:29,730 --> 00:43:34,780 so the quarter Nina has to get updated 844 00:43:32,530 --> 00:43:37,840 in a kind of a funky way I put a 845 00:43:34,780 --> 00:43:42,610 reference for this update formula in the 846 00:43:37,840 --> 00:43:43,720 references but nevertheless that that's 847 00:43:42,610 --> 00:43:48,070 how you update your car turning 848 00:43:43,720 --> 00:43:50,230 according to the angular velocity and 849 00:43:48,070 --> 00:43:52,090 notice so I'm using the new angular 850 00:43:50,230 --> 00:43:57,070 velocity to update the quaternion 851 00:43:52,090 --> 00:44:00,190 okay so how does that work out so here 852 00:43:57,070 --> 00:44:02,440 what I'm going to do is I have a testbed 853 00:44:00,190 --> 00:44:04,960 for physics engine I use it lizard 854 00:44:02,440 --> 00:44:06,940 called Domino and I'm gonna drop a 855 00:44:04,960 --> 00:44:12,030 capsule that's spinning about the long 856 00:44:06,940 --> 00:44:12,030 axis and here's here's the simulation 857 00:44:12,810 --> 00:44:16,740 all right it just went unstable right 858 00:44:15,010 --> 00:44:24,940 away 859 00:44:16,740 --> 00:44:26,260 it's not working out so what happened I 860 00:44:24,940 --> 00:44:29,470 thought semi implicit where they were 861 00:44:26,260 --> 00:44:36,460 supposed to take care of me and make all 862 00:44:29,470 --> 00:44:38,590 my stuff stable but the problem is the 863 00:44:36,460 --> 00:44:42,310 velocity update the angular velocity 864 00:44:38,590 --> 00:44:44,200 update it's basically like by itself 865 00:44:42,310 --> 00:44:46,450 it's basically an explicit update it's 866 00:44:44,200 --> 00:44:50,470 basically its extrapolation it's like 867 00:44:46,450 --> 00:44:54,670 explicit Euler and if you look there 868 00:44:50,470 --> 00:44:57,310 you can see that I'm updating the 869 00:44:54,670 --> 00:45:01,350 angular velocity using a force the 870 00:44:57,310 --> 00:45:05,620 gyroscopic force which is quadratic in 871 00:45:01,350 --> 00:45:07,990 the current angular velocity so that's 872 00:45:05,620 --> 00:45:12,820 kind of a really you know steep function 873 00:45:07,990 --> 00:45:19,300 and I'm trying to extrapolate and it's 874 00:45:12,820 --> 00:45:22,180 it's not working out so exactly it's 875 00:45:19,300 --> 00:45:27,700 that term there that is quadratic and 876 00:45:22,180 --> 00:45:30,880 causing the instability so what I did 877 00:45:27,700 --> 00:45:33,040 for many years and a lot of physics 878 00:45:30,880 --> 00:45:36,460 programmers do this is we just throw out 879 00:45:33,040 --> 00:45:39,630 this gyroscopic torque it's just causing 880 00:45:36,460 --> 00:45:44,370 me problems I got a ship get rid of it 881 00:45:39,630 --> 00:45:44,370 so that's how it looks on the top there 882 00:45:45,570 --> 00:45:51,630 some people also like to try to do 883 00:45:49,960 --> 00:45:56,290 something a little bit more clever is 884 00:45:51,630 --> 00:45:59,380 well if I don't use Omega angular 885 00:45:56,290 --> 00:46:02,890 velocity as my state I could use the 886 00:45:59,380 --> 00:46:05,380 angular momentum L is my state and then 887 00:46:02,890 --> 00:46:08,680 the rigid body rotation equation just 888 00:46:05,380 --> 00:46:10,540 becomes L dot equals T and the 889 00:46:08,680 --> 00:46:14,350 gyroscopic torque doesn't actually 890 00:46:10,540 --> 00:46:18,160 appear and and that becomes that's a 891 00:46:14,350 --> 00:46:21,490 much simpler equation so here's how that 892 00:46:18,160 --> 00:46:25,570 looks if I'm gonna use angular momentum 893 00:46:21,490 --> 00:46:28,090 as my state so my velocity update now 894 00:46:25,570 --> 00:46:30,090 becomes my angular momentum update I 895 00:46:28,090 --> 00:46:32,170 compute L to the Nu 896 00:46:30,090 --> 00:46:37,810 based on the current hanging momentum 897 00:46:32,170 --> 00:46:42,580 and the torque that should be t1 all 898 00:46:37,810 --> 00:46:44,740 right but to update my rotation I got to 899 00:46:42,580 --> 00:46:47,350 get I got to get the angular velocity to 900 00:46:44,740 --> 00:46:51,010 update the quaternion so to get the 901 00:46:47,350 --> 00:46:54,990 angular velocity I got to multiply the 902 00:46:51,010 --> 00:46:54,990 inverse inertia tensor times the 903 00:46:55,470 --> 00:47:04,570 velocity but there's a problem here 904 00:46:58,840 --> 00:47:08,350 I mean I'm using the old inertia tensor 905 00:47:04,570 --> 00:47:11,650 I want I don't I can't I can't get the 906 00:47:08,350 --> 00:47:14,710 new inertia tensor I can't get I to the 907 00:47:11,650 --> 00:47:17,680 reason is the inertia tensor depends on 908 00:47:14,710 --> 00:47:22,510 orientation so if I had it you know this 909 00:47:17,680 --> 00:47:24,640 object and spinning about that axis it's 910 00:47:22,510 --> 00:47:26,800 called that the Z well if I rotate it 911 00:47:24,640 --> 00:47:30,810 now the inertia around that axis has 912 00:47:26,800 --> 00:47:33,700 changed so near short answer depends on 913 00:47:30,810 --> 00:47:35,740 orientation the bent so usually what you 914 00:47:33,700 --> 00:47:38,470 do is you have you have the inertia 915 00:47:35,740 --> 00:47:40,480 tensor you you construct it in local 916 00:47:38,470 --> 00:47:42,730 coordinates like here's my box and it's 917 00:47:40,480 --> 00:47:44,920 aligned with the exact season there's 918 00:47:42,730 --> 00:47:47,110 there's my inertia about those axes and 919 00:47:44,920 --> 00:47:49,480 then to get that into world coordinates 920 00:47:47,110 --> 00:47:53,020 you have to multiply like on both sides 921 00:47:49,480 --> 00:47:54,790 by the rotation matrix not gonna go into 922 00:47:53,020 --> 00:47:59,350 detail about that but that's basically 923 00:47:54,790 --> 00:48:00,700 the gist of it ok so I'm not getting 924 00:47:59,350 --> 00:48:03,160 talking more about using angular 925 00:48:00,700 --> 00:48:05,500 momentum at state let's look at the 926 00:48:03,160 --> 00:48:10,890 repercussions of just dropping the 927 00:48:05,500 --> 00:48:17,680 gyroscopic term so this is from Diablo 928 00:48:10,890 --> 00:48:20,130 so this I would I didn't have audio no 929 00:48:17,680 --> 00:48:24,880 oh well let me play that again 930 00:48:20,130 --> 00:48:26,890 oops all right so this character is 931 00:48:24,880 --> 00:48:30,070 going ragdoll on he has a staff the 932 00:48:26,890 --> 00:48:31,780 staff is like a thin box and so the tech 933 00:48:30,070 --> 00:48:35,140 artist came to me and say hey this this 934 00:48:31,780 --> 00:48:37,930 the staff is just spinning forever I'm 935 00:48:35,140 --> 00:48:41,020 like why isn't it falling down and I 936 00:48:37,930 --> 00:48:41,860 knew it's because of that gyroscopic 937 00:48:41,020 --> 00:48:45,400 term that I was 938 00:48:41,860 --> 00:48:49,030 I was ignoring and like what am I gonna 939 00:48:45,400 --> 00:48:51,460 do so let's look at that and more 940 00:48:49,030 --> 00:48:53,170 precisely here here is the the capsule 941 00:48:51,460 --> 00:48:57,640 the same initial conditions as I showed 942 00:48:53,170 --> 00:49:00,160 you before went unstable so this is 943 00:48:57,640 --> 00:49:02,050 dropping that gyroscopic torque it's 944 00:49:00,160 --> 00:49:03,910 it's almost too stable now it's just 945 00:49:02,050 --> 00:49:05,710 spinning forever and ever I want this 946 00:49:03,910 --> 00:49:07,240 thing to fall down and it actually 947 00:49:05,710 --> 00:49:08,590 dropped it a little bit tilted it wasn't 948 00:49:07,240 --> 00:49:10,630 I didn't drop it exactly vertically 949 00:49:08,590 --> 00:49:14,260 actually I tried to help it out I gave 950 00:49:10,630 --> 00:49:18,010 it a little little tilt it didn't fall 951 00:49:14,260 --> 00:49:19,510 down so okay I'm gonna I'm gonna add 952 00:49:18,010 --> 00:49:23,950 some rolling resistance I had this nice 953 00:49:19,510 --> 00:49:28,600 rolling resistance thing so it got 954 00:49:23,950 --> 00:49:33,210 really strange it's completely 955 00:49:28,600 --> 00:49:39,910 unphysical so this is not working out I 956 00:49:33,210 --> 00:49:42,010 can't I can't ship that alright so I 957 00:49:39,910 --> 00:49:45,040 scratched my head for a while I talked 958 00:49:42,010 --> 00:49:47,440 to my friend Dirk for quite some time 959 00:49:45,040 --> 00:49:51,940 about this problem and so I had this 960 00:49:47,440 --> 00:49:53,770 idea like well I want to update the 961 00:49:51,940 --> 00:49:55,090 velocity I have all these forces or 962 00:49:53,770 --> 00:49:57,250 these torques I have the gyroscopic 963 00:49:55,090 --> 00:50:00,340 torque I have maybe some external force 964 00:49:57,250 --> 00:50:04,030 I got constraints each one of these 965 00:50:00,340 --> 00:50:06,370 torques is kind of wanting to modify the 966 00:50:04,030 --> 00:50:11,400 angular velocity it's creating a delta 967 00:50:06,370 --> 00:50:14,950 for the angular velocity what if I just 968 00:50:11,400 --> 00:50:18,550 treat each one of those extorts 969 00:50:14,950 --> 00:50:22,180 separately so I'll take this one torque 970 00:50:18,550 --> 00:50:25,300 and I'll use explicit solution for that 971 00:50:22,180 --> 00:50:29,050 or I'll take another one out and I'll 972 00:50:25,300 --> 00:50:32,470 use implicit so that's what I did so I 973 00:50:29,050 --> 00:50:36,700 broke it apart my applied torques I'm 974 00:50:32,470 --> 00:50:38,920 just using explicit but for the the 975 00:50:36,700 --> 00:50:40,840 gyroscopes of torque I'm gonna do 976 00:50:38,920 --> 00:50:45,400 implicit order because I know 977 00:50:40,840 --> 00:50:47,230 implicit order is super stable all right 978 00:50:45,400 --> 00:50:51,230 and then so I get those deltas and then 979 00:50:47,230 --> 00:50:55,820 I add those Delta's in to get the new 980 00:50:51,230 --> 00:50:57,410 angular velocity all right but I set 981 00:50:55,820 --> 00:51:01,430 myself up for a kind of difficult 982 00:50:57,410 --> 00:51:04,339 problem I have this implicit equation to 983 00:51:01,430 --> 00:51:05,900 solve because remember I you apply the 984 00:51:04,339 --> 00:51:07,579 backwards difference and that converts 985 00:51:05,900 --> 00:51:09,320 the difference equation into an 986 00:51:07,579 --> 00:51:13,839 algebraic equation but the algebra 987 00:51:09,320 --> 00:51:16,190 equation is implicit there's this 988 00:51:13,839 --> 00:51:20,270 there's something to be solved here so I 989 00:51:16,190 --> 00:51:24,560 have to solve for Omega 2 but Omega 2 is 990 00:51:20,270 --> 00:51:26,450 it's a period in this cross product so 991 00:51:24,560 --> 00:51:29,180 this is actually a nonlinear algebraic 992 00:51:26,450 --> 00:51:31,910 equation so I have some function of 993 00:51:29,180 --> 00:51:34,280 Omega 2 that's equal 0 that's kind of 994 00:51:31,910 --> 00:51:36,020 like a root finding problem now there's 995 00:51:34,280 --> 00:51:39,170 only a root finding problem it's three 996 00:51:36,020 --> 00:51:41,300 dimensional remember this that I the 997 00:51:39,170 --> 00:51:47,990 inertia tensor is a three by three 998 00:51:41,300 --> 00:51:52,099 matrix Omega is a three vector another 999 00:51:47,990 --> 00:51:54,320 problem I don't have I - I don't have 1000 00:51:52,099 --> 00:51:58,520 the next inertia tensor we talked about 1001 00:51:54,320 --> 00:52:02,000 this already there's a solution remember 1002 00:51:58,520 --> 00:52:03,770 I said in the you compute the inertia 1003 00:52:02,000 --> 00:52:05,060 tensor kind of in local coordinates well 1004 00:52:03,770 --> 00:52:07,550 you can actually solve this whole 1005 00:52:05,060 --> 00:52:10,010 problem in local coordinates and in 1006 00:52:07,550 --> 00:52:15,140 local coordinates the inertia tensor is 1007 00:52:10,010 --> 00:52:16,970 constant another thing I need to do is I 1008 00:52:15,140 --> 00:52:20,450 need to solve a three dimensional 1009 00:52:16,970 --> 00:52:22,430 nonlinear equation well for 1010 00:52:20,450 --> 00:52:23,720 one-dimensional equations there's this 1011 00:52:22,430 --> 00:52:27,109 well-known method called newton-raphson 1012 00:52:23,720 --> 00:52:29,619 which is a route-finding technique and 1013 00:52:27,109 --> 00:52:31,900 it involves computing the derivative and 1014 00:52:29,619 --> 00:52:34,490 dividing by the derivative been 1015 00:52:31,900 --> 00:52:37,730 iterating basically until things 1016 00:52:34,490 --> 00:52:38,660 converge you can do this also in 1017 00:52:37,730 --> 00:52:41,599 multi-dimensions 1018 00:52:38,660 --> 00:52:45,800 so instead this the derivative I have 1019 00:52:41,599 --> 00:52:48,670 this J which is the Jacobian and that's 1020 00:52:45,800 --> 00:52:50,720 basically the function 1021 00:52:48,670 --> 00:52:54,790 differentiate with respect to all the 1022 00:52:50,720 --> 00:52:54,790 different variables they exist so 1023 00:52:57,690 --> 00:53:03,940 anyway so if I work in in local 1024 00:53:01,510 --> 00:53:07,330 coordinates I can I can derive that 1025 00:53:03,940 --> 00:53:09,790 Jacobian so I did that and here's the 1026 00:53:07,330 --> 00:53:13,450 exact formula for that Jacobian using 1027 00:53:09,790 --> 00:53:14,590 the inertia tensor in body or local corn 1028 00:53:13,450 --> 00:53:18,270 so that's isube 1029 00:53:14,590 --> 00:53:21,790 and that involves this this Q matrix 1030 00:53:18,270 --> 00:53:24,670 which takes a vector and then returns a 1031 00:53:21,790 --> 00:53:30,760 skew symmetric matrix anyway that's the 1032 00:53:24,670 --> 00:53:32,890 formula so now I have a gyroscopic 1033 00:53:30,760 --> 00:53:36,750 torque solver using an implicit order 1034 00:53:32,890 --> 00:53:40,480 and here's the pseudocode for that so 1035 00:53:36,750 --> 00:53:46,420 first step convert angular velocity into 1036 00:53:40,480 --> 00:53:49,570 local coordinates second step compute 1037 00:53:46,420 --> 00:53:54,150 the air in the nonlinear equation I want 1038 00:53:49,570 --> 00:53:57,130 to solve then compute the Jacobian and 1039 00:53:54,150 --> 00:54:00,610 then so I the Jacobian is a three by 1040 00:53:57,130 --> 00:54:02,830 three matrix and I have to basically 1041 00:54:00,610 --> 00:54:05,290 invert by the Jacobian or just do a 1042 00:54:02,830 --> 00:54:07,570 solve directly and then that updates 1043 00:54:05,290 --> 00:54:10,600 that anger velocity and body coordinates 1044 00:54:07,570 --> 00:54:12,070 and then I convert the the angular 1045 00:54:10,600 --> 00:54:14,500 velocity back out in the world 1046 00:54:12,070 --> 00:54:16,710 coordinates and notice there's no kind 1047 00:54:14,500 --> 00:54:19,510 of iteration loop I found that just one 1048 00:54:16,710 --> 00:54:21,700 newton-raphson iteration was enough so 1049 00:54:19,510 --> 00:54:23,940 you could integrate more but I didn't 1050 00:54:21,700 --> 00:54:27,960 see any any and you gain from that 1051 00:54:23,940 --> 00:54:27,960 alright so here's the capsule 1052 00:54:30,660 --> 00:54:35,029 yeah so it falls down that's exactly 1053 00:54:34,289 --> 00:54:43,230 what I wanted 1054 00:54:35,029 --> 00:54:44,880 here it is with rolling resistance so 1055 00:54:43,230 --> 00:54:47,130 yeah you can see you even stops rolling 1056 00:54:44,880 --> 00:54:49,710 that's really really nice solution 1057 00:54:47,130 --> 00:54:53,220 that's the kind of stuff we like all 1058 00:54:49,710 --> 00:54:54,750 right so the gyroscopic torque is not 1059 00:54:53,220 --> 00:54:58,049 only important for getting kind of 1060 00:54:54,750 --> 00:55:00,779 correct behavior and you know stability 1061 00:54:58,049 --> 00:55:03,029 problems there's also some pretty cool 1062 00:55:00,779 --> 00:55:04,950 effects that you can pick up when you 1063 00:55:03,029 --> 00:55:08,670 simulate the gyroscopic torque correctly 1064 00:55:04,950 --> 00:55:10,259 you know we know that when things rotate 1065 00:55:08,670 --> 00:55:12,720 they wobble and they move all kinds of 1066 00:55:10,259 --> 00:55:15,269 crazy ways things can kinda get even 1067 00:55:12,720 --> 00:55:19,700 crazier in space so this is the nipa 1068 00:55:15,269 --> 00:55:22,710 cough effect reproduced on Space Station 1069 00:55:19,700 --> 00:55:24,869 so look at this object is spinning about 1070 00:55:22,710 --> 00:55:30,000 one axis flipping around and spinning 1071 00:55:24,869 --> 00:55:31,710 about the same axis so I wanted to see 1072 00:55:30,000 --> 00:55:33,450 if I could recreate this effect because 1073 00:55:31,710 --> 00:55:35,549 like if I can recreate these kind of 1074 00:55:33,450 --> 00:55:37,680 cool wobbly effects you know maybe I'll 1075 00:55:35,549 --> 00:55:44,700 have that are looking explosions in my 1076 00:55:37,680 --> 00:55:47,819 game all right so here's with explicit 1077 00:55:44,700 --> 00:55:50,130 integration so this is the the unstable 1078 00:55:47,819 --> 00:55:53,309 solution so you can see in the upper 1079 00:55:50,130 --> 00:55:55,740 left I have the kinetic energy and every 1080 00:55:53,309 --> 00:55:57,690 time the object flips around its it's 1081 00:55:55,740 --> 00:56:01,759 increasing you know it's reproducing the 1082 00:55:57,690 --> 00:56:01,759 effect very well quite it's unstable 1083 00:56:04,640 --> 00:56:12,790 all right just goes faster and faster 1084 00:56:09,010 --> 00:56:16,280 all right if I drop that gyroscopic term 1085 00:56:12,790 --> 00:56:17,720 it's very boring not know I don't have a 1086 00:56:16,280 --> 00:56:21,200 cool explosion now or whatever in my 1087 00:56:17,720 --> 00:56:22,940 game it's just when you drop the 1088 00:56:21,200 --> 00:56:24,440 gyroscopic term would you actually end 1089 00:56:22,940 --> 00:56:26,390 up doing is you're conserving angular 1090 00:56:24,440 --> 00:56:29,570 velocity instead of conserving angular 1091 00:56:26,390 --> 00:56:32,840 momentum and then here's the the 1092 00:56:29,570 --> 00:56:34,460 implicit solution so you can see the 1093 00:56:32,840 --> 00:56:36,470 kinetic energy up there it might be 1094 00:56:34,460 --> 00:56:38,900 sorry that it's too small but it's 1095 00:56:36,470 --> 00:56:40,580 actually it's actually damping so the 1096 00:56:38,900 --> 00:56:42,290 kinetic energy is going down but I think 1097 00:56:40,580 --> 00:56:46,210 this is OK for games I'm getting kind of 1098 00:56:42,290 --> 00:56:52,400 cool effects and things are damping down 1099 00:56:46,210 --> 00:56:53,960 that that's okay all right so that's 1100 00:56:52,400 --> 00:56:59,150 about all I have I might go through a 1101 00:56:53,960 --> 00:57:00,350 couple conclusions and things you want 1102 00:56:59,150 --> 00:57:04,160 to look for when you're looking at 1103 00:57:00,350 --> 00:57:05,600 numerical methods so first if you have a 1104 00:57:04,160 --> 00:57:08,030 numerical method 1105 00:57:05,600 --> 00:57:10,820 maybe you invented a new one right you 1106 00:57:08,030 --> 00:57:13,550 want to make sure that as the time step 1107 00:57:10,820 --> 00:57:15,320 gets smaller that you're approaching an 1108 00:57:13,550 --> 00:57:16,520 exact solution and it's good to have 1109 00:57:15,320 --> 00:57:18,860 something like the mass spring system 1110 00:57:16,520 --> 00:57:23,360 where you know the exact solution so you 1111 00:57:18,860 --> 00:57:26,440 can compare so if your numerical method 1112 00:57:23,360 --> 00:57:28,310 always returns 42 no matter what 1113 00:57:26,440 --> 00:57:33,710 differential equation you put in it's 1114 00:57:28,310 --> 00:57:36,440 it's not converging accuracy so accuracy 1115 00:57:33,710 --> 00:57:40,190 is important when you're comparing two 1116 00:57:36,440 --> 00:57:43,490 numerical methods for the same time step 1117 00:57:40,190 --> 00:57:45,140 is one method closer to the exact 1118 00:57:43,490 --> 00:57:47,600 solution than the other that's accuracy 1119 00:57:45,140 --> 00:57:50,030 for games we don't need a whole lot 1120 00:57:47,600 --> 00:57:53,840 accuracy first order accuracy is fine 1121 00:57:50,030 --> 00:57:55,520 and then of course stability we have to 1122 00:57:53,840 --> 00:57:57,980 make sure that our simulations don't 1123 00:57:55,520 --> 00:58:00,140 blow up so this is where it really pays 1124 00:57:57,980 --> 00:58:03,050 to have an understanding of what what 1125 00:58:00,140 --> 00:58:05,560 you have what's implemented and what are 1126 00:58:03,050 --> 00:58:08,810 its limitations and I encourage you to 1127 00:58:05,560 --> 00:58:10,310 test those limitations you know know 1128 00:58:08,810 --> 00:58:13,900 know how far you can push your 1129 00:58:10,310 --> 00:58:13,900 simulation before before it blows up 1130 00:58:14,740 --> 00:58:18,780 all right so here's some references and 1131 00:58:19,020 --> 00:58:23,140 yeah you'll be able to get this slides 1132 00:58:21,099 --> 00:58:25,480 at box2d org and you can reach me on 1133 00:58:23,140 --> 00:58:26,980 Twitter there's lots of references out 1134 00:58:25,480 --> 00:58:31,510 there for this stuff sorry if I don't 1135 00:58:26,980 --> 00:58:32,680 have everything also when I upload the 1136 00:58:31,510 --> 00:58:34,150 slides there's gonna be a bunch of 1137 00:58:32,680 --> 00:58:35,740 stuffs I had to skip over I just 1138 00:58:34,150 --> 00:58:37,150 couldn't fit it all in so you'll be able 1139 00:58:35,740 --> 00:58:40,060 to see that stuff too if you download 1140 00:58:37,150 --> 00:58:42,160 the slides and feel free if you want to 1141 00:58:40,060 --> 00:58:43,810 go to the next talk now go ahead if you 1142 00:58:42,160 --> 00:58:47,680 want to stick around ask a question 1143 00:58:43,810 --> 00:58:54,870 please come up to the mic thank you 1144 00:58:47,680 --> 00:58:54,870 [Applause] 1145 00:58:59,220 --> 00:59:01,280 you