Consider two-body central force problem in polar co-ordinates $r, \theta$. Corresponding 2nd order differential equation is obtained by using conservation of angular momentum. This equation is :

$$ \frac{d^2 r}{dt^2} = \frac{l^2}{m^2r^3} - \frac{GM}{r^2}$$

where $r(t)$ is radial position of particle (of mass $m$) as a function of time $t$, $l$ is angular momentum which is constant, $G$ is gravitational constant & $M$ is mass of the heavier body, assumed to be at rest at the origin of co-ordinate system i.e. at $(r,\theta)=(0,0)$.

I want to solve above non-linear differential equation; it is non-linear since dependent variable $r$ has powers -3 and -2 on RHS.

Can I use 4-order Runge-Kutta method to solve this equation ?

  • 2
    $\begingroup$ I'm voting to close this question as off-topic because it's about implementation details of a computational problem. $\endgroup$
    – David Z
    Mar 7, 2017 at 4:36
  • 1
    $\begingroup$ @OfekGillon: $\theta$ is not ill-defined. Actually we have two different 2-order differential equations (coupled) : one for $r$ and another for $\theta$. Conservation of angular momentum de-couples them and reduces to one equation given above. ...Also if we try to solve above 1-Dim equation analytically, we end up in the solution of the form $t(r)$ i.e. time is function of $r$. So we have to invert that into $r(t)$. And this inversion process can be extremely difficult in practice. Please see standard textbook on classical mechanics e.g. by Goldstein (Chapter 3). $\endgroup$
    – atom
    Mar 7, 2017 at 4:48
  • 3
    $\begingroup$ @DavidZ: Seems to me this falls under point #1: We can advise on numerical methods when tied to a physical problem and not a generic P/ODE. $\endgroup$
    – Kyle Kanos
    Mar 7, 2017 at 11:06
  • 2
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Mar 7, 2017 at 12:44
  • 2
    $\begingroup$ What makes you think it couldn't be used? Have you looked at what other people simulating bodies in orbit use? $\endgroup$
    – tpg2114
    Mar 7, 2017 at 17:06

2 Answers 2


For that kind of problem a symplectic integrator will give more physically accurate results, since solutions will preserve the energy and angular momenta. You will get a better long term predictions.




Yes you can. The only caveat with numerically solving non-linear ODE's exists when the derivative function is discontinuous in some way. For example a sliding block with dry friction. The direction of friction can switch at an instant and care must be taken not to step over such event. Similarly with collisions.

In your case, I think the derivative function is smooth and continuous if the radial distances $r$ are large enough. To get good results I would focus a lot on time step management. Remember the RK4 scheme is equivalent to fitting cubic splines into the motion and so care must be taken to set the time steps so the results have sufficient geometrical accuracy.


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.