True Physics

One of the vital organs of a physics engine is the code that integrates Newton's laws of motion. There are several techniques which employ approximations of calculus in order to numerically solve this problem. I am proposing a scheme which can solve this problem without approximation.

In this article you will learn

• How the stepping equations are derived from Newton's laws.
• A good method for doing this approximation, which is very common in many third party physics engines.
• The Kinematic Integrator method, which allows you to use unapproximated calculus in the force calculations. This enables extremely large time steps.
  

In this article I will just be talking about the linear motion, not the angular motion. If you want to see a pretty decent way to integrate the angular portion check out the last section of Quaternions: How.

Newton's Laws

Given in short:

1. Without external influences, objects preserve their state of motion. If the object is at rest, it remains so. If it is moving it will continue to move at that same speed.
2. If a force acts on an object, the object will change its state of motion. The amount of change induced by the force is related to the mass of the object. i.e.
   F = ma
3. If one body applies a force on another body, then the other body applies an equal and opposite force on the first body

The third law is taken care of during resolution of collisions, and applications of forces, etc. The first two laws are handled in the integration stage.

The net force acting on a body is directly proportional to acceleration. Acceleration is defined as the time derivative of velocity. In order to get the velocity of an object, given the acceleration, we need to integrate.

Velocity is an important portion of the physical state of an object. Several forces require knowledge of the velocity. However, the most important quantity is the position, since this determines where we place the object in the world. The velocity is defined as the time derivative of position, so in order to get the position we must integrate the velocity.

Since the information we start with is acceleration, and the information we desire is velocity and position, then we must integrate the acceleration two times. For this reason the equations of motion are usually second order equations. In order to fully solve a second order equation, we need to know two quantities. These two quantities are usually the initial position, and the initial velocity.

How do we integrate?

When it comes to physics engines I always start my integration by taking a derivative. This may seem anti-intuitive, since an integral seems to be the opposite of a derivative. But with physics, the thing we are integrating can be expressed in terms of derivatives, and nothing is as easy to integrate as a derivative.

The equations of motion can be expressed in terms of two differential equations

a = ∂v/∂t

v = ∂x/∂t

In order to see how we solve these equations numerically, first consider the standard mechanism for taking a derivative.

∂f/∂t = (f₁ - f₀) / (t₁ - t₀) ; t₁ → t₀

Which defines the value of the derivative at t₀.

On a computer, we have finite precision arithmetic, and so we cannot allow t₁ to get arbitrarily close to t₀. Therefore, we can only approximate derivatives as finite differences.

∂f/∂t ≈ (f₁ - f₀) / (t₁ - t₀) ; t₁ = t₀ + Δt

Here Δt represents our time step. Obviously, the smaller the time step is, the closer we get to representing a true derivative. Making this simple approximation to the representation of the calculus has introduced an ambiguity. Is this the value of the derivative at t₀ or t₁, or at some point in between? If we choose t₀, then the approximation is called a forward, or explicit derivative. If we choose t₁, then the approximation is called a reverse, or implicit derivative. The quality of the representation of the equations of motion ends up being tied to our choice.

If we choose to use only explicit derivatives to represent the equations of motion, then we have

a₀ = (v₁ - v₀) / Δt

v₀ = (x₁ - x₀) / Δt

In these expressions, we know we are using a forward derivative because the a and the v on the left hand sides are evaluated at t₀. We can rewrite these equations, so that we can find the values of x and v at time t₁, given that we know what the values are at time t₀.

v₁ = v₀ + a₀Δt

x₁ = x₀ + v₀Δt

These equations are called the explicit Euler stepping equations. We see that if the value of a is zero, then the equations will maintain the motion of an object, satisfying the first law. The solitary occurrence of a is the contribution from the second law, given that you calculate a as the net Force divided by mass.

As it turns out, the explicit Euler stepping equations are very crappy. The acceleration ends up introducing energy into the system, causing it to be very volatile. This volatility can be managed by making the time step smaller, but this increases the computational expense.

The current gold standard for integration in high profile physics engines is the semi-implicit, or symplectic, Euler stepping equations. These equations represent the acceleration with an explicit derivative, but the velocity is represented with an implicit derivative.

a₀ = (v₁ - v₀) / Δt

v₁ = (x₁ - x₀) / Δt

Since on the left hand side, the acceleration is evaluated at t₀, it is a forward or explicit derivative. Since the velocity on the left hand side is evaluated at t₁ it is a reverse, or explicit derivative. The combination of these two gives us the name semi-implicit. These equations can be rearranged so that x and v at t₁ are functions of their values at t₀

v₁ = v₀ + a₀Δt

x₁ = x₀ + v₁Δt

These equations are solvable, so long as we evaluate the velocity equation first, so that we have the result to plug into the position equation.

So how accurate are these stepping equations?

The semi-implict Euler technique is a second order method, which means that the error is tied to the second power of the time step. If you cut the time step in half, you quadruple the accuracy of the method.

