Equation of motion for Central Force problem with repulsive inverse square law

Question

This seems to be a very simple question, but it has been bugging me for days. I have found many material related to this, including the great Goldstein book, but none get quite were I want. I'm starting to think it's impossible.

I have a classic one body problem on a central force field, repulsive, $1/r^2$ (let's say coulomb). I know the motion is in the plane, and I have $r0(r, theta)$ and $v0(r, theta)$. So I can get $E$ and $l$ as Goldstein and many others suggest and I get to the integral he presents for $$t = \int_{r_0}^r \frac{dr}{\sqrt{\frac{2E}{m} - \frac{2k}{mr} - \frac{e^2}{m^2r^2}}}.$$ Ok, I solved it with Sage (no problem for me), got:

$$ t(r) = \frac{\sqrt{2} \sqrt{E m} k \log\left(4 \, E m r - 2 \, k m + 2 \, \sqrt{2} \sqrt{2 \, E m r^{2} - 2 \, k m r - e^{2}} \sqrt{E m}\right)}{4 \, E^{2}} - \frac{\sqrt{2} \sqrt{E m} k \log\left(4 \, E m r_{0} - 2 \, k m + 2 \, \sqrt{2} \sqrt{2 \, E m r_{0}^{2} - 2 \, k m r_{0} - e^{2}} \sqrt{E m}\right)}{4 \, E^{2}} + \frac{\sqrt{2 \, E m r^{2} - 2 \, k m r - e^{2}}}{2 \, E} - \frac{\sqrt{2 \, E m r_{0}^{2} - 2 \, k m r_{0} - e^{2}}}{2 \, E} $$

Ughh. That looks bad. But ok, no problem, I'm going to compute it with a program. But I need to finally to invert this monster to get $r(t)$, to finally get $\theta(t)$. Again, I just want the equations of motion (not the orbit equation!). When did that became this hard? I really thought it should be easier, it's a simple problem (apparently). I tried to use solve from Sage to work it out, but it simply doesn't work. It seems Sage's solve doesn't solve when he can't find a way. And I don't blame him.

My question is, how to find $r(t)$ and $\theta(t)$ for a single particle on a static central force $F(\textbf{r}) = + \frac{k}{r^2} \hat{\textbf{r}}$ (repulsive), given initial conditions $\textbf{r}_0$, $\textbf{v}_0$, $E$, $l$ in polar coordinates $\textbf{r} = (r, \theta)$ in the plane defined by $\textbf{r}_0$ and $\textbf{L}$?

Background: I'm writing a simulation for this scenario. I'm doing it numerically, I iterate over very small $dt$'s and assume in then the acceleration is constant, and that works OK to showing the simulation with a single body as a video to the user. But my end goal is to simulate thousands of independent particles, and I'm interested in some characteristics of the final outcome. So had I had a formula for $r(t)$, I could skip straight to larger $t$'s! Also, I'd like to render tons of particles on the screen, making a quicker video, because I would not need to update countless times each render passage. As of for now, I need to set $dt$ very small to get accurate results (it changes drastically on $dt$ chosen), and that gives hundreds of update ticks before each render, dropping the frame rate considerably for 1000+ particles. A single formula would be a bliss.

Answer 1

There isn't a simple formula for $r(t)$ but there is a very simple parameterization $(r(\eta),t(\eta))$. This is the repulsive version of the Kepler equation that Ben Crowell referred to. There is also a simpler formula for $t(r)$ than you got. Here's how you derive these.

The integral for $t(r)$ is

$$t=\int_{r_0}^r\frac{dr}{\sqrt{\frac{2E}{m}-\frac{2k}{mr}-\frac{L^2}{m^2 r^2}}}$$

where $E$ is the conserved energy (positive), $L$ is the conserved angular momentum (either positive or negative), and $k$ is the constant for the potential, $V(r)=k/r$ (which is positive for a repulsive inverse-square force).

We can rewrite this as

$$Ct=\int_{r_0}^r\frac{r\,dr}{\sqrt{r^2-2Ar-B^2}}$$

where $A=k/2E$, $B=L/\sqrt{2mE}$, and $C=\sqrt{2E/m}$; $A$, $B^2$, and $C$ are all positive constants.

As you've found, a brute-force integration produces a mess, so the trick is to change variables to make the integral look nicer.

Let $$r=A+\sqrt{A^2+B^2}\cosh\eta.$$ Then the integral becomes

$$Ct=\int_{\eta_0}^\eta(A+\sqrt{A^2+B^2}\cosh\eta)\,d\eta$$

so

$$Ct=A\eta+\sqrt{A^2+B^2}\sinh\eta$$

Thus $r(\eta)$ and $t(\eta)$ are simple formulas giving a nice parameterization of $r(t)$ in terms of a parameter $\eta$.

You can invert $r(\eta)$ to get

$$\eta=\cosh^{-1}{\frac{r-A}{\sqrt{A^2+B^2}}}$$

and then substitute this into $t(\eta)$ to get

$$t(r)=\frac{1}{C}\left(A\cosh^{-1}{\frac{r-A}{\sqrt{A^2+B^2}}}+\sqrt{A^2+B^2}\sinh\cosh^{-1}{\frac{r-A}{\sqrt{A^2+B^2}}}\right)$$

I think your best approach will be to use the parametric formulas. For a given $t$, you can numerically solve $t(\eta)$ to get $\eta$ to the precision you need, and then substitute it into $r(\eta)$.

Answer 2

The attractive case is hundreds of years old and is called the Kepler problem. I believe the repulsive case works out similarly, maybe requiring tweaks here and there to account for the different sign. It's a problem that has been of interest in understanding cases like scattering of alpha particles by nuclei, which is how the nucleus was originally discovered. Getting the position as a function of time requires the solution of a transcendental equation called Kepler's equation, which can be efficiently solved by iteration.

It may also be quite accurate and efficient to simply solve the equations of motion using Runge-Kutta. Using the Euler method as you describe is likely to give a poor trade-off of precision versus efficiency. There are lots of free implementation of Runge-Kutta out there -- it's not a good idea to roll your own unless there is some reason why you can't use one of those.

Answer 3

Numerical analysis is better suited for this problem if all you want is to form a simulation. Euler's method or runge kutta method can be used to solve the equations of motion numerically and thus find a solution for r and $\theta$ as functions of time.

$\ddot r=k/r^2$

$\ddot\theta=-2\dot r\dot\theta/r$

Euler's method says

$ \dot r(t+\delta t) =\dot r(t) +h*\ddot r $ $r(t+\delta t) =r(t) + h*\dot r(t+\delta t) $

Where h is a small number say

h=$\delta t$=0.001

You can use this algorithm and a similar algorithm for $\theta (t) $ in a loop to find the trajectory of motion numerically.