I want to simulate a simple orbit of a planet moving around a star which is fixed in position. I have formulated the ODEs for this problem using Lagrangian Mechanics and have found the equations of motion for a 2D orbit in polar coordinates.

I have found these to be:

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


$$r\ddot\theta+2\dot r\dot \theta = 0.$$

I would like some help formulating the problem in a way in which I can solve for $(r,\theta)$ using the RK4 method (I am aware this method is non-symplectic but I want to use this method first and will probably update to use velocity vertlet or something later on)

I have never used RK4 before to solve for two variables simultaneously and this is where my confusion in how to proceed arises.

  • $\begingroup$ I'm voting to close this question as off-topic because it's a request for codes rather than physics $\endgroup$ – Kyle Kanos Jul 28 at 18:40
  • $\begingroup$ @KyleKanos I disagree. The OP didn't ask for code. He only needed the concept how to transform second-order ODEs to a first-order ODE, as required for Runge-Kutta. $\endgroup$ – Thomas Fritsch Jul 28 at 18:49
  • $\begingroup$ Solving the two-body problem using cartesian coordinates in an inertial frame and without using relative coodinates is by far the most instructive way to do it. My students learn a lot of mechanics doing that way. Velocity Verlet algoritm is recommended, not oly because it is symplectic, which is important for a problem where angular momentum should be conserved, but also because it is very simple to implement and very efficient (one evaluation of force per step). $\endgroup$ – GiorgioP Jul 29 at 6:42
  • $\begingroup$ @ThomasFritsch then he should have opened any book on numerical methods where this surely is discussed. $\endgroup$ – Kyle Kanos Jul 29 at 11:26

To use numerical simulation program , you have to transformed the differential equations to first order differential equations (1)

$$\vec{\dot{y}}=\vec{f}(\vec{y},t)\tag 1$$


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

$$r\ddot{\theta}+2\,\dot{r}\dot{\theta}=0\tag 3$$

new variables :





equation (2):

$$\dot{y}_1-y_2\,y_3=-\frac{G M}{y_2^2}$$

equation (3):




$$\begin{bmatrix} \dot{y}_1 \\ \dot{y}_2\\ \dot{y}_3 \\ \dot{y}_4\\ \end{bmatrix}=\begin{bmatrix} -\frac{G M}{y_2^2}+y_2\,y_3 \\ y_1\\ -\frac{2\,y_1\,y_3}{y_2}\\ y_3 \end{bmatrix}$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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