Also, the semi-implicit Euler method obeys a different property, which makes it useful. It conserves energy on average. What this means is that, although the error you accumulate over the course of the time step might increase or decrease the effective energy of the system, over time the energy will average out to a constant. This property is where the symplectic term comes from. I'm not going to explain what symplectic means, but if you are interested, search for Hamiltonian dynamics, the equations of which are used to prove that the semi-implicit Euler method conserves energy.

We see in this instance that there is a defining difference between accuracy and stability. The semi-implicit Euler method may not be the most accurate method ever, but it is stable since it conserves energy. On the user end, we don't tend to see the tiny inaccuracies, but we do tend to notice if the physics wigs out and everything explodes. Therefore, we would like our inaccuracies to NOT lead to instabilities.

Can you have perfect integration?

With numerical integration schemes, such as Euler based techniques, the only way we can get perfect integration is if we let the time step become zero, which we know we can't do. However, it is possible to achieve perfect integration in a computational simulation, and I will tell you how, but we will need to make a slight departure from the standard development of numerical integration.

Let's start with a very simple example of perfect integration.

Consider that you are writing the game Lunar Lander where you must use thrusters on a spacecraft to guide it safely to the landing pad. In this game, the only forces that come into play are the gravitational force, and the thrusters on the spacecraft.

Over the course of a single time step, we know which thrusters are on. We consider that the force due to a thruster is constant over the course of the time step. We also use a constant force for the gravitational force. Thus we have a full knowledge of the value of the force at all intermediate points during the simulation time step.

Add all of the forces together to get the total force. Dividing the total force by the mass will give you the acceleration.

a = (F_grav + F_thruster + ...) / m

Now, we need not resort to any approximation of calculus in order to solve this problem. Since we know the form of the forces (constant) at all points during the simulation time step, we can solve the problem directly - analytically. The true solution and stepping equations without approximation are given by

v₁ = v₀ + aΔt

x₁ = x₀ + v₀Δt + (1/2)aΔt²

Anyone who has taken a physics class, even a beginning physics class, should recognize these equations. Restate: if you have taken a beginning physics classes you have seen these equations before, even if you don't recognize them. They are called the kinematic equations of uniform acceleration, and are the primary feature of almost all beginning physics problems, especially problems involving a baseball, basketball, or any other type of sports related sphere.

Here the a is constant, so we don't need to use a subscript to specify the point in time where we are evaluating it. These equations look almost like the stepping equations we were using before, except the position equation has an extra term involving acceleration. These equations are perfect, which means that there is no upper bound of accuracy that is related to the size of the time step. In a real application the upper bound on the size of the time step would not be related to the accuracy of the integration, but rather to the ability to perform collision with the landscape. Collision issues aside however, this method will give you the right answer, no matter how big the time step is.

If these stepping equations are perfect, why don't we use them? The answer is that these equations are only perfect if the forces are really constant. If we use these equations to simulate forces that are not constant, e.g. the spring force, we will not get very good results. If we decide to use these equations to step the physics of our system forward in time, we have essentially moved the approximation out of the calculus and into the representation of the forces. And the approximation applied to the forces ends up having worse results than applying the approximation to the calculus.

In short: the kinematic equations rock for forces that are constant, but really suck for forces that aren't.

If we can't use the kinematic equations all of the time, should we ditch them all together? Is it possible for constant forces to use the kinematic equations, and non-constant forces to use the semi-implicit Euler equations? Then our integration would be perfect with respect to the constant forces, at least. Under such a scheme, how would we go about mixing constant and non-constant forces?

Consider once again the semi-implicit Euler equation. We can plug the velocity equation into the position equation, since they both involve v₁. The result is similar in appearance to the kinematic equations

v₁ = v₀ + a₀Δt

x₁ = x₀ + v₀Δt + a₀Δt²

We now see a very striking resemblance to the kinematic equations. We actually have three terms that match exactly. The first term in the velocity equation, and the first two terms in the position equation. Only the terms that involve acceleration are different. It can be shown in a general case, that the three matching terms will always match, regardless of the form of acceleration. The form of the acceleration only affects the terms which involve acceleration.

We now have a clue as to how we might combine a constant force with a non-constant force.

v₁ = v₀ + ( a_CΔt + a_NCΔt )

x₁ = x₀ + v₀Δt + ( (1/2) a_CΔt² + a_NCΔt² )

Thus, the constant force gets the benefit of using the kinematic equations, and the non-constant force can use the semi-implicit Euler technique. You may see some simplification that can take place, since the acceleration terms contain common factors. Now, it is possible to analytically solve for the acceleration terms given other forms of acceleration. These terms may not have common factors, so I chose to keep the acceleration terms separated.

The Kinematic Integrator:

The kinematic integration method is given as follows

v₁ = v₀ + dv

x₁ = x₀ + v₀Δt + dx

Which will give us the new position and velocity, given we supply the value of the previous position and velocity as well as the two integral parameters dx and dv. So what are dx and dv? They are generalizations of the terms which involve integrals of acceleration.

In Euler based integration techniques, you usually have a unique method to calculate each different kind of force. There might be a method for a spring force, and a method for a pulse, etc.. This resulting force is divided by the mass, and accumulated with whatever existing acceleration is already acting on the object.

