Free fall

This is just elementary calculus, but it is harder than I thought at first glance. A friend told me that there is no formula expressing the altitude of a falling body if you is the inverse square law for the force, rather than making the approximation of a constant gravitational field. As far as I am able to tell, he is right. One can however compute some things.

We study the motion of a mass m, moving along an axis, whose position is denoted by x(t). At the time t=0, the mass is at rest in the position x(0)=x_0. There is a mass M, at the position x=0, that does not move and attracts m by a force f given by the inverse square law:
\displaystyle f=-\frac{GMm}{x^2},
G being the gravitational constant.

Let us write the conservation of the mechanical energy for the mass m: at any time t\ge 0,
\displaystyle \frac{1}{2}mx'(t)^2-\frac{GMm}{x(t)}=-\frac{GMm}{x_0}.
We consider a falling body, therefore x'(t)\le 0, and thus
This yields
and therefore
\displaystyle  -\frac{x'(t)x(t)}{\sqrt{x(t)(x_0-x(t))}}=k, (1)
\displaystyle k=\sqrt{\frac{2GM}{x_0}}
To find t as a function of x(t), we compute
\displaystyle \int\frac{x\,dx}{\sqrt{x(x_0-x)}}.
We have
\displaystyle \frac{x\,dx}{\sqrt{x(x_0-x)}}=\frac{-(-2x+x_0)}{2\sqrt{x(x_0-x)}}+\frac{x_0}{2\sqrt{x(x_0-x)}},
and therefore
\displaystyle \int\frac{x\,dx}{\sqrt{x(x_0-x)}}=-\sqrt{x(x_0-x)}+\frac{x_0}{2}\int \frac{dx}{\sqrt{x(x_0-x)}}.
We make the change of variable x=\frac{x_0}{1+u^2}, which gives
\displaystyle dx=\frac{2x_0u\,du}{(1+u^2)^2}
\displaystyle \sqrt{x(x_0-x)}-x_0\arctan\left(\sqrt{\frac{x_0}{x}-1}\right).
We obtain
\displaystyle \int\frac{x\,dx}{\sqrt{x(x_0-x)}}=-\sqrt{x(x_0-x)}-\frac{x_0}{2}\int_{u=\sqrt{\frac{x_{0}}{x}-1}} \frac{2\,du}{1+u^2}
and finally
\displaystyle \int\frac{x\,dx}{\sqrt{x(x_0-x)}}=-\sqrt{x(x_0-x)}-x_0\arctan\left(\sqrt{\frac{x_0}{x}-1}\right)+C.
We now integrate  (1) between 0 and t:
\displaystyle  k\,t=\sqrt{x(x_0-x)}+x_0\arctan\left(\sqrt{\frac{x_0}{x}-1}\right). (2)
This gives us the duration of the fall as a function of the position of the falling mass m.

We can find another formula for the antiderivative. Let us make the change of variable x=x_0\frac{1+\cos(u)}{2}=x_0\cos^2\left(\frac{u}{2}\right). Then
dx=-x_0\sin\left(\frac{u}{2}\right)\cos\left(\frac{u}{2}\right), and
\displaystyle \sqrt{x(x_0-x)}=x_0\sqrt{\cos^2\left(\frac{u}{2}\right)\left(1-\cos^2\left(\frac{u}{2}\right)\right)}=x_0\sin\left(\frac{u}{2}\right)\cos\left(\frac{u}{2}\right),
u being taken between 0 and \pi. We obtain
\displaystyle\int\frac{x\,dx}{\sqrt{x(x_0-x)}}=-\sqrt{x(x_0-x)}-\frac{x_0}{2}\int_{u=2\arccos\left(\sqrt{\frac{x}{x_0}}\right)} du,
and therefore
\displaystyle  k\,t=\sqrt{x(x_0-x)}+x_0\arccos\left(\sqrt{\frac{x}{x_0}}\right), (3)
which is the answer that can be found in Wikipedia. Taking x=0, we obtain the free-fall time t_f, the time it takes for the falling mass to reach the center of attraction:
\displaystyle t_f=\frac{\pi}{2}\sqrt{\frac{x_0^{3}}{2GM}}.

How does we find these changes of variable (this is the only tricky part in these calculations)? Well, if we set y=\sqrt{x(x_0-x)}, the points (x,y), with x in [0,x_0], are on the circle of radius x_0/2 centered at (x_0/2,0). Now if (x(u),y(u)) is a parametric representation of this circle, the change of variable x=x(u) allows us to write
\displaystyle \int \frac{dx}{\sqrt{x(x_0-x)}}=\int_{u=x^{-1}(x)}\frac{x'(u)\,du}{y(u)}.
In the first change of variable, we used a parametrization of the circle by rational functions:
\displaystyle\left\{\begin{array}{lcl}  x(u)&=&\frac{x_0}{1+u^2};\\  y(u)&=&\frac{x_0\,u}{1+u^2}.  \end{array}  \right.
In the second, we used a parametrization by trigonometric functions:
\displaystyle\left\{\begin{array}{lcl}  x(u)&=&x_0\frac{1+\cos(u)}{2};\\  y(u)&=&\frac{x_0}{2}\sin(u).  \end{array}  \right.
Both are natural choices (the second maybe more so than the first), and both produce a simple result.

If we set
\displaystyle F(\xi)=\xi\sqrt{1-\xi^2}+\arccos(\xi),
we obtain
\displaystyle k\,t=x_0\,F\left(\sqrt{\frac{x}{x_0}}\right)
and thus
\displaystyle  x(t)=x_0\,\left(F^{-1}\right)^2\left(\frac{k\,t}{x_0}\right). (4)
Indeed, F is continuous on [0,1], differentiable on [0,1[, and
\displaystyle F'(\xi)=\sqrt{1-\xi^2}-\frac{\xi^2}{\sqrt{1-\xi^2}}-\frac{1}{\sqrt{1-\xi^2}}=\frac{-2\,\xi^2}{\sqrt{1-\xi^2}}.
Therefore F is decreasing, and thus a bijection from [0,1] to [0,\pi/2]. Its inverse F^{-1} is continuous on [0,\pi/2] and differentiable on ]0,\pi/2[.

On the other hand, if we set
\displaystyle G(\xi)=\frac{\xi}{1+\xi^2}+\arctan(\xi),
we obtain
and thus
\displaystyle x(t)=\frac{x_0}{1+\left(G^{-1}\right)^{2}\left(\frac{k\,t}{x_0}\right)}.
The function G is continuous and differentiable on \mathbb{R} and
Its is therefore a diffeomorphism from \mathbb{R} to ]-\pi/2,\pi/2[. Let us denote by g the inverse G^{-1}. It satisfies the differential equation
\displaystyle g'(\zeta)=\frac{1}{2}+g(\zeta)^2+\frac{1}{2}g(\zeta)^4,
with g(0)=0. In particular, g'(0)=1/2. This, taken with Equation (4), implies that
\displaystyle x(t)=x_0\left(1-\frac{k^2t^2}{4x_0^2}+O\left(\frac{k^3t^3}{x_0^3}\right)\right)
for t close to 0. We have approximately
\displaystyle x(t)\simeq x_0-\frac{1}{2}\gamma_0t^2,
\displaystyle \gamma_0=\frac{GM}{x_0^2}.
The right-hand side is the result that we would obtain for a constant gravitational field \gamma_0.

We have
\displaystyle \lim_{\zeta\to \pi/2}g(\zeta)=+\infty.
\displaystyle G(\xi)=\frac{\xi}{1+\xi^2}+\arctan(\xi)=\frac{\pi}{2}-\frac{1/\xi}{1+(1/\xi)^2}-\arctan(1/\xi),
we have
\displaystyle \zeta=\frac{\pi}{2}-\frac{1/g(\zeta)}{1+(1/g(\zeta))^2}-\arctan(1/g(\zeta))\\    = \frac{\pi}{2}-\frac{1}{g(\zeta)}-\frac{1}{g(\zeta)^3}-\frac{1}{g(\zeta)}+\frac{1}{3\,g(\zeta)^3}+o\left(\frac{1}{g(\zeta)^3}\right)\sim -\frac{2}{3\,g(\zeta)^3},
and therefore
\displaystyle g(\zeta)\sim \left(\frac{3}{2}\left(\frac{\pi}{2}-\zeta\right)\right)^{-1/3}.
Therefore, when t tends to t_f,
\displaystyle x(t)\sim x_0\left(\frac{3}{2}\left(\frac{\pi}{2}-\frac{k\,t}{x_0}\right)\right)^{2/3}=\left(\frac{9GM}{2}\right)^{1/3}(t_f-t)^{2/3}.
Using the differential equation, we also find
\displaystyle  x'(t)\sim-\sqrt{2GM}\left(\frac{9GM}{2}\right)^{-1/6}(t_f-t)^{-1/3}\\  =-\frac{2}{3}\left(\frac{9}{2}\right)^{1/3}(GM)^{1/3}(t_f-t)^{-1/3},
which is what we would obtain by formally differentiating the equivalent.

The free-fall time can be rewritten
\displaystyle t_f=\frac{\pi}{2}\sqrt{\frac{x^{2}}{2GM}}\,\frac{x_0}{x}\sqrt{x_0}=\frac{\pi x_0}{2x}\sqrt{\frac{x_0}{2\gamma}}
with \gamma the accelaration of gravity at the altitude x. How many time does it take to fall from the altitude of the Moon? We take x to be  the radius of the Earth, and x_0 is the distance form the Earth to the Moon. Wikipedia gives
x_0=384399 kilometers, x=6378 kilometers and \gamma=9.780327 meters per second squared. We find roughly 116 hours, which is to say four days and twenty hours.

The Wikipedia article I linked above gives two references for the free fall:

1. From Moon-fall to motions under inverse square laws, by S. K. Foong, European Journal of Physics, 2008 29: 987–1003;
2. Radial motion of Two mutually attracting particles, by Carl E. Mungan, The Physics Teacher, 2009, 47: 502-507.

I skimmed the first article (unfortunately it is behind a paywall, but I could get access with my university subscription). The formula for x \ll x_0 and the numerical value for the free fall time from the orbit of the moon agree with what I found. There is also a series expansion, and probably other stuff. I will read in detail later.


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