Ordinary Differential Equations and Numerical Integration

[Martin Rampersad]

I have been trying to wrap my head around physics, kinematics, and dynamics for quite some time by reading all the papers I can find. However I am a stuck at this idea of different types of numerical integration being necessary for each time step taken.

The common equations I have found for motion are as follows:

Code: Select all

r_final = r_initial + v_initial * timestep + 0.5 * a * timestep * timestep
v_final = v_initial + a * timestep

However, I constantly see references to Euler Integration and Runge-Kutta Methods.

Euler Integration works as follows (I believe):

Code: Select all

r_final = r_initial + v_initial * timestep
v_final = v_initial + a * timestep

Why is it that the exact equation is not used in physics engines, or is something missing in my understanding of what's going on?

[h4tt3n]

Hi Martin,

Welcome to the frustrating world of numerical integration! First, there is no such thing as an exact equation. All numerical integration schemes have an error margin called the order of magnitude. Some integrators are good in some situations and poor in other. Good all-round integrators are the Velocity Verlet 2nd order (low order but symplectic) and Runge-Kutta 4th order (high order but not symplectic).

The Euler equation is the mother of all integration, but actually quite poor. Noone actually uses it.

Code: Select all

(accelerate masses)
r_final = r_initial + v_initial * timestep
v_final = v_initial + a * timestep

Instead use symplectic Euler. Notice how small the difference is.

Code: Select all

(accelerate masses)
v_final = v_initial + a * timestep
r_final = r_initial + v_initial * timestep

Cheers,
Michael

[Martin Rampersad]

Thanks for your reply, Michael.

I understand that different types of numerical integration come with different levels of applicability (capacity to deal with stiffness, etc.), but what I don't understand is why the equation for motion is not used directly in the integrator?

Martin

[Dirk Gregorius]

When you solve a constrained dynamic system you don't have an ordinary differential equation (ODE), but a system of algebraic and differential equations (DAE). This is *not* the same! Actually in the most simple form you need to solve :

dv/dt = M^-1 * ( f_ext + f_c ) = M^-1 * ( f_ext + JT * lambda )
dx/dt = v

C(x) = 0

The constraint forces resulting from the algebraic constraints are *not* know. The classical solution is to differentiate the constraint equations twice, plug them into the equations of motion and solve for the constraint forces (see e.g. Baraff). Once you have the constraint forces you could of course use any integrator you like. This approach has several disadvantages like e.g. infinite friction forces (falling rod example). In literature this is sometimes called solving on the acceleration level.

Physic engines use a different approach and solve on the so called velocity level. In short you basically update the velocities without considering the constraints and then you filter out all velocities violating the constraints in a second pass. This is what a Projected Gauss-Seidel (PGS) or Sequential Impulse (SI) solver do. For a more mathematical formal description see either the PhD of A. Anitescu or the poster presentation of Erin on the GDC.


HTH,
-Dirk

[Martin Rampersad]

Ok, I think the picture is becoming more clear. If you were simulating unconstrained rigid bodies (no interaction, no friction), you could directly use the equation for motion.

Ignoring contact and joints, if you wish to simulate drag (where the force of drag is proportional to the velocity, aka the rate of change of velocity depends on the current velocity), you have an ODE:

So the problem with drag is that it is very difficult to analytically integrate. However, if you use Euler's method (with small enough timesteps) we can numerically integrate the function and generate an approximation.

Code: Select all

a = 0; // start the frame with no acceleration
a += g; // add the force of gravity -9.80
a += C * v; // add the velocity times the drag coefficient (C)
v += a * t; // euler integrate for velocity over time (t)
r += v * t; // euler integrate for the position (r) over time (t)

I can see why Euler's method is not good for large timesteps because while the drag is continuously changing (drag is getting smaller at the same time as velocity is getting smaller), Euler's method treats it as a constant over the timestep. This error is compounded with the fact that the velocity is also continuously changing and is also treated as constant.

My calculus skills are pretty weak, I would understand more clearly if someone could help me derive the general solution for this simple drag scenario I have depicted above. Then I could compare the analytical solution with an Euler approximation for shrinking values of t.

Ah, I see even for this simple example Wikipedia reveals the difficulty of finding a general solution.

After that, it becomes even more complicated (joints and contact) because computing those forces depends on the positions and physical characteristics of the other bodies?

[bone]

Martin-

The simple answer to your initial question is:

Your "equations of motion" are only correct for a constant force during the entire timestep.

Any force that depends on time, position, or velocity will NOT be constant. And that's just about every force that I can think of, with the possible exception of gravity (which for typical Earth-based games can often be assumed to be constant).

Using those equations of motions is often referred to as Parabolic Integration, IIRC. Seeing as how most forces are not constant, then, even that integrator has errors (and significant ones compared to some of the other possible integrators).

Finally, throw in the integration of rotational components and the Parabolic Integrator performs quite poorly (and even for constant torques if you are properly simulating gyroscopic a.k.a. Coriolis effects, but few game engines do).

-Bone