Numerically Integrating Equation of Motion for an Object

[The_Engineer]

I'm trying to integrate the equations of motion for a object.

F + mg = ma

where F is the drag force, g gravity, a is acceleration, etc...

I'm trying to do it numerically and I'm confused about one thing:

Since this is a 2nd order vector differential equation, should it be equivalent to a system of two 1st order ODEs or four 1st order ODEs?

 

[AlephZero]

It is a system of two 1st order vector equations. if you call the unknown variables x and v (position and velocity vectors), a = dv/dt, so your equation becomes

dv/dt = F/m + g

and the other equation is just

dx/dt = v

Iif the motion is in two dimensions, the vectors x and v will each have two components, and you need to resolve F and g into their components as well, so in a sense you end up with four 1st order scalar equations.

 

[Filip Larsen]

Welcome to PF!

You will get one (simple) ODE for each dimension in your position state, and one ODE for each dimension in the velocity state. So, for a 2D problem you'd have 4 scalar ODE's, and in 3D you'd have 6.

 

[The_Engineer]

It is a system of two 1st order vector equations. if you call the unknown variables x and v (position and velocity vectors), a = dv/dt, so your equation becomes

dv/dt = F/m + g

and the other equation is just

dx/dt = v

Iif the motion is in two dimensions, the vectors x and v will each have two components, and you need to resolve F and g into their components as well, so in a sense you end up with four 1st order scalar equations.

Thank you that was very informative. If I would like to include an initial angle (from where the projectile is taking off) in my simulation, should I convert to polar coordinates?

 

[The_Engineer]

Thank you that was very informative. If I would like to include an initial angle (from where the projectile is taking off) in my simulation, should I convert to polar coordinates?

I just realized that the initial velocity vector is the same thing as the initial angle. Nevermind and thank you!

 

[D H]

Since this is a 2nd order vector differential equation, should it be equivalent to a system of two 1st order ODEs or four 1st order ODEs?

In theory, the answer is yes. In practice, when using numerical techniques, the answer is a resounding no.

Here's the theory. Create a phase vector that contains the position and velocity vectors. This new vector contains 2N elements, where N is the dimensionality of the space. For your 2D problem, this combined position+velocity vector is $\vec u = [x,y,\dot x, \dot y]$. The derivative of this 4-vector is another 4-vector: $\dot{\vec u} = [\dot x, \dot y, \ddot x, \ddot y]$. You already know the first two elements, and you also know the next two via $F=ma$ (presumably you know the net force).

You have a 4-vector and you a way to calculate its derivative, so just use your favorite numerical first order ODE technique to integrate. Easy!

Not so easy. The problem is that this approach ignores the geometry of the problem. If you can find a 2nd order ODE solver that is of the same degree as that 1st order ODE solver, the 2nd order technique is almost always going to be better because it pays at least some attention to the geometry.

A simple example: The simplest 1st order ODE technique is the explicit Euler method for 1st order ODEs. A much improved (but still lousy) technique that pays attention to the 2nd order nature of F=ma is the symplectic Euler method. Explicit Euler doesn't conserve any of the conserved quantities (energy, angular momentum, and linear momentum). Symplectic Euler comes much closer to doing so.

Don't use either of the Euler methods. Symplectic Euler is lousy; explicit Euler is extremely lousy. You still need to understand the Euler methods because all but the most advanced techniques are variations on the Euler methods. A slight step up in complexity is the leapfrog / verlet family of integrators. Done right, these require only one computation of acceleration per time step, but the accuracy is vastly improved. There is no 1st order ODE equivalent to any of the leapfrog / verlet integrators.

You can do even better than this, but now things get complex. There is a problem with higher degree geometric integrators: They get hairy.

 

[K^2]

Let it make very simple for you. For most practical purposes, integrating EoM should be done using Velocity Verlet algorithm. Why? Because it's very simple, and it gives you second order method which happens to be extremely stable for most practically encountered interactions. In particular, it does extremely well with Hook's potential, which is what you'd typically use for collisions. (This method belongs to the family of methods D H mentions. I'm not disagreeing with him, just saying, go for this particular method. Velocity Verlet will make your life easier.)

Here is one notable exception to keep in mind. Verlet, as well as any other explicit Runge-Kutta method, is really pants at doing gravity problems. I don't mean F + mg kind of gravity. That you can do with Verlet. I'm talking plotting space-ship trajectory through a solar system. If you ever need to do that level of simulation, Verlet will betray you. For problems involving gravity, when you need good precision, you have to use implicit methods. 2nd or 3rd order Gauss-Legendre should work much better for these kinds of problems, but that's a whole different level of complexity.