Unity impulses

As I covered in my previous blog post, the Unity physics simulation falls somewhat short when simulating classical mechanics. Having analysed the problem, it seems that the error is bounded. Unfortunately, this error is fairly significant for distances and forces in the typical range for a video game.

A typical use for game physics is to simulate a jump. The average human jump height is about 40cm and the current world record for a standing jump is 170cm. Calculating the velocities required using the standard equation and running the simulation, the jump heights fall short by 2.8cm for a 40cm jump and 5.7cm for a 170cm jump.

A small part of that error (less than 1mm) is due to the actual peak of the jump not landing on an exact physics frame. According to Unity, the majority of the error comes from PhysX and its use of the symplectic Euler method to perform physics updates.

You might think from the above graph, that the error is increasing as the size of the jump increases. And that is true in terms of the absolute error. However, the purple line, representing the peak height reached at each jump height, is actually slowly curving towards a line parallel to the target dotted line. So as you go out to 100m, the errors become relatively less significant.

Conversely, as you go down to smaller distances, the shortfall becomes so large that attempting a jump of less than 2mm you would not even get off the ground.

There go my plans for Insect Adventure.

As symplectic as possible

The symplectic Euler method applied to physics updates goes something like this:

void PhysicsUpdate()
{
	velocity = velocity + acceleration * dt;
	position = position + velocity * dt;
}

The order of the operations is important here. The position is updated with the already modified velocity.

To see what effect this has on the displacement of objects over time, we’ll have to find a function which applies these modifications for each physics update. Mathematically, the displacement after n updates is:

s(n) = \sum_{k = 1}^n (v(k) \cdot \Delta{t})

Where v(n) is the velocity after n updates and can be calculated with:

v(n) = u + a \cdot \Delta{t} \cdot n

So if we expand v(n) into the sum we get:

s(n) = \sum_{k = 1}^n ((u + a \cdot \Delta{t} \cdot k) \cdot \Delta{t})

To resolve an equation like this we can use the rules for sums. For example, we can pull out the delta-t that multiplies each step:

s(n) = \Delta{t} \cdot \sum_{k=1}^n (u + a \cdot \Delta{t} \cdot k)

And where each step is itself a sum, we can split the result into two simpler sums:

s(n) = \Delta{t} \cdot \left(\sum_{k=1}^n u + \sum_{k=1}^n (a \cdot \Delta{t} \cdot k)\right)

Again, we can pull out the common multipliers:

s(n) = \Delta{t} \cdot \left(\sum_{k=1}^n u + a \cdot \Delta{t} \cdot\sum_{k=1}^n k\right)

In the first sum, u is a constant, so the sum of n instances of u is simply u times n:

s(n) = \Delta{t} \cdot \left(u \cdot n + a \cdot \Delta{t} \cdot\sum_{k=1}^n k\right)

The sum of the numbers between 1 and n is trickier, but it has a known solution.

Naughty Gauss

There is a story that the famous mathematician Gauss was notoriously badly behaved in his mathematics classes. To keep him quiet, his teacher set him the task of adding all the numbers between 1 and 100. But to his teacher’s surprise Gauss came back quickly with the correct answer of 5,050.

His trick was to pair up the numbers from opposite ends of the list:

\begin{align*}
1+100=101\\
2+99=101\\
3+98=101\\
4+97=101\\
\text{etc.}\\
\end{align*}

Notice that if we do this all the way to 100, we are counting each pair twice. So the total answer is half of 100 times 101.

In the general case, for any sum from 1 to n, the sum of each pair will be n + 1, and so the total will be half of n times (n + 1).

Back to the math

We can therefore substitute the last sum term in our earlier equation:

s(n) = \Delta{t} \cdot \left(u \cdot n + a \cdot \Delta{t} \cdot \frac{1}{2} \cdot n \cdot (n+1)\right)

Multiply through by delta-t:

s(n) = u \cdot n \cdot \Delta{t} + \frac{1}{2} \cdot a \cdot n \cdot (n + 1) \cdot \Delta{t}^2 

Multiply through by (n + 1):

s(n) = u \cdot n \cdot \Delta{t} + \frac{1}{2} \cdot a \cdot n^2 \cdot \Delta{t}^2 + \frac{1}{2} \cdot a \cdot n \cdot \Delta{t}^2

Gather the n and delta-t terms:

s(n) = u \cdot (n \cdot \Delta{t}) + \frac{1}{2} \cdot a \cdot (n \cdot \Delta{t})^2 + \frac{1}{2} \cdot a \cdot \Delta{t} \cdot (n \cdot \Delta{t})

And look at s(1), that is, the first physics update frame:

s(1) = u \cdot \Delta{t} + \frac{1}{2} \cdot a \cdot \Delta{t}^2 + \frac{1}{2} \cdot a \cdot \Delta{t} \cdot \Delta{t}

It collapses to:

s(1) = u \cdot \Delta{t} + a \cdot \Delta{t}^2

I think this is what the PhysX developer, that Unity spoke to, meant by “PhysX uses a semi-implicit Euler integrator which leads to omission of the 0.5 factor.” Of course, this only holds for one updated frame.

In the more general case, the expected displacement e for any frame n is:

e(n) = u \cdot (n \cdot \Delta{t}) + \frac{1}{2} \cdot a \cdot (n \cdot \Delta{t})^2

So if we subtract the symplectic version s from this, we can see that the shortfall after n frames is:

e(n) - s(n) = -\frac{1}{2} \cdot a \cdot (n \cdot \Delta{t}) \cdot \Delta{t}

In this equation, n times delta-t is just another way to say the time t at frame n. If we write it in terms of time t instead, we get:

e(t) - s(t) = - \frac{1}{2} \cdot a \cdot t \cdot \Delta{t}

If we pull out the variable t, we can see that the shortfall is always some constant amount multiplied by t:

e(t) - s(t) = \left(- \frac{1}{2} \cdot a \cdot \Delta{t} \right) \cdot t

And a constant amount times t is just another way of saying a velocity. So what this is telling us, is that our velocity is short by this amount. Which implies we can fix it, by simply adding it to our initial velocity.

Therefore, the error isn’t in the calculation of the physics update, but in the application of the impulse.

A little extra push

I concluded my previous post with a lament that there was no accurate way to calculate the force required for a given jump height in Unity. After some time, thought, and experimentation, I conclude that is not entirely accurate.

The force required is accurately calculated by the method I described in the previous post. But to apply that force, you cannot just add it using the impulse force mode in Unity. If there is a bug in the Unity / PhysX physics system, it is at this point.

The corrected jump function is therefore:

if (Input.GetButton("Jump"))
{
	float g = Physics.gravity.y;
	float dt = Time.fixedDeltaTime;
	float bodge = -0.5f * g * dt;
	float v = Mathf.Sqrt(-2.0f * g * height);
	float force = body.mass * (v + bodge);
	body.AddForce(Vector3.up * force, ForceMode.Impulse);
}