Solving the Free Fall Differential Equation

[Meithan]

I've been trying to derive the full solution to the problem of free fall when the distances considered are large (and hence, the acceleration is no constant, but a function of height). This isn't homework, I was just curious.

The problem is innocently simple, but yields a differential equation which I cannot solve (and it's in principle solvable).

Free fall in a gravitational field

Consider a large mass M, which we place at the origin of a 1D coordinate system. We then drop a test mass m from a certain height r=h_0[/itex], from rest. The problem is to calculate its position as a function of time, r(t)[/itex]. We assume that M \gg m so the larger mass remains fixed at the origin and does not move.<br /> <br /> An example of this would be a high-orbit satellite free-falling onto the Earth in a radial trajectory.<br /> <br /> <b>My attempt at a solution</b><br /> <br /> The physics of the problem is very straightforward. At a distance r from the origin (where the larger mass is located), the local acceleration due to the gravitational force is given by<br /> <br /> g(r) = -\frac{G M}{r^2}<br /> <br /> And so the (differential) equation of motion of the test mass is<br /> <br /> \frac{d^2r}{dt^2} = -\frac{G M}{r^2}<br /> <br /> or, put simply,<br /> <br /> r^2 \frac{d^2r}{dt^2} = -\mu<br /> <br /> where \mu \equiv GM is the standard gravitational parameter of the system. This second-order ODE is not linear, so its solution might be complicated, but it's not impossible.<br /> <br /> We define v \equiv dr/dt so that d^2r/dt^2 = dv/dt, and using the chain rule we write the ODE as:<br /> <br /> r^2 \frac{dv}{dr} \frac{dr}{dt} = -\mu<br /> <br /> or, since dr/dt = v,<br /> <br /> v r^2 \frac{dv}{dr} = -\mu<br /> <br /> In this form, the ODE is separable:<br /> <br /> v dv = -\frac{\mu}{r^2} dr<br /> <br /> Integrating, we get<br /> <br /> \frac{v^2}{2} = \frac{\mu}{r} + C<br /> <br /> where C is an integration constant. Since we assume the object is dropped from rest from a height h_0, we can set v(r=h_0)=0 to obtain C[/itex]:<br /> <br /> 0 = \frac{\mu}{h_0} + C<br /> <br /> and hence<br /> <br /> C = -\frac{\mu}{h_0}<br /> <br /> The solution of the ODE for the velocity is therefore (<b>edit</b>: I just found out that this can be obtained very easily by an energy consideration):<br /> <br /> \frac{v^2}{2} = 2\mu \left( \frac{1}{r} - \frac{1}{h_0} \right)<br /> <br /> which gives the velocity as a (ugly) function of height:<br /> <br /> v(r) = -\sqrt{2\mu \left( \frac{1}{r} - \frac{1}{h_0} \right)}<br /> <br /> (we take the negative root because we expect the velocity to be negative as the particle falls). The problem I encounter is when I try to integrate this equation:<br /> <br /> -\int \frac{dr}{\sqrt{2\mu \left( \frac{1}{r} - \frac{1}{h_0} \right)}} = \int dt<br /> <br /> The right-hand side is trivial, and we can wright the left-hand side as follows:<br /> <br /> -\frac{1}{\sqrt{2\mu}} \int \frac{dr}{\sqrt{\frac{\displaystyle 1}{r} - A}} = t + C<br /> <br /> where A \equiv 1/h_0 is a positive constant. <br /> <br /> Here's where I'm stuck. I've no clue how to solve that integral. I gave it to Maxima, and I get a pretty ugly solution. Can anyone help? is there another way of solving the original differential equation?

 

[zeebek]

expression number 2, where you write the Newtons second law in shperical coordinates is wrong

<br /> \frac{d^2r}{dt^2} = -\frac{G M}{r^2}<br />

acceleration is not just

<br /> \frac{d^2r}{dt^2}<br />

it should include angular derivative term.

<br /> \frac{d^2r}{dt^2} - r (\frac{d \theta}{dt})^2<br />

So this is not "very straightforward"...

 

[Meithan]

I was writing it in cartesian coordinates! I just chose to call the linear distance, measured from the center of mass of the Earth, r. I could've used y or any other letter, it's entirely arbitrary. Regardless, if you choose to write it in spherical coordinates, it would still be correct, because then you'd be describing a purely radial trajectory, for which d\theta / dt = 0 (the angular position doesn't change with time).

 

[Meithan]

I've done some advances towards solving this problem (though not the kind I wanted).

1) I used Mathematica to calculate the integral that I can't solve. The solution is ugly, but it's a solution. However, this gives time as a function of position, but we want the other way around. Since the solution t(x) is so complicated, I don't think there's a way to invert it to get x(t). When I get some more time, I'll post Mathematica's solution. I'd still like to know exactly how one calculates the integral.

2) There's an http://en.wikipedia.org/wiki/Free_fall#Inverse-square_law_gravitational_field" explaining this as the special case of pure radial motion in a central, inverse-square field. If the particle doesn't exceed the escape energy, its motion is described by an "elliptical radial orbit" -a degenerated ellipse with a semi-minor axis of length zero-, a special case of the general solution (which, you might recall, is a conic section). The full solution for t(x) is in that Wikipedia article, and it's ugly as hell. Likewise, another article (which I can't recall) states that this function is not algebraically invertible: its inverse can only be written as a (very ugly) power series.

 

[D H]

<br /> -\int \frac{dr}{\sqrt{2\mu \left( \frac{1}{r} - \frac{1}{h_0} \right)}} = \int dt<br />

Here's where I'm stuck. I've no clue how to solve that integral. I gave it to Maxima, and I get a pretty ugly solution. Can anyone help? is there another way of solving the original differential equation?

Maxima, Mathematica, and Maple can be pretty stupid at times.

First a change of variables to normalize things. Define s=r/h_0. Then

-\int \frac{dr}{\sqrt{2\mu \left( \frac{1}{r} - \frac{1}{h_0} \right)}}<br /> = \sqrt{\frac{h_0^3}{2\mu}}<br /> \int \frac{-ds}{\sqrt{\frac 1 s - 1}}

Mathematica does indeed come up with something ugly for that integral.

That radical suggests a trigonometric substitution. Being of the form (something)-1, making that (something) be \sec^2\theta will get rid of the radical. So, define \theta such that \sec^2\theta = 1/s, or

\aligned<br /> s &amp;= \cos^2\theta \\<br /> ds &amp;= -2\cos\theta\sin\theta\,d\theta<br /> \endaligned

The radical becomes

\sqrt{\frac 1 s - 1} = \sqrt{\sec^2\theta - 1} = \sqrt{\tan^2 \theta} = \tan\theta

That final step should be \pm \tan\theta in general. That nasty ± becomes + upon restricting theta to 0 to pi/2. That corresponds to r falling from the initial height to zero. With this restriction, the integral becomes

\aligned<br /> \int \frac{-ds}{\sqrt{\frac 1 s - 1}} &amp;=<br /> \int \frac{2 \sin\theta \cos\theta \, d\theta}{\tan\theta} \\<br /> &amp;= \int 2\cos^2\theta\,d\theta \\<br /> &amp;= \sin\theta\cos\theta + \theta \\<br /> &amp;= \sqrt{s(1-s)}+ \arccos\sqrt s<br /> \endaligned

I left the constant of integration out of the above; the constant of integration is zero here. Backsubstituting the original variable r yields

t(r) = \sqrt{\frac{h_0^3}{2\mu}}<br /> \Biggl(\sqrt{\frac r{h_0}\left(1-\frac r{h_0}\right)}+<br /> \arccos\left(\sqrt{\frac r{h_0}}\,\right)\Biggr)

 

[Meithan]

Fantastic! I'm ashamed I did not realize a trig substitution would handle that integral! Thanks for the enlightenment. I suppose, as the Wikipedia suggests, that this result is not invertible in terms of elementary functions. So that would be as far as one can go in this problem.

Now I'll try to solve it for the more general case, where one is given arbitrary initial conditions (position and velocity). It should be similar.