Is there a closed-form solution r(t) regarding universal law of gravitation?

[JDoolin]

Newton's universal law of gravitation:

$$F=-G \frac{m_1 m_2}{r^2}$$

I'd like to set up the problem so the particle begins at t=0 at radius r=r₀ and radial velocity v_r=v₀. And there is only a component of velocity, in the radial direction. (The particle is going straight toward the gravitating body)

The above equation becomes

$$F=m a_r= m_2 \frac{\mathrm{d}^2 r}{\mathrm{d} t^2}= G \frac{m_1 m_2}{r^2}$$

And so the most simply I can give the differential equation is:

$$\left (\frac{\mathrm{d}^2 }{\mathrm{d} t^2} \right )r = -G \frac{ M}{r^2}$$

Does this equation have a closed-form solution?

If not, is there a well-known technique of solving it numerically? Here is what I did numerically.

(1) determine r0 and v0.
(2) determine the acceleration at your current radius a=-GM/r^2
(3) choose a small Δt during which the acceleration is approximately constant.
(4) Find a final velocity by multiplying the acceleration times the Δt.
(5) Find the average velocity from (initial velocity + final velocity)/2
(6) Multiply average velocity by Δt to find Δr.
(7) Add Δt and Δr and Δv to your previous values of t, r, and v, and start over at step (2)

One particular set-up is to determine the initial velocity at r₀ to be the negative of the escape velocity at that same r₀. It stands to reason that if a particle is traveling at a speed equal and opposite to its escape velocity as it falls in, that conservation of energy will maintain this, so that the particle will fall in at its escape velocity the entire way.

I tried out this numerical method and the velocities calculated numerically did match the escape velocities pretty well. But this gives us a function v(r) when I am wanting a function r(t), and it also is a solution that requires a specific pre-defined initial radial velocity, where I would like a function where the initial velocity can be decided arbitrarily.

I appreciate any help anyone might think of. Thanks.

 

[Mute]

There's not really a nice closed form for r(t). You can get an expression for t(r), though, by multiplying

$$\ddot{r} = -GM\frac{1}{r^2}$$

by $\dot{r}$ to get

$$\ddot{r}\dot{r} = -GM\frac{\dot{r}}{r^2}$$

Then notice that $\ddot{r}\dot{r} = (1/2) d(\dot{r}^2)/dt$ and $\dot{r}/r^2 = -d(1/r)/dt$.

You can then integrate the equation once, and then you will be able to separate the resulting equation and integrate it, but the best you will be able to get is $f(r(t)) = t$ for some function f.

 

[JDoolin]

There's not really a nice closed form for r(t). You can get an expression for t(r), though, by multiplying

$$\ddot{r} = -GM\frac{1}{r^2}$$

by $\dot{r}$ to get

$$\ddot{r}\dot{r} = -GM\frac{\dot{r}}{r^2}$$

Then notice that $\ddot{r}\dot{r} = (1/2) d(\dot{r}^2)/dt$ and $\dot{r}/r^2 = -d(1/r)/dt$.

You can then integrate the equation once, and then you will be able to separate the resulting equation and integrate it, but the best you will be able to get is $f(r(t)) = t$ for some function f.

Thank you.

That is helpful. From your replacements I derived:

$$\frac{1}{2} \left (v_f^2 - v_i^2 \right ) = G M\left ( \frac{1}{r_f}-\frac{1}{r_i} \right )$$

This satisfies one of my criteria, because it is a solution that let's you set your initial location and velocity arbitrarily, and determines final velocity as a function of final location.

But was there something more you had in mind? Because I didn't come up with t as a function of r or t=f(r(t)). When I made the replacements you suggested, the dt's canceled out.

 

[Mute]

Thank you.

That is helpful. From your replacements I derived:

$$\frac{1}{2} \left (v_f^2 - v_i^2 \right ) = G M\left ( \frac{1}{r_f}-\frac{1}{r_i} \right )$$

You could also shortcut your way to this by using conservation of energy. Note that you have a minus sign error on the right hand side. It should be 1/r_0 - 1/r_f.

This satisfies one of my criteria, because it is a solution that let's you set your initial location and velocity arbitrarily, and determines final velocity as a function of final location.

But was there something more you had in mind? Because I didn't come up with t as a function of r or t=f(r(t)). When I made the replacements you suggested, the dt's canceled out.

The equation you derived, for arbitrary final time, is also

$$\frac{1}{2}(\dot{r}^2(t) - \dot{r}(0)) = GM\left(\frac{1}{r(0)} - \frac{1}{r(t)}\right)$$

If you solve for $\dot{r}(t)$, you have another equation to solve that in principle get you r(t), but it looks like $dr/dt = g(r)$, for some function g(r), so when you solve it you get f(r(t)) = t, where f(r) is the integral of 1/g(r), and cannot really be inverted explicitly (in terms of elementary functions) to get $r(t) = f^{-1}(t)$.

 

[JDoolin]

You could also shortcut your way to this by using conservation of energy. Note that you have a minus sign error on the right hand side. It should be 1/r_0 - 1/r_f.