Using the kinematic integration technique, the entire point of the force calculation method is not to calculate the force, but to calculate dv and dx. These are, likewise, accumulated with whatever existing dv's and dx's that were previously acting on an object.

The point of this generalization is to move the calculus out of the stepping equations, and into the force calculation, so that each different kind of force can give the true result of the acceleration integrals. This leaves the stepping equations in a form that is not dependent on the form of the input forces.

A method which calculates a constant force will return the exact dv and dx of a constant force. A method which calculates a spring force will return the exact dv and dx of a spring force. This is really great, but what do we do if we don't know how to calculate the true integral of a given force? Simple, we just go back to using the semi-implicit Euler technique. We have already expanded this technique out into kinematic form, so we know what dv and dx are.

dv = a₀Δt

dx = a₀Δt² = dvΔt

The spring force

So how do we calculate the dv and dx of a spring? As an example, consider a spring which connects a mass to the origin. The spring has a spring constant k, which determines how tight it is. The force acting on the mass has the form

F = -kx = -mω²x

Here ω is the angular frequency of the spring, and the reason for defining it is because the form of motion of the spring is given by a sinusoidal function of the following form

x = A sin(ωt) + B cos(ωt)

We can take the derivative of this expression to find the functional form of the velocity. Then we can use x₀ and v₀ as our initial conditions in order to fix the values of A and B.

v₁ = v₀ cos(ωΔt) - x₀ω sin(ωΔt)

x₁ = v₀/ω sin(ωΔt) + x₀ cos(ωΔt)

These are the true, un-approximated results of the integration. We can use these to find the values of dv and dx, merely by equating these results with the left hand side of the kinematic integrator equations. We find that

dv = v₀ ( cos(ωΔt) - 1 ) - x₀ ( ω sin(ωΔt) )

dx = v₀ ( sin(ωΔt)/ω - Δt ) + x₀ ( cos(ωΔt) - 1 )

If the parameters of your spring do not change, and you are using a fixed time step, then the coefficients of x₀ and v₀ do not change, and should only need to be calculated once.

Again, the accuracy of the stepping equations does NOT depend on the size of your time step, because we are using the exact result. I have compared the numerical result of the kinematic integrator with the exact result in the case of a spring, and there is literally zero error - no matter what the time step is, or how long the system runs.

Multiple forces:

Now, there is a problem that arises when you try to apply two forces that depend on position or velocity. They both are integrated without a knowledge of the other. Therefore they are oblivious to the changes in intermediate position or velocity caused by the other force. These changes would have altered the final result, had they been taken into account.

A good example of this is the example of having two springs.

If one spring acts alone, the result of the kinematic integrator is exact, and so the result is always correct, regardless of the size of the time step. The story is different when two springs act simultaneously.

The system begins to accumulate error, due to the misrepresentation of the forces. However, the error is very small, almost insignificant. In a test of such a system, I calculated that the error of the kinematic integrator was 12,000 times smaller than the error from the semi-implicit Euler method. Unlike the semi-implicit Euler method, however, the kinematic integrator is not guaranteed to preserve energy. Therefore these very small inaccuracies build up, and over a very long time, the system diverges. Now, the introduction of even the smallest amount of friction to the system would conceal these tiny inaccuracies, however, there ends up being a better way.

The error due to multiple forces is actually an error in the calculus. We can play the trick where we move the error out of the calculus and into the representation of the forces, by making an approximation of the force. The approximation I have in mind ends up introducing less error than the error due to multiple forces. Also, this approximation tends to remove energy from the system, which makes the system stable. With this method, you would still calculate dv exactly, like before, but when you calculate dx, you do it like this

dx = (1/2) dvΔt

Using this method to calculate dx actually beats out calculating it exactly in cases where multiple forces that depend on position are in play. This method is called the average acceleration method. On a test performed on a system with two springs, the average acceleration method removed about 5% of the energy from the system. After about 3 hours. In the test I ran of this sample system, the error introduced by the average acceleration method was 58 million times smaller than the error introduced by the semi-implicit Euler method.

It is recommended that forces which depend on position or velocity use the average acceleration method of calculating dx.

Conclusion:

Using the kinematic integrator allows you to remove the process of integration from the stepping equations, and push it down into the force calculation. This enables you to taylor the integration method to each different type of force. Having this freedom allows you to provide exact solutions, if they are attainable. If the solutions are not attainable, you can still use the semi-implicit Euler method as a standby.

This method introduces a trivial amount of extra computation, i.e. two or three operations per force applied. In the spectrum of numerical methods, this method has roughly the same computational expense as a first order method, and acheives results that far exceed the more expensive fourth order methods.

Obviously this method only has utility if you know the form of the forces acting during a given simulation step, and also if the forces are integrable. Happily, most if not all of the forces that are used in run-time simulation are integrable.

With some systems it is possible to completely eliminate error. In others, the error is a very small fraction to the error of other methods that have comparable computational expense. Thus, with the kinematic integrator you can use much larger time steps, and expect higher accuracy from integration.