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

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.

• Well $\theta(t) = \theta_0$ because the particle would just move radially outwards so $\theta$ should not change. Nov 13, 2018 at 15:18
• That's not true if the particle has a starting velocity $\textbf{v}_0$ with $\theta$ component Nov 13, 2018 at 15:20
• Oh! I thought you just left the particle there. Anyways I really don't think someone will answer this question on this site because it does not ask anything conceptual and therefore will not be useful for the community. Nov 13, 2018 at 15:27
• If you're doing a simulation, why are you analyticially solving it? Why not just use the equations of motion, and some sort of Euler method to do the time-steps? It probably will be less computationally intensive at the end of the day. Nov 13, 2018 at 16:17

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)$$.

• Thanks! I think that is indeed the best way. Given t I find $\eta$ and then find r! Nov 13, 2018 at 22:00

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.

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.

• Though you will need to specify radial and angular component of the velocity and the initial position. Nov 13, 2018 at 16:37