all right so I'm gonna be talking about
numerical methods I'm Erin kado
I work at Blizzard Entertainment I do a
lot of physics work there you've
probably seen something you'll see a
video later I'll give you a little hint
about that
so numerical methods in game physics
were often dealing with really complex
problems we have lots of moving objects
friction contact constraints and these
situations are very complex and we have
to use numerical methods to solve these
situations in particular we need to use
numerical methods to get approximate
solutions to all the differential
equations that pop up as a result of
that so there are many numerical methods
available but it isn't always clear what
is the best numerical method to use what
is the right accuracy performance
trade-off and also once you've chosen
method you know what are the stability
limits of that method you know is this
method gonna blow up you know you don't
want to be about to ship your game and
all some things are blowing up and you
didn't know why so it's it's good to
have a good understanding of numerical
methods and that's kind of the goal of
this presentation is to give you a
better kind of background and way to
think about these things to help you
make better games all right I'm gonna
start out first with a little bit of
review of differential calculus so
differential calculus starts with the
derivative so here's the the definition
of the derivative so say I have some
function f that is a function of the x
coordinate and so that's this curve here
and I want to get the slope at any point
along that curve well you can imagine I
can approximate that by drawing the
secant through there sample of two
points and then I have a step with
eh and I draw a line through there and
that's kind of the approximate slope
well if we take that step size H and
make that smaller and smaller and the
limit as that goes to zero our
approximation of the slope at that point
X is improved and so that's kind of the
definition of the derivative and lots of
formulas for computing derivatives and
giving you the exact slope at any point
along a curve if those have been
developed and so I'll go through a few
examples of that so the derivative of a
line is well that's that's the slope of
the line so the slope of the line is the
same everywhere along the line so if I
have my my function for the line you
know y equals MX plus B the derivative
of that is just M the slope or you know
the rise over the over the run so what
if we have a like a curve like a
parabola well there that the slope is
not the same at every point so so we
need to use you know the ideas from
differential calculus to compute that
that slope so if I have this parabola
that's y equals ax squared where a is
just an arbitrary constant I take the
derivative of that and that's to ax so
the the slope varies linearly along the
parabola when it gets when we get two
differential equations the the
exponential function is very interesting
so the exponential function has this
really interesting property that if I
take its derivative I get back another
accident exponential so here I have y
equals e a X where a is yet again an
arbitrary constant and then I compute
the derivative of that
get a e to the ax so friend of mine
found this cute little cartoon so a wild
exponential function appeared you use
differentiate it is not very effective
because you get back another exponential
when we look at you know simulation the
sine wave is another function that's
very important if I compute the
derivative of a sine wave I get back
basically another sine wave shifted by
90 degrees so when I compute the
derivative of sine I get cosine which is
basically signed as shifted 90 degrees
all right so in videogame physics were
often dealing with functions of time so
I'm going to move now to having time is
my independent variable and then I'm
going to use X a lot as the dependent
variable so X you can think of it as the
position and then if I compute the time
derivative of the position I get the
velocity and then the time derivative of
velocity
I get the acceleration well since we're
doing this a lot writing all this dx/dt
stuff gets a little busy so it's a
standard notation to use the dot to
represent the derivative so X dot is DX
DT and I can have for the acceleration X
double dot so I'm gonna be using that
notation going forward all right a
little bit more review let's talk about
differential equations so the main
premise here is okay what if I had the
slope and I want to get back to the
original function why is that
interesting because you know Newton's
law as you'll see that kind of that
gives us the slope and when I get back
to the original function so in
particular we're going to be looking at
what's called a first-order differential
equation so here's the standard form for
that
to X dot so that's DX DT is equal to
some function of time and X itself so
what's this function well there's all
kinds of functions that we might come
upon but we can kind of talk about this
in a generic way but important point
here is that it is a first order
differential equation because it's just
involving X not X double dot so this is
kind of what we want to do we want to
have this formula for the derivative of
X and then get back X so that's that is
what it means to solve a differential
equation if that function is just a
function of time it's kind of
straightforward it's just like reverse
differentiation or integration if that
function depends on X then things get a
little bit more complex I'll show some
examples all right so I want to create a
bunch of examples and I was thinking
well what's the simplest differential
equation I could have well how about X
dot equals 0 so that means that the
derivative is zero everywhere that means
the slope is constant and it's zero and
so then the solution is just that the X
is just constant it's equal to the
initial value of X x0 pretty exciting
all right so let's look at something a
little bit more complex so X dot is now
equal to a nonzero constant a so the
solution to that is well that's a line
because before we saw well if you have a
constant slope and that unless that's
what a line has so x equals 80 plus x 0
where X 0 again is the initial value of
X at T equals 0
okay but as I said things see a little
bit more complex if X is on the
right-hand side so here X dot equals a X
so that means that the slope of X
depends on the value of X well it turns
out the solution to this is the
exponential so x equals x 0 e to the a T
so X 0 is the initial value that makes
sense
if I put in T equals 0 exponential e to
the 0 is 1 so X at T equals 0 is X 0 and
this is where for me it gets a little
bit weird that you don't really derive
the solution you kind of guess the
solution and prove that it's true this
is how calculus works so you can see
down here on the the bottom is the proof
I'll substitute in my solution X 0 into
that top equation take its derivative
and then kind of regroup things and say
oh that is that is equal to a times X so
that's the proof all right of course in
game physics were we're interested in
Newton's second law so mass times
acceleration equals force and in
particular if we look in one dimension x
is just positioned along one the
coordinate axis and X double dot equals
force so that's a scalar equation and
this is called a second order
differential equation because it
involves the second derivative of X all
right so let's look at some examples of
that so gravity this is a constant force
so the force is negative mg and pulling
this point mass down alright let's let's
try to solve that well since the right
hand side is constant I can just
basically get the solution through
reverse differentiation
so I divide through by the mass so that
shows you that the acceleration doesn't
depend on the mass so y double dot
equals negative mg and then I integrate
that once and I get y dot equals v-0
minus GT V zero is the initial velocity
and then I integrate that one more time
and I get this quadratic function for
the vertical position Y and that
involves the the initial position there
as well and so that's how you get
parabolic motion because it's this
parabolic quadratic function
alright so things get a little bit more
interesting with Newton's law when we
start to involve forces that depend on
position and I'm gonna be using this
example throughout the rest of the
presentation this mass spring system
so I'll gonna go over it and in a bit of
detail so I have this point mass m and
it can only move along the horizontal
axis and it's connected to ground by
this spring K and the way it works is
well the more I stretch that spring the
harder the force pulls back so the force
is actually negative KX convention if
you look on Wikipedia you're gonna see
it written this form so the the force
has been move over to the left hand side
so there's your KX there
all right so here's how you might expect
that system to move so you would expect
some nice oscillation all right so let's
try to solve this thing so first thing
I'm going to do is reduce the number of
parameters so if I divide through by the
mass then I can reduce down to one
parameter and I'm going to use this
parameter Omega which is the angular
frequency
so that's has units of radians per
second and that is defined as the square
root of K over m so notice that as the
spring stiffness K increases the angular
frequency also increases but as the mass
increases the angular frequency goes
down so here is the solution to that
it's a combination of a sine and cosine
and we have the initial position and
velocity in there I noticed that the the
angular frequency shows up inside the
the sine and cosine functions so that
kind of scales the number of
oscillations per second that you get so
here's a plot of the solution for the
position and for the velocity versus
time for some specific initial
conditions and angular frequency of one
so you can see it's just you know it's a
cosine and a sine so what's interesting
is like well if I'm going to be solving
lots of these things maybe I could have
another way to look at this data and I
it an interesting way to look at the
data is actually plotting the velocity
versus the position and in this case
time becomes a parameter on that curve
so in this case the solution starts at
that green circle and then goes around
in a clockwise fashion and the angular
frequency determines how quickly it goes
around that circle as time progresses
all right so I want to make this a
little bit clearer so I I created an
animation here and so what you're gonna
see is on the upper left there that's
going to be position versus time the
bottom left is velocity versus time and
then over on the right that's the phase
portrait which is velocity versus
position I'm sorry it's going it's going
counterclockwise
all right we're gonna be seeing more of
these all right so like I said in the
beginning and game physics we're dealing
with lots of moving objects lots of
complicated effects like friction and
constraints and we can't come up with
these nice pretty exact solutions so we
have to turn to numerical methods so
right off the bat when you when you
switch to numerical methods there's a
few things that you have to to deal with
first is you're gonna be sacrificing
accuracy you're not going to get the
exact sine wave of more
you may not get the exact parabolic
motion anymore for video games this is
almost always okay it's okay to
sacrifice accuracy for video games we
can sacrifice a lot of accuracy but a
nice thing that you pick up is with
numerical methods that they work very
well with time stepping so if you have a
game loop that is you know moving the
game forward by the time step delta-t
every frame well that's how are the
numerical methods work they just like to
move from from one solution to the next
solution by a small amount of time so
that works out well another thing that's
nice about numerical methods is they
they let you treat everything in a kind
of generic way so if I have a whole
bunch of differential equations I can
just slam them into an array and then
feed that off to my numerical method and
then it fits back a result so it can
kind of view your system that's kind of
a black box but it does require your
differential equations to be set up in a
specific way so numerical methods like
the differential equations to be set up
in a first-order form so you can have
just an array of first-order
differential equations they can be
mutually dependent they can be super
complex in terms of the the forces or
whatever there
but it should be in this format so
that's usually not a problem at all we
can't we can assemble this and this is
you know computers are good at this sort
of thing alright but there is a problem
so we're usually dealing with Newton's
law so mass times X double dot equals
have that's at least for one dimension
you could view that maybe as X is a
three dimensional quantity and forces
three dimensional that would work too
but the problem is that's a second-order
differential equation and we want to get
that into first-order form so that we
can and hand it off to a numerical
method so it's actually not that bad you
can you can do that and and and the way
you do it is by introducing another
state so instead of just having X my
position is my state I'm also gonna have
velocity as my state and so then I can
say okay the V dot the velocity the
derivative of the velocity which is
acceleration that's equal to F divided
by M and then I can introduce this kind
of trivial differential equation X dot
equals V it is trivial but is it
actually a valid differential equation
and then I can just stack those two
first order differential equations
together and so what I am able to do is
convert that one second order
differential equation into two first
order differential equations so it's
pretty simple to make this happen all
right so when we're getting into
numerical methods a good place to start
is looking at the finite difference so
what I want to do here is show you a way
to approximate the derivative with a
basically a secant like we saw before so
this is pretty much going back to that
first slide about calculus so I can
approximate the derivative by
say I have my current state x1 and then
I take a future state x2 compute the
difference and divide that by the time
step and that's an approximation of the
the slope at X 1 now you can see there
I've drawn with the dashed line the the
actual slope at X 1 and then the green
line that's my approximation so you can
see it's kind of crude it gets better as
the time step gets smaller obviously and
this is called the forward difference
you can also have the backwards
difference so say I'm gonna store my old
state x0 and then I'll take my current
state - my old state and divide that by
the time step and that's called the
backward difference you're gonna see
later how that how these impact our
solution ok but let's let's go ahead and
apply that forward difference
approximation so what I'm doing here is
I have my generic first order
differential equation X dot equals F of
T and X and I'm just gonna get rid of X
dot by by inserting my approximation for
X dot using the forward difference x2
minus x1 over H and what that does is
that turns that differential equation
into an algebraic equation so down on
the bottom there I just rearranged the
terms so now I can compute X to my new
state in terms of my current state x1
plus the function
you know multiplied by the time step H
and that function is evaluated at the
current state t1 and x1 and this method
is called explicit Euler it's called
explicit Euler because one boiler I
think he invented it but it's also
called it's it called explicit because I
can explicitly just plug in values I
already have and get the result back
there's nothing enough to solve at this
point it's just plug and chug
so here's explicit Euler applied to
Newton's law so remember I converted
Newton's law into two first order
differential equations so I have a first
type of velocity update so I compute the
new velocity in terms of the current
velocity and the current force and then
I also update the position in terms of
the current position and the current
velocity and if you look at this you can
see this is kind of like extrapolation
so that's gonna be important but anyhow
here here's how that looks in pseudocode
so I just formed some state structure it
has time position velocity and then I
have my current state I called the
explicit order function I pass in the
time step H as well you know then I use
that state to compute the force whatever
that force may be it may be springs
dampers friction who-knows-what wind
whatever and it just shows that all into
that flow here I'm just doing one
dimensional motion but you could have
three dimensional motion it's two
dimensional motion it would all work you
just add vectors instead of adding
scalars and then I apply the explicit
euler so compute the new time compute
the new velocity and then the new
position and then return the new state
so it's pretty straightforward very easy
to implement in practice all right so we
had that mass spring system
so let's look explicit order for that
mass spring system so here the force is
coming in as a mega squared x1 divided
by the mass
all right so that's the only change I
made so I'm just now I have an explicit
force it's not generic anymore so here
is the phase portrait for that so the
exact solution is that red dashed circle
that's that's that phase portrait I
showed earlier and then if I did a
numerical simulation with the explicit
Euler and you can see the solution is
that blue line spiraling outwards all
right this is not stable this thing is
blowing up actually it's just getting
moving further and further as time goes
on why is it unstable
well I hinted this before you know
explicit order is kind of like
extrapolation if I if I'm at equal zero
there you can see I have that dashed
line that's my extrapolation but then
the curve is shooting off in another
direction well the force is well the
sine wave is kind of a curve and if I'm
trying to approximate the sine wave that
the mass spring system would generate by
extrapolating I'm gonna I'm gonna mean
that a curve and so the actually the
bigger the time step the worse it gets
but we can we can actually analyze the
stability of explicit or there a bit
more precisely so here's how I'm gonna
do that
so I've assembled the the two first
order differential equations kind of in
an array here and if you look at it for
a second on top on top there you can see
whoa that's actually just like taking
the current state and then multiplying
by a matrix and then I get the new state
so there and that in the bottom left I
rearranged that into kind of a matrix
vector form so I have this this matrix
that just depends on the time step and
Omega and then I just wrote that in
vector notation X to the new state
equals the matrix a times the current
state x1 all right
and then if you think about it well it's
time goes on I'm doing time step after
time step
that's just like multiplying by that
matrix a over and over again so I can
get X 2 is a times X 1 X 3 is a squared
times X 1 and I just keep multiplying
mother that a matrix alright now if you
think about a matrixes well that matrix
can be have a kind of a magnifying
influence or scaling down influence so
if you think about the magnitude of the
matrix well if the magnitude of that
matrix is greater than 1 that means I'm
kind of like scaling up that vector
every time I multiply and so that would
be unstable but if the magnitude is less
than or equal to 1 then this solution
would be stable so we can actually go
further so the you can express the the
magnitude of a matrix in terms of eigen
values so there's you know there's lots
of references for eigenvalues and
eigenvectors I don't have time to go
through this in detail but the
eigenvalue and eigenvector have a very
special property so if I if I find a
specific direction so if I'm 2d I find
some specific direction and then I
multiply the matrix times that vector
that is equivalent to multiplying that
vector by a scalar the eigen value
lambda so U is the eigen vector and
lambda there that's nice I don't really
care about the eigenvectors going
forward I just care about the eigenvalue
lambda so what I want to do is find this
eigen value and then I know and I like
what the magnitude of a matrix is so
using the definition of the eigen value
I can form this a linear algebra problem
and then I know that well if so so let
me say that that landed there that's a
scalar
and then I that's the identity matrix so
in for backward I showed you that be the
two by two I did any matrix and then
there's my a matrix so I know that thing
times U is zero so that means that the
the rows must be linearly dependent so
that must mean that the determinant is
zero so that's the determinant down
there on the bottom
alright so bringing back in the mass
spring system matrix there I can
actually go ahead and compute that
determinate and that results in a
quadratic polynomial and that that
polynomial is called the characteristic
polynomial so I have a quadratic
polynomial in lambda now and it depends
on these constants I have the time step
and Omega so I can actually compute
lambda in terms of the time step in
Omega using you know the quadratic
formula and so I went ahead and do that
and then I found that the magnitude of
lambda is equal to square root of 1 plus
h squared Omega squared well if you look
at that for a second you can see well if
H and Omega are greater than 1 then it's
like we're sorry greater than zero then
it's like 1 plus something square root
of that well I know that's always bigger
than 1 so that shows you that explicit
Euler is always unstable for a mass
spring system you can add some damping
in your simulation try to get explicit
order working but I'm going to show you
a better method you don't even have to
bother doing that all right
let's return back to the original
problem of solving differential
equations numerically so instead of
using the forward difference what if I
use that backwards difference I showed
you so there I put in that backwards
difference x1 minus x0 divided by H and
the functions Y evaluated at x1
all right so that that succeeds in
converting that differential equation
into an algebraic equation but let's
just say I want to I want to get to X -
I don't want to I already have X 0 I'm
at X 1 I want to get to X 2 from X 1
well I can do that just by shifting the
indices forward 1 so I add 1 to each of
those indices and then I get this update
formula in the bottom excuse X 1 plus
the function evaluated at X 2 so that's
different
now the the function depends on the new
state not the current state so I can't
just plug and chug here I can't just
compute X 2 and and move on because the
X 2 appears on both sides of the
equation so it it's actually I have to
solve something now I have to solve an
algebraic equation it could be linear it
could be nonlinear we'll see examples of
both so implicit order for the
mass-spring system here's how it looks
so I have V 2 depends on V 1 plus the
influence of the spring at X 2 and then
I have an update for X 2 that depends on
the new V 2 so this is not too bad
this is actually just a couple linear
equations I can solve this I can get the
new X 2 and V 2 it's not too bad but if
you had a huge system this would be
really expensive to do in a game so we
generally don't apply this to large
systems but anyhow here's how the
solution looks with the phase portrait
so it starts at the green circle and the
solution de spirals inwards so it's
stable this is great you know explicitly
where's unstable now we have something
it's stable and also there's some like
there's some damping there's there's
actually no damping in the mass spring
system it should just go on that red -
circle forever but here this thing is
actually
it's damned so this is called numerical
damping so implicit order it kind of
brings along its own damping and this is
generally not a problem for games it's
okay if things dampen games off most of
the time all right so I went ahead and
and did the eigen value computation and
here's what I found for implicit order
the eigen value is 1 over square root of
1 plus h squared Omega squared so
positive H and Omega that means that the
denominator is greater than 1 so that
means the magnitude that lambda there is
is always less than 1 so that shows the
implicit where they're stable all right
so I had this really nice cheap method
explicit Euler completely unstable it's
not usable then I had an implicit order
it's like super stable but also it's
super expensive to compute I don't put
that he's like my general solver in my
game so what do we do it'd be nice to
have kind of a balance this can I have
like the cheapness of explicit order but
some of the stability maybe not all the
stability of implicit were there and
there it turns out there's a way to do
that so here I'm gonna make a small
modification to explicit Euler so if you
look at explicit order there first you
see there's the velocity update I get
that new velocity in terms of their
current velocity and current Forest so
that's this plug in charge and then when
I go do the position update well I could
use the current velocity or I could use
I could use that new velocity and I just
computed well that's over and what I'm
doing there over on the right is I'm
gonna use that new velocity in my
position update and that actually makes
a world of difference it seems it seems
kludgy it seems hacky but there's
actually a lot of good theory and
developed around this and you can find
that in the references
so I went ahead and I did the the eigen
value computation it's kind of tricky
but I worked it out and it turns out
that semi implicit Euler is
conditionally stable so the the eigen
value is equal to one the magnitude of
it is equal to 1 as long as H times
Omega is less than equal to 2 that might
seem really restrictive but keep in mind
games H is often a thirtieth of a second
or a sixtieth of a second so Omega can
actually get to be fairly significant so
you can't handle Springs with this and
so here are some examples of some
implicit order so here's here's the
phase portrait you can see it's tracking
that exact solution very well it's a
little bit off of it but what is really
cool is that it is actually coming back
on itself it's not damped and it's not
blowing up it's actually conserving
energy on average and here's an example
where I increased Omega from 1 to 4 you
can see it's still stable the the
accuracy is diminishing so here's an
animation of that with omega equals 4 so
it is you know kind of a nice looking
sine wave increase the angular frequency
to 10 which puts me at half the
stability limit so H Omega equals 1 and
I got this really strange result it's
like this very geometric shape there so
I did the animation for that is pretty
wild you know increasing I make it even
more no no I'm a goes 16 I get it's kind
of getting kind of crazy it's kind of
like some sort of Spira graph but it's
still stable
it may not be the solution you want your
games maybe it's a little wobbly or
something okay but to show you that same
implicit warrior can blow up I increased
a mega to 24 so now HM it's 2.4 yeah it
blows up it does have its limits
alright so that's all the numerical
methods I'm going to show there's lots
of numerical methods out there that you
can use in games and you may be tempted
like oh I can use Rangga CUDA because
it's more accurate most of the time 99%
of the time you don't need the accuracy
that those other methods might offer and
they're more to be potentially much more
expensive to implement in terms of
processor and memory so and then of
course you should just completely avoid
explicit order these are their methods
they're you know it's there's obviously
no harm in knowing about them and maybe
they'll be useful at some point but I I
recommend your first choice should be
some implicit order it's fast and it's
pretty stable and it you know has some
nice properties like that it kind of
conserves energy all right so now I'm
going to look at a more advanced topic
and I'm gonna show you how I applied
these methods to a kind of difficult
problem that I ran into in in games and
that's under a gyroscopic motion alright
so getting a little bit advanced here so
rigid body rotation I don't have time to
go into all the details of this but I'm
gonna make a short attempt at that
so this equation up here on the left
that's basically like Newton's law for
rotation so there's also some notational
challenges here because rigid body
dynamics has a really bad notation so I
here is not the identity matrix it's the
inertia tensor the inertia tensor is a
three by three matrix so this equation
up here is a
vector equation so I the 3x3 matrix it's
symmetric and it represents the
rotational inertia of an object of a
rigid body then Omega here is no longer
than any other frequency it's the
angular velocity and it is a three
vector so you can and then there's this
T this applied torque okay so you could
kind of see that like as mass times
acceleration equals force except there's
this weird term in the middle
Omega across I Omega and that is the
gyroscopic torque that's something that
only shows up and rotational dynamics
you don't get this weird kind of term in
linear motion and so I'm gonna try try
to explain briefly where that gyroscopic
torque comes from so in linear motion
you have linear momentum which is mass
times velocity okay
likewise in angular motion you have
angular momentum which is the inertia
tensor times the angular velocity so
that's L down there on the bottom right
that's the angular momentum so one thing
that's interesting that shows up in
rigid body rotation is that my angular
velocity my point in a different
direction than my angular momentum and
this creates discharge scopic torque
because the object is moving in a
different way than it wants to be kind
of like that so let's look at an example
of that so once I had this box and every
side of that box is it different with so
if I compute the rotational inertia
about each of those three axes I'm gonna
get a different value and those so I
have ixx iyy is easy and when I build
the inertia tensor that's a diagonal
matrix
but each element on the diagonal is
different
and then when I compute the angular
momentum that's I times Omega that is
basically like applying non-uniform
scale to Omega so then that sends that
creates a vector that points in a
different direction than Omega so when I
go compute Omega cross Omega which is
basically the same as Omega cross L well
that cross product is nonzero so the
only case where that cross product turns
out to be 0 is if you have a sphere and
the rotational inertia is the same about
any axis all right so let's go ahead and
solve that rigid body rotation using
semi implicit order so the first first
equation there is the velocity update
and it's it's basically what you've seen
before I have the current angular
velocity the time step the kind of the
inverse mass which is the inverse
inertia tensor and then I have my forces
the second equation there that is when
we get into rotation we don't have such
a simple way of representing rotation we
have we're using the quarter knee in
here
so the quarter Nina has to get updated
in a kind of a funky way I put a
reference for this update formula in the
references but nevertheless that that's
how you update your car turning
according to the angular velocity and
notice so I'm using the new angular
velocity to update the quaternion
okay so how does that work out so here
what I'm going to do is I have a testbed
for physics engine I use it lizard
called Domino and I'm gonna drop a
capsule that's spinning about the long
axis and here's here's the simulation
all right it just went unstable right
away
it's not working out so what happened I
thought semi implicit where they were
supposed to take care of me and make all
my stuff stable but the problem is the
velocity update the angular velocity
update it's basically like by itself
it's basically an explicit update it's
basically its extrapolation it's like
explicit Euler and if you look there
you can see that I'm updating the
angular velocity using a force the
gyroscopic force which is quadratic in
the current angular velocity so that's
kind of a really you know steep function
and I'm trying to extrapolate and it's
it's not working out so exactly it's
that term there that is quadratic and
causing the instability so what I did
for many years and a lot of physics
programmers do this is we just throw out
this gyroscopic torque it's just causing
me problems I got a ship get rid of it
so that's how it looks on the top there
some people also like to try to do
something a little bit more clever is
well if I don't use Omega angular
velocity as my state I could use the
angular momentum L is my state and then
the rigid body rotation equation just
becomes L dot equals T and the
gyroscopic torque doesn't actually
appear and and that becomes that's a
much simpler equation so here's how that
looks if I'm gonna use angular momentum
as my state so my velocity update now
becomes my angular momentum update I
compute L to the Nu
based on the current hanging momentum
and the torque that should be t1 all
right but to update my rotation I got to
get I got to get the angular velocity to
update the quaternion so to get the
angular velocity I got to multiply the
inverse inertia tensor times the
velocity but there's a problem here
I mean I'm using the old inertia tensor
I want I don't I can't I can't get the
new inertia tensor I can't get I to the
reason is the inertia tensor depends on
orientation so if I had it you know this
object and spinning about that axis it's
called that the Z well if I rotate it
now the inertia around that axis has
changed so near short answer depends on
orientation the bent so usually what you
do is you have you have the inertia
tensor you you construct it in local
coordinates like here's my box and it's
aligned with the exact season there's
there's my inertia about those axes and
then to get that into world coordinates
you have to multiply like on both sides
by the rotation matrix not gonna go into
detail about that but that's basically
the gist of it ok so I'm not getting
talking more about using angular
momentum at state let's look at the
repercussions of just dropping the
gyroscopic term so this is from Diablo
so this I would I didn't have audio no
oh well let me play that again
oops all right so this character is
going ragdoll on he has a staff the
staff is like a thin box and so the tech
artist came to me and say hey this this
the staff is just spinning forever I'm
like why isn't it falling down and I
knew it's because of that gyroscopic
term that I was
I was ignoring and like what am I gonna
do so let's look at that and more
precisely here here is the the capsule
the same initial conditions as I showed
you before went unstable so this is
dropping that gyroscopic torque it's
it's almost too stable now it's just
spinning forever and ever I want this
thing to fall down and it actually
dropped it a little bit tilted it wasn't
I didn't drop it exactly vertically
actually I tried to help it out I gave
it a little little tilt it didn't fall
down so okay I'm gonna I'm gonna add
some rolling resistance I had this nice
rolling resistance thing so it got
really strange it's completely
unphysical so this is not working out I
can't I can't ship that alright so I
scratched my head for a while I talked
to my friend Dirk for quite some time
about this problem and so I had this
idea like well I want to update the
velocity I have all these forces or
these torques I have the gyroscopic
torque I have maybe some external force
I got constraints each one of these
torques is kind of wanting to modify the
angular velocity it's creating a delta
for the angular velocity what if I just
treat each one of those extorts
separately so I'll take this one torque
and I'll use explicit solution for that
or I'll take another one out and I'll
use implicit so that's what I did so I
broke it apart my applied torques I'm
just using explicit but for the the
gyroscopes of torque I'm gonna do
implicit order because I know
implicit order is super stable all right
and then so I get those deltas and then
I add those Delta's in to get the new
angular velocity all right but I set
myself up for a kind of difficult
problem I have this implicit equation to
solve because remember I you apply the
backwards difference and that converts
the difference equation into an
algebraic equation but the algebra
equation is implicit there's this
there's something to be solved here so I
have to solve for Omega 2 but Omega 2 is
it's a period in this cross product so
this is actually a nonlinear algebraic
equation so I have some function of
Omega 2 that's equal 0 that's kind of
like a root finding problem now there's
only a root finding problem it's three
dimensional remember this that I the
inertia tensor is a three by three
matrix Omega is a three vector another
problem I don't have I - I don't have
the next inertia tensor we talked about
this already there's a solution remember
I said in the you compute the inertia
tensor kind of in local coordinates well
you can actually solve this whole
problem in local coordinates and in
local coordinates the inertia tensor is
constant another thing I need to do is I
need to solve a three dimensional
nonlinear equation well for
one-dimensional equations there's this
well-known method called newton-raphson
which is a route-finding technique and
it involves computing the derivative and
dividing by the derivative been
iterating basically until things
converge you can do this also in
multi-dimensions
so instead this the derivative I have
this J which is the Jacobian and that's
basically the function
differentiate with respect to all the
different variables they exist so
anyway so if I work in in local
coordinates I can I can derive that
Jacobian so I did that and here's the
exact formula for that Jacobian using
the inertia tensor in body or local corn
so that's isube
and that involves this this Q matrix
which takes a vector and then returns a
skew symmetric matrix anyway that's the
formula so now I have a gyroscopic
torque solver using an implicit order
and here's the pseudocode for that so
first step convert angular velocity into
local coordinates second step compute
the air in the nonlinear equation I want
to solve then compute the Jacobian and
then so I the Jacobian is a three by
three matrix and I have to basically
invert by the Jacobian or just do a
solve directly and then that updates
that anger velocity and body coordinates
and then I convert the the angular
velocity back out in the world
coordinates and notice there's no kind
of iteration loop I found that just one
newton-raphson iteration was enough so
you could integrate more but I didn't
see any any and you gain from that
alright so here's the capsule
yeah so it falls down that's exactly
what I wanted
here it is with rolling resistance so
yeah you can see you even stops rolling
that's really really nice solution
that's the kind of stuff we like all
right so the gyroscopic torque is not
only important for getting kind of
correct behavior and you know stability
problems there's also some pretty cool
effects that you can pick up when you
simulate the gyroscopic torque correctly
you know we know that when things rotate
they wobble and they move all kinds of
crazy ways things can kinda get even
crazier in space so this is the nipa
cough effect reproduced on Space Station
so look at this object is spinning about
one axis flipping around and spinning
about the same axis so I wanted to see
if I could recreate this effect because
like if I can recreate these kind of
cool wobbly effects you know maybe I'll
have that are looking explosions in my
game all right so here's with explicit
integration so this is the the unstable
solution so you can see in the upper
left I have the kinetic energy and every
time the object flips around its it's
increasing you know it's reproducing the
effect very well quite it's unstable
all right just goes faster and faster
all right if I drop that gyroscopic term
it's very boring not know I don't have a
cool explosion now or whatever in my
game it's just when you drop the
gyroscopic term would you actually end
up doing is you're conserving angular
velocity instead of conserving angular
momentum and then here's the the
implicit solution so you can see the
kinetic energy up there it might be
sorry that it's too small but it's
actually it's actually damping so the
kinetic energy is going down but I think
this is OK for games I'm getting kind of
cool effects and things are damping down
that that's okay all right so that's
about all I have I might go through a
couple conclusions and things you want
to look for when you're looking at
numerical methods so first if you have a
numerical method
maybe you invented a new one right you
want to make sure that as the time step
gets smaller that you're approaching an
exact solution and it's good to have
something like the mass spring system
where you know the exact solution so you
can compare so if your numerical method
always returns 42 no matter what
differential equation you put in it's
it's not converging accuracy so accuracy
is important when you're comparing two
numerical methods for the same time step
is one method closer to the exact
solution than the other that's accuracy
for games we don't need a whole lot
accuracy first order accuracy is fine
and then of course stability we have to
make sure that our simulations don't
blow up so this is where it really pays
to have an understanding of what what
you have what's implemented and what are
its limitations and I encourage you to
test those limitations you know know
know how far you can push your
simulation before before it blows up
all right so here's some references and
yeah you'll be able to get this slides
at box2d org and you can reach me on
Twitter there's lots of references out
there for this stuff sorry if I don't
have everything also when I upload the
slides there's gonna be a bunch of
stuffs I had to skip over I just
couldn't fit it all in so you'll be able
to see that stuff too if you download
the slides and feel free if you want to
go to the next talk now go ahead if you
want to stick around ask a question
please come up to the mic thank you
[Applause]
you