In 1D the force of gravity does not diminish with distance however here I use 1D to mean the point mass only moves in a line in 3D space, with the other mass stationary at the origin.

I started by putting the equation $y''(t)=-y/abs(y)^3$ into a numerical ODE solver with the starting velocity set to 0 and at first the trajectory falls as expected but when the point mass hits the origin it flies off into the sky... Presumably because of the singularity at $y=0$. So I wrote a Julia script removing the singularity:

for t in 1:100
  if abs(r)>1
    vd=a(r) #acceleration is treated as constant wrt time in a single second and then the velocity difference should be this

Now the particle passes through the origin with an escape velocity and charges onward on the other side indefinitely, albeit with a decreasing velocity which should be enough to turn it around yet doesn't.

If I set $abs(r)>0.1$ the particle dips below zero and then shoots up into the sky like with the ODE solver...

A particle that starts at rest can't gain an escape velocity from the gravitational acceleration, so what is wrong with the simulation?

  • $\begingroup$ You need to use a symplectic integrator en.wikipedia.org/wiki/Symplectic_integrator and treat the step through the center with an exact solution (use the symmetry of the problem). OTOH, this problem has an exact solution, so there is no need for numerics to begin with. $\endgroup$
    – FlatterMann
    Apr 24, 2023 at 8:12
  • $\begingroup$ @FlatterMann The wikipedia entry is complicated, any concrete code example? $\endgroup$
    – Qni
    Apr 24, 2023 at 8:31
  • $\begingroup$ @FlatterMann Before turning to numerical methods I tried inputting the equation into WolframAlpha which failed to give a solution. What is the exact solution? $\endgroup$
    – Qni
    Apr 24, 2023 at 8:33
  • $\begingroup$ See this related question on radial trajectories: physics.stackexchange.com/q/19388/123208 and links therein. $\endgroup$
    – PM 2Ring
    Apr 24, 2023 at 9:12
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Apr 24, 2023 at 12:01

1 Answer 1


what is wrong with the simulation ?

Your problem is that you are using the Euler method to derive your numerical solution. Although this method is simple to implement, it is inaccurate and unstable. In your example, as $|r|$ decreases, the inaccuracy created when the method assumes that $a$ is constant between one time step and the next becomes increasingly large. You can reduce these inaccuracies by using a time step that is smaller than one second, but then you need more steps to simulate a given time period. Alternatively you could use a more complex but more accurate numerical integration method, such as shortening the length of your time step as $|r|$ gets smaller, or using one of the Runge-Kutta methods.

There is also a theoretical problem, which is that the behaviour of your point mass at $r=0$ is undefined, since its acceleration becomes "infinite" at that point - although your code addresses this in a rather ad-hoc fashion by forcing the acceleration to be zero if $|r| <= 1$.

  • 1
    $\begingroup$ Thanks! When I decreased to milliseconds it produced a back and forth orbit $\endgroup$
    – Qni
    Apr 24, 2023 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.