I'm fairly certain the original sign is right. You want it so if the initial radius is larger than the smaller radius, you come out with a positive value for the change in velocity.
Also, the conservation of energy law is $$\Delta KE + \Delta PE = \Delta E_{non-cons}$$ where $$\begin{align*} PE&=-G\frac{m_1 m_2}{r} \\ KE &= \frac{1}{2}m v^2 \end{align*}$$
( Potential Energy's zero point can be set wherever you want, but by common convention, it is set so PE->0 as r-> infinity)
(and of course the non-conserved energy in this case is zero)

In any case, the overall process seemed to be an elegant way to derive the conservation of energy law.

The equation you derived, for arbitrary final time, is also

$$\frac{1}{2}(\dot{r}^2(t) - \dot{r}(0)) = GM\left(\frac{1}{r(0)} - \frac{1}{r(t)}\right)$$

If you solve for $\dot{r}(t)$, you have another equation to solve that in principle get you r(t), but it looks like $dr/dt = g(r)$, for some function g(r), so when you solve it you get f(r(t)) = t, where f(r) is the integral of 1/g(r), and cannot really be inverted explicitly (in terms of elementary functions) to get $r(t) = f^{-1}(t)$.

That's good. Once this is done, (and I decided to let the initial velocity be zero for simplicity) I have a differential equation of the form

$$\left ( \frac{dr}{dt} \right )^2-\frac{2 G M}{r}+\frac{G M}{r_0} =0$$


$$\left ( \frac{dr}{dt} \right )^2-\frac{c_1}{r}+c_0=0$$

I'm still not sure how to solve this thing, but that's one step closer; at least now it's just a first derivative.

 

[jackmell]

You take the square root and separate variables and obtain two solutions

 

[Ben Niehoff]

The central motion problem in a 1/r potential is solved in graduate texts on classical mechanics such as Goldstein. You can get a nice solution t(r) for t as a function of r.

 

[eztum]

Who is more deeply interested in the problem should study the exact solution based on Kepler's famous transcendental equation. Compared with this, the present approach is a 'toy solution'. Comparing the poor physics textbook method with Kepler's method allows one to appreciate what genius in science means.

 

[dextercioby]

Who is more deeply interested in the problem should study the exact solution based on Kepler's famous transcendental equation. [...]

Do you have a reference for this ? Book or freely available article.

 

[eztum]

Do you have a reference for this ? Book or freely available article.

The Wikipedia article 'Kepler orbit' has this and much more. Notice that solving Kepler's equation solves the problem where on his orbit the planet is at a given instant of time. Solving this equation by some iteration scheme is easy and fast and the computational burden is independent of how many centuries the this point in time is in the future or in the past. So there is no fundamental difference to the situation that you express the planetary position by an explicit 'closed formula'. Solving it by a series expansion is the historical origin of the Bessel functions (Bessel was an astronomer). I just read Kepler's original text on the matter in a facsimile edition of selected works of Kepler. Really not easy to read (although his language, German, is mine).

Physics textbooks normally are content with deducing the geometry of the orbits (conic sections) from Newton's law and neglect the question how to find the position of the body on its orbit efficiently. Texts on Celestial Mechanics and Astronomy are better sources for this problem.

 

[eztum]

If not, is there a well-known technique of solving it numerically? Here is what I did numerically.

(1) determine r0 and v0.
(2) determine the acceleration at your current radius a=-GM/r^2
(3) choose a small Δt during which the acceleration is approximately constant.
(4) Find a final velocity by multiplying the acceleration times the Δt.
(5) Find the average velocity from (initial velocity + final velocity)/2
(6) Multiply average velocity by Δt to find Δr.
(7) Add Δt and Δr and Δv to your previous values of t, r, and v, and start over at step (2)

You can considerably increase the power of such a numerical method by
following Verlet (actually the method is much older and is often associated with
Störmer, sometimes it is called explicit midpoint method)

Treat r (position) and v (velocity) as vectors (in 3-space, or 2-space).
(1) determine r0 and v0.
(2) determine the acceleration am at position rm = r0 + (Δt/2)v0
(3) find the final velocity vf = v0 + Δt am
(4) find the final position rf = rm + (Δt/2)vf
(5) restart the loop with r0 = rf and v0 = vf

 

[JDoolin]

You take the square root and separate variables and obtain two solutions

The central motion problem in a 1/r potential is solved in graduate texts on classical mechanics such as Goldstein. You can get a nice solution t(r) for t as a function of r.

For completion:


$$\begin{align*} \left ( \frac{dr}{dt} \right )^2&=\frac{c_1-c_0 r}{r} \\ \left ( \frac{dr}{dt} \right )&=\sqrt{\frac{c_1-c_0 r}{r}}\\ \frac{dt}{dr}&=\sqrt{\frac{r}{c_1-c_0 r}}\\ \int dt&=\int_{r_0}^{r_f}\sqrt{\frac{r}{c_1-c_0 r}}dr \end{align*}$$

where c1=2 G M, c0 = G M /r0

The last line, evaluated via Mathematica was:


$$\frac{\sqrt{\frac{r}{\text{c1}-\text{c0} r}} \left(\sqrt{\text{c0}} \sqrt{r} (\text{c0} r-\text{c1})+\text{c1} \sqrt{\text{c1}-\text{c0} r} \tan ^{-1}\left(\frac{\sqrt{\text{c0}} \sqrt{r}}{\sqrt{\text{c1}-\text{c0} r}}\right)\right)}{\text{c0}^{3/2} \sqrt{r}}$$

A closed form, but not a pretty one.

Edit: Just a note about the domain of the arctan function.

$$\sqrt{\frac{c_0 r}{c_1-c_0 r}} = \sqrt{\frac{\frac{G M}{r_o}}{2 G M -\frac{G M}{r_o}r}}=\sqrt{\frac{r}{2 r_0 - r}}$$

This will be real-valued so long as $0<r<r_0$. That makes sense, because based on the assumption of v=0 when r=$r_0$, there would be no time such that $r>r_0$.

Also, you would get two times for each radius, a negative time on the way up, and a positive time on the way down, where t=0 would represent the time where the object reaches the apex.

Also, the form as I've given it isn't quite general, because by setting the velocity to zero at ANY given radius, I've left out all the particles that are traveling at speeds greater than the escape velocity.