## A ball with air friction (part 3)

In the previous parts of this series, I gave a gentle introduction to simulating a ball with air friction using the traditional Euler integration approach. It’s fast and it works well enough for games, but I have an interest in seeing how to make a ball simulation that does not suffer from the limitations of stepwise integration schemes, and perhaps work out formulae to help with AI, such as working out the collision time with the ground, or equally useful, the velocity required to send a ball to a particular point in space at a precise time. If you consider a game such as Worms, the AI needs to be able to work out the correct angle to fire a weapon to destroy a target, even under the influence of wind and so on. In my early games I could not use mathematics to solve these problems, at least not pure mathematics. Instead I relied on finely tuned engineered algorithms, fiddling with them until they worked well enough. Now of course, I seek to do better than that.

So this takes us to the problem of how to create a formula for the height of a ball under gravity and with air friction in order to be able to work out the height from the last moment the ball was set in motion.

We saw how this is done  in the case without air friction, but now we have an inconvenience. The function for acceleration is no longer constant, and in fact depends on the velocity. This means that we have a circular dependency: the acceleration depends on the velocity which depends on the acceleration:

$\mathrm{a}\left(t\right) :=g+\frac{\mathrm{v}\left(t\right)\,k}{m}$

To those who’s eyes are glazing over, it really is not as bad as it seems. The first term, g, is simply the constant of acceleration. The second term is the velocity times a constant k divided by the mass of the object, m. From the most basic physics, you should remember that force = mass x acceleration. This means that acceleration is equal force divided by mass, and that is what we see: k is a constant representing how strong the effect of friction is. If k is 0, there is no air friction, otherwise the force due to friction is proportional to the velocity, in accordance to k (which should be a negative value, to ensure that friction slows the ball down instead of speeding it up). So this formula should be quite easy to understand. We have to know what v(t) is however, and to find it we need to integrate a(t), which we can’t do without knowing v(t), so we have a catch 22 situation.

To resolve this, we need to look at differential equations. These are equations that relate functions that are derivatives in order to find a solution (which is usually a function, rather than a value). It’s just another form of algebra, with variables representing functions instead of numbers.

What we know is that there is a function a(t) that is the derivative of v(t). So we can say that if f(t) = v(t), then f'(t)  = a(t). The f’ is pronounced ‘f prime’ and just means an equation that is the differential of f. If you integrate f’ you get f. if you differentiate f you get f’. So, for fun, we can try to make an equation that specifies f'(t) and f(t) instead of a(t) and v(t), and hope that some mathematician can work out what the function f(t) would have to be  in order to make that equation valid. Something like this:

$\mathrm{f'}\left(t\right) :=g+\frac{\mathrm{f}\left(t\right)\,k}{m}$

So, there we go. All we need to do is find some kind of function for f that would fit. At this point, I don’t attempt to do this myself. I am not a mathematician, I am a video game programmer. So I outsource that work to a tool in fact.

### Symbolic mathematics packages

These things are a dream come true. No longer do you have to spend hours solving equations by hand, you can just ask any of a number of software packages to do this for you. There are three that I use commonly these days: Mathcad (commercial), Maxima (open source) and the Wolfram Alpha website. Here I am going to start with Maxima.

First of all I need to enter the differential equation into Maxima, for which I have to type it in in a kind of C style manner, with * for multiply and so on:

d1: 'diff(f(t),t) = g + f(t)*k/m;
$\frac{d}{d\,t}\,\mathrm{f}\left( t\right) =\frac{k\,\mathrm{f}\left( t\right) }{m}+g$

Now we can use the Maxima function desolve which solves differential equations for you (if it can find the answer, there isn’t always a solution!):

desolve([d1],[f(t)]);
$\mathrm{f}\left( t\right) =\frac{\left( g\,{m}^{2}+\mathrm{f}\left( 0\right) \,k\,m\right) \,{e}^{\frac{k\,t}{m}}}{k\,m}-\frac{g\,m}{k}$

The above formula is a little bit complicated, but it works. For those of you who are not familiar with the exponential function, it is a magic number, e, to the power of x. Maxima shows e as %e. The expression involving e should be read as e to the power of (k times t over m). When coding this you can use the exp() function.

Having found a formula for velocity at time t, we need to integrate that to give us a formula for position at time t, and differentiate it to get the acceleration at time t. One thing to note is f(0). This is the velocity at time 0 (or initial velocity), so we can replace it with v0. Thus:

$\mathrm{v}\left( t\right):=\frac{\left( g\,m+v0\,k\right)\,{e}^{\frac{k\,t}{m}}-g\,m}{k}$

We can use the Maxima diff and integrate functions to find the derivative and integral of our function, giving us all three formulae for acceleration, velocity and position:

a(t):=diff(((g*m+v0*k)*%e^((k*t)/m)-g*m)/k,t,1);

$\mathrm{a}\left(t\right):=\mathrm{diff}\left(\frac{\left(g\,m+v0\,k\right)\,{e}^{\frac{k\,t}{m}}-g\,m}{k},t,1\right)$

p(t):=integrate(((g*m+v0*k)*%e^((k*t)/m)-g*m)/k,t,1);

$\mathrm{p}\left( t\right):=\int_{1}^{false}\frac{\left(g\,m+v0\,k\right)\,{e}^{\frac{k\,t}{m}}-g\,m}{k}dt$

This gives us the following results:

a(t) = $\frac{{e}^{\frac{k\,t}{m}}\,\left(k\,v0+g\,m\right)}{m}$

pp(t) = $\frac{\frac{m\,{e}^{\frac{k\,t}{m}}\,\left(k\,v0+g\,m\right)}{k}-g\,m\,t}{k}$

On testing the pp(t), I found that it was basically correct, but with an offset (i.e. not giving zero at zero t). When integrating velocity, we need to add a term for the initial position (p0). It should be true that our function p(t) return zero for p(0). Other than this, the curve was correct, so to fix this I simply subtract p(0) away from the result. This is why I have named it pp instead of just p as an intermediate step. The final formula for p(t) is as follows:

p(t) = pp(t) - pp(0)
p(t) = $\frac{\frac{m\,{e}^{\frac{k\,t}{m}}\,\left(k\,v0+g\,m\right)}{k}-g\,m\,t}{k}-\frac{m\,\left(k\,v0+g\,m\right)}{{k}^{2}}+p0$
So, our final three functions, all together, are:
$\mathrm{v}\left( t\right):=\frac{\left( g\,m+v0\,k\right)\,{e}^{\frac{k\,t}{m}}-g\,m}{k}$
a(t) = $\frac{{e}^{\frac{k\,t}{m}}\,\left(k\,v0+g\,m\right)}{m}$
p(t) = $\frac{\frac{m\,{e}^{\frac{k\,t}{m}}\,\left(k\,v0+g\,m\right)}{k}-g\,m\,t}{k}-\frac{m\,\left(k\,v0+g\,m\right)}{{k}^{2}}+p0$

The above three functions are what we have been looking for: functions that give us the acceleration, velocity and position of our ball (vertically) for any time t, given a starting position (p0), a starting velocity (v0), gravity (g), mass of the object (m) and a coefficient of friction (k). Phew! Here are some graphs to show it off:

acceleration at time t

Velocity at time t

Position (height) at time t

And that’s it for part 3. I hope you are hanging in there, because in the final part we will look at how to solve the position equation p(t). In other words, given the initial conditions p0,v0 and constants g,k,m, can we find a formula that will return t, the time to impact of the ball with the ground? In the case plotted above, this function would return around 15. It turns out it can be done, but the answer requires a lesser known function that is not generally provided in math libraries: the Lambert W function.