Geometric Algebra, First Course, Episode 15: Dynamics

Transcript

hello and welcome to the stem c studio channel and another video in our series geometric algebra first course and of course we'll be continuing the little mini series that we're doing called flatland where we're doing the physics of a rigid body moving in the plane up till now our rigid body is able to move at a constant speed in a certain direction a constant velocity or it can be at rest which is also a constant velocity or our body can rotate at a constant speed in in the same plane which is a constant angular velocity and of course it can not rotate but that's also a constant angular velocity of 0. what we're going to do in this video is we're going to bring in interactions in the picture that is the body's angular velocity or linear velocity are going to be are going to change that's evidence of an interaction of some sort and those interactions are brought about by forces and talks there'll be a little bit of theory to start with and then we will implement it in code okay let's get started we'll start with the uh the linear momentum and how we update the linear momentum of our rigid body the linear momentum of course will translate into change into the linear velocity by the relationship that the linear momentum has to uh the velocity which is usually that the momentum equals the mass times velocity in non-relativistic physics okay so the the important formula that we need um is this fact which is newton's second law which is that the force is equal to the rate of change momentum or in other words if you apply a force then you will get a change in momentum of what you're applying it to which is p dot is the rate of change in momentum with respect to time so let's take this formula this very kind of compact formula and expand it out get its true meaning the first thing that we'll do is we'll make the time parameter an explicit argument for each of these things so we'll treat them as like functions and the second thing we'll do is we'll expand out the dot definition into dp by dt and then we'll expand out the definition of the derivative with respect to time to be this this difference dt here is uh an infinitesimally small amount of time a piece of time that's so small it's it's smaller than every other finite amount of time and so that makes this equation exact now we rearrange the equation we rearrange this term and we keep this term and we get this okay so this is an exact equation because dt is infinitesimal and then we can go to our update rule by turning this into an approximate when we replace dt by a finite amount of time so the final momentum is approximately equal to the initial momentum plus the force times the time interval so what does that look like in code quite simply the body's momentum the linear momentum uh the final momentum is the initial momentum plus f multiplied by delta t okay let's go and implement that so i'll hide the theory here's our running program and what we want to do in this program is we want to introduce some kind of a force so i'm going to bring it at this point and let's just make the force constant force to the right this is going to update the or cause an update of the momentum of the body [Music] our rule is that the and don't worry about the fact that it flew off our rule is that um the final momentum is equal to the initial momentum plus the force times delta t and you can see the the body is accelerating away there and i'll just clean up a little bit here okay if i make this step a little bit smaller 0.005 for example then things will be slower at the start but this is a constant force so it's getting faster and faster and faster as time goes on so let's stop that and we will bring back our definitions let's now take a look at the angular momentum well the total angular momentum of the system is defined to be the sum over all the different points in the body which are indexed by i of the position vector of that point from the origin wedged with the momentum of that point what we can do now is we can differentiate this expression l dot which is the rate of change of the total angular momentum with respect to time is then just the the sum on this sum is applies to both terms the sum of x dot which with p i plus x i wedge with p i dot okay now the first term is going to vanish because x i dot is parallel to p i and since this is a wedge product and the two vectors are parallel you're gonna get zero so this term is going to vanish and then this one we can replace the p dot i with just the force on that eighth particle the rate of change of the total n momentum is the sum of x i wedged with f i now remember that if we introduce the center of mass we can refer to a particular point in the body by two vectors one is the direct vector the direct position vector from the origin to that particle and then the other is by taking this relative offset from the center of mass so if we plug that in this equation here then we notice that this x part can factor out because it's a constant and then we get x wedged with the sum of f i and the sum of f i is just total force on the body and then for the final part the ri wedged fi we just leave that as it as is so this total angular momentum we talk about it as being composed of two parts this part is what we call the orbital uh the orbital uh contribution to the angular momentum or the rate of change of angular momentum rather and this part is the intrinsic contribution to the rate of change okay so the part that we're interested in is is the part that relates to the rate of change of the angular momentum of the body because this reference refers to the spinning body and we'll just we'll just drop the intrinsic and and call that l dot for now so the final thing is we we define torque um to be uh ri wedged fi um so this is basically the force times the perpendicular distance from the the point of the center of mass and then that gives us the equation that the rate of change of the angular momentum of the body is equal to the applied torque so once again we just go through the motions of expanding the definitions here adding the time arguments and we get that do a little bit of rearrangement of that and this is the exact formula for for an update an infinite decimal update but since we're going to do finite calculations our end angular momentum is equal to the initial angular momentum plus n times delta t the finite change in time so it looks very similar to the the update rule that we actually have for the momentum except the force has been replaced by the torque so here's our update rule the fi the final uh angular momentum is equal to the initial angular momentum plus n delta c and that translates in code to body.l equals dot l plus n times delta t so this is an assignment of the right hand side to the left hand side causing an update let's try that in code okay so recall that n is a bivector now okay as a torque it's a bi-vector because n is the wedge product it's the sum of all of the little torques applied to each particle so let's do n equals let's keep it really simple e1 wedge e2 is the sort of the simplest bi vector we know the time being i'm just going to set our force to be equal to zero and let's do the update i'll move this down a little bit let's do the update to the angular momentum the body's angular momentum is updated [Music] by adding n times delta t okay so let's give that a go so here we go we've actually got a slightly smaller time interval but you can see that our body is starting to rotate if i change the torque and put a minus sign in there then you can see that the body starts to rotate in the other direction all right now i did observe that when i ran this for a little while something interesting happens so i'm just going to do that and then we'll take a look and see how ways we can fix this so everything looks pretty good you know the square's basically staying the same size but it does appear to be getting a little bit bigger i'll just let it go for a little longer i'm going to turn off the documentation so as it's getting times going on it definitely seems to be getting bigger yeah so it's definitely a lot bigger than it was before getting faster and faster rotating all right so what's going on here um what could be going on well let's go and look at our go and look at our square and let's take a look at how the points were updated so you can see that what we're doing is we are applying a transform which is basically converting from the local position to the to the rotated and translated position and the way that we do that is we apply the rotation in the first part and then we apply the translation in the second part the points are getting further away from the origin um in this in this case you know body.x is is not changing it's just zero so what seems to be happening is that r appears to be getting larger and larger in each step and that would kind of cause the point to move further away from the origin so that's my belief and of course the rotor should have magnitude 1 overall so let's um let's kinda like see if we can see that going on so we'll come back to here and uh here's a place where we are we've just updated the uh ah so i'm going to print print line r equals dollar ah okay and that's body.r by the way just to be clear or call it body.r okay okay let's let's run this then now we're obviously going to get a lot of a lot of stuff keeps on increasing it's not obvious that this is getting bigger because it's changing the the amount of uh of each part is changing but what we probably will see here is that at some point um different parts of this will be greater than one there we go okay we can see the different parts greater than one um that shouldn't be the case because each part like the scalar part is supposed to go like cosine theta or cosine theta over two and the uh by vector part here goes like sine theta so that should also not be getting too big and you can see it's getting getting too big here so let's see how we can how we can stop that so what we what we would like to do perhaps is just to keep the magnitude of our uh rigid body close to one so let's write a function that can compute the magnitude of our rotor in fact for that matter let's write a function that can compute the magnitude of any geometric number okay so here we go we'll just kind of like create the prototype we're going to return like i'm going to return a number first now recall that the squared norm which is the square of the magnitude of a multi-vector is defined to be the multi-vector scalar product with the reverse of the multivector so this is what we call the squared norm okay so i'm going to write can't say quad equals that okay but we don't want the square of this multi vector so we need to take the square root of this quaditude if you like i'm going to call it quantitude okay it's not the magnitude it's the square of the magnitude uh we need to take the square root of it and return that so what we need to return is something more like math dot square root and then because the quantity is a scalar then the component that we're interested in is the a part okay so there we go magnitude is declared but its value is never read all right so finally now what we want to do is we want to take this magnitude and we want to kind of like keep the um [Music] keep the magnitude of the rotor to be one so we want to do something like this body dot r equals body dot r divided by the magnitude of the body dot so we'd like to do something like this now unfortunately it's not going to work just at this moment in time because recall that when we do operator overloading division we are only dividing by something that on the right hand side is a geometric so we can fix that by changing this to return not a number but instead we're going to return a new geometric and this will be the scalar part and then the other parts will be zero okay all right so now let's print out the magnitude of the body [Music] and we'll run our simulation okay so there it is it's going round and round it seems to be staying close to one if not exactly sometimes it's exactly one sometimes it's close to one and the body is picking up speed don't seeing it get don't see it getting any bigger in this moment in time let's check down here you should print a few less times certainly seems to be certainly seems to be staying close to one as we'd hoped i didn't think this is going to produce so much output let's uh just hide this now okay and we'll let it go for a second there we go speeding up we could of course reduce our time interval make that smaller and that would be another way of having it can like go faster initially or we could increase the the torque as well but we definitely seem to be uh we seem to have managed the uh the problem i guess one thing i could kind of like do here is like it's just sort of demonstrate that you know you could have different you could have different talks here's one we've made a little bit bigger going the other direction i'm going to set to zero uh we could also simulate say a spring so a spring would have a restoring force proportional to how far you are away so if the force was minus the body dot x and i've left out spring constant at this point and i just got like the minus body dot dot x well nothing's happening because we haven't disturbed it from equilibrium but let's say that we move the body over to 2e1 as our initial position and you can see that the body travels over we should change our time interval maybe and then it goes back or we can uh what we can do is we can increase the strength of the spring so let's say use 10 times the body dot x and what we should get is simulation of the motion of the body under a spring force and so you can put in any force that you want here you can put in your gravitational forces constant forces um electromagnetic forces um there are some kind of like i'll have to explain how that we do a magnetic uh force but it could be done essentially if you if you want to play around with it it's going to be the contraction of the velocity with the magnetic field and the magnetic field will be a bivector so you'll be taking the magnetic you'll be taking the velocity out of the magnetic field it's time to wrap up i hope you enjoyed that take care happy coding bye