Ball with air friction (part 4)

As we saw in the previous installment, we were able to get a formula that provides the future or past height of the ball, i.e. the height at time t. Now we will look at how to solve this equation to find the time when it will impact the ground. First of all, before going through the mathematics, I should point out that there is a reasonable way of doing this that does not require trying to solve the equation. What we can do is search for the answer, by trying different values until we find a value that returns zero. Let’s look at the example curve of the p(t) function again:

Position (height) at time t

Clearly we cannot just search all possible values of t, that would be extremely inefficient. However we can do a binary search, by starting with a point where the ball is above ground and a point where it is below, then then in a manner similar to the “think of a number” game, try half way and keep splitting until we get the answer. For example, if we started with 0 and 20, the first guess would  be 10. This result would be high, so we would split the difference between 10 and 20 and try 15 for the next guess. This result would be low, so we would split again between 10 and 15 and try 12.5… and so on. This works very well and converges to a solution very quickly. For a lot of purposes this is sufficient, and it is a good fallback. However, I am trying to see if we can avoid iterating (that is, looping) at all and just have an equation. This would potentially be faster.

So onward. We want to find the solution to this equation:


This look a bit overwhelming, but that’s ok, we have software to solve it (that is rearrange the equation so that it means the same thing, but states t = <something>). Let’s ask Maxima:


This is asking Maxima to solve our equation for t. Great, let’s take a look at the answer:


Maxima has let us down. It seems to have given an answer, but it gives the answer t in terms of an expression containing t: clearly we can’t use it unless t does not appear on the right (it’s no use having an answer for t that depends on knowing t!). But let’s not give up just yet. If we look at the answer, we can see that it is possible to reduce the equation to a much easier to manage one, with constants made up of the v0,p0,m,g and k values:





The values of c1,c2 and c3 are all generated from the constants v0,p0,m,g and k, leaving a much simpler equation to solve. We should not expect Maxima to be able to solve it though; if it could it surely would. But maybe we can try a different package that may have the answer.

The reason why the equation is hard to solve is because the t on the right of the equation is the power of a constant (c3). This actually turns up a lot in applied mathematics, but the consideration of this problem goes back hundreds of years to Lambert, who’s work was credited by Euler in studies he made of it. Lambert proposed equations involving the variable as a power, (see Lambert’s Transcendental Equation). A result of this was a function that is useful in solving such equations, in recent times named “Lambert’s W function”. It is amazing to see that in my little journey into solving such an apparently simple problem, I arrive at Lambert’s function. It is even more amazing that the true importance of the function was not appreciated until the 1990’s when it was used in quantum mechanics. According to Wikipedia “developers of the Maple Computer algebra system made a library search to find that this function was in fact ubiquitous to nature”. Pretty amazing then that it was not actually appreciated properly until the 1990’s, and explains why it almost certainly does not exist in your favorite math library. Sometimes I wonder if the old mathematicians were cleverer than we are today. They figured out all this stuff with nothing more than pen and paper. We would seem to be decadent in our nature in comparison, with all our computer power at our disposal. Anyway I digress.

I don’t have Maple, but I do have access to Mathcad 15 and Wolfram Alpha. Both give answers to t=c1\,{c3}^{t}+c2, but Wolfram’s is prettier:

Solution to t = c3^t*c1+c2

So, this means that since we can work out c1,c2 and c3 from our constants p0,v0,g,k and m, we have our answer, in one step! Yay! Except of course… no one implements the Lambert W function in math libraries yet. It’s a bit of a shame, because it is not different from other functions we take for granted, such as sin, cos, exp, log. It would be really nice if this function was implemented in floating point hardware, but it isn’t. So, the only way is to create our own implementation, and that requires some iteration. However the nice thing about having the Lamber W function is that it will help with solving all kinds of equations that are unsolvable otherwise; and this iteration to work out the W function is likely to be more efficient than the binary search (although there may not be a lot in it… time will tell).

I searched for an implementation of the W function, which took some digging, and eventually found Keith Briggs’ site: You will find there two equations for the Lambert W function, for what are called ‘branches’. In the above example, we are using the ‘principal branch’. There are only two branches in the real domain (i.e. not using imaginary numbers). Each branch gives a solution: remember that even in this case with air friction there is a theoretical solution backwards in time. For the case we are studying here, only the principal branch is needed (often written as W with a 0 subscript; the other branch is signified with a -1 subscript).

All well and good, but does it work? Here is the expanded out final function from Maxima for the time of impact:

\frac{-k\,m\,v0+{k}^{2}\,p0-g\,{m}^{2}}{g\,k\,m}-\frac{m\,\mathrm{lambert\_w}\left(-\frac{\left(k\,m\,v0+g\,{m}^{2}\right) \,{e}^{\frac{-k\,m\,v0+{k}^{2}\,p0-g\,{m}^{2}}{g\,{m}^{2}}}}{g\,{m}^{2}}\right)}{k}

A lot of the terms duplicate, so in code it would not be as complex as it appears here. With it, you can find in one step (although an admitedly big step) the impact time. Oh provided you have an implementation of the Lambert W function.

With the values that created the graph at the start of this installment (p0=10; v0=20; g=-2; m=2; k=-0.5), our mega formula gives us the answer 14.91355897825015. Awesome.

One Response to “Ball with air friction (part 4)”

  1. Sébastien Says:

    Very nice serie of articles. Awesome. For the history, may it possible that you give some details about how ball is implemented in Kick Off 🙂

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: