I'm supposed to integrate differential equations for $r$ and $\theta$ in order to simulate orbital motion. The differential equation I used for $r$ is second order and $\theta$ is first order. The end goal is to graph $r(t)$, $\theta(t)$, which I'll probably do with matplotlib in python. I was wondering if anyone with prior knowledge of orbits and Runge-Kutta could let me know where I went wrong with my code? It compiles without error but isn't returning values that I expected. This is my code with parameters describing Earth's orbit:

#include <stdio.h>
#include <math.h> 
#define N 3         /* number of first order equations */
#define dist 0.00001        /* stepsize in t*/
#define MAX 1       /* max for t */
#define MAXTHETA 2*M_PI

FILE *output;           /* internal filename */

void runge4(double x, double y[], double step);
double f(double x, double y[], int i);

int main()
double t, y[N];
int j;

void runge4(double x, double y[], double step); /* Runge-Kutta function */

double f(double x, double y[], int i);      /* function for derivatives */

output=fopen("kepler.rtf", "w");                   /* external filename */

y[0]= 1.02343552273569;                                       /* initial position */
y[1]= 6.18068911;                                       /* initial velocity */
y[2] = 0;                                               /* initial angle */

fprintf(output, "0\t%f\t%f\n", y[0], y[2]); /* prints time and position */

for (j=1; j*dist<=MAX ;j++)         /* time loop */
    runge4(t, y, dist);

    fprintf(output, "%f\t%f\t%f\n", t, y[0], y[2]); /* prints time, position, theta */


void runge4(double x, double y[], double step)
double h=step/2.0,          /* the midpoint */
    t1[N], t2[N], t3[N],        /* temporary storage arrays */
    k1[N], k2[N], k3[N],k4[N];  /* for Runge-Kutta */
int i;

for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
for (i=0;i<N;i++) t3[i]=y[i]+    (k3[i]=step*f(x+h, t2, i));
for (i=0;i<N;i++) k4[i]=        step*f(x+step, t3, i);

for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;

double f(double x, double y[], int i)
    double k = 39.4784176044, m = 0.000003003, h = 0.00000300299, l = 0.0000188781;
    if (i==0) return(y[1]);         /* derivative of first equation (dr/dt)*/
    if (i==1) return((-k/(m*y[0]*y[0])) + (l*l)/(m*y[0]*y[0]*y[0])); /* derivative of second equation (d^2r/dt^2) */
    else return(h/(y[0]*y[0])); /* derivative of theta (d(theta)/dt) */

The second order differential equation I used for the position, $r$, is:

$$\ddot{r} = \frac{-k}{mr^2} + \frac{l^2}{mr^3}$$

where the parameter $l$ is the angular momentum and $k$ is the gravitational constant multiplied my the mass of the object that the planet is orbiting, in this case the sun. The first order differential equation I used for $\theta$ is:

$$\dot{\theta} = \frac{h}{r^2}$$

$$h = l\left(\frac{m_1 m_2}{m_1 + m_2}\right)^{-1}$$

where $m_1$ and $m_2$ are the masses of the two bodies in the system, and the mass term in the equation for $h$ is the reduced mass. All of the units of parameters are in AU, years, and solar masses. If anybody could point me in the right direction that would be awesome.

  • 5
    $\begingroup$ I cannot comment on your code, but want to point out that the choice of time integration scheme for such problems is very important. You need an RK scheme that conserves the energy, otherwise you will get very wrong solutions. See for example en.wikipedia.org/wiki/Symplectic_integrator $\endgroup$
    – cfdlab
    Jul 19, 2016 at 14:05
  • 2
    $\begingroup$ What's the issue? You said you didn't get the numbers you were expecting. What happens if you use a lower Dt? Isn't blowing up? This needs more information. $\endgroup$ Jul 19, 2016 at 17:34
  • $\begingroup$ The issue is that the position keeps on decreasing into negative numbers, and theta stays 0. I just can't see what's giving these incorrect values. $\endgroup$
    – mathrox
    Jul 19, 2016 at 17:39
  • $\begingroup$ The second order equation for $r$ is decoupled from the equation for $\theta$. Why don't you multiply the equation for $r$ by $\dot{r}$ on both sides, integrate it once, take a sqare root (r is always positive, so no ambiguity here ) and get to a first order equation for $r$. Now you can use a first order RK scheme... $\endgroup$ Jul 19, 2016 at 19:21
  • 2
    $\begingroup$ @scidaddy I want to note that I implemented this dynamical system in matlab and when I integrated it, it too got negative values for $r$ and an essentially $0$ value for $\theta$. So I am pretty sure your equations aren't correct, if this is a result you don't expect. $\endgroup$
    – spektr
    Jul 20, 2016 at 16:49

1 Answer 1


As all the commentators have already pointed out: Something is fishy about your ODE.

  • why is there just $m$ for $\ddot r$ and then $m_1$ and $m_2$ for $\dot \theta$?
  • why is there no equation involving $\ddot \theta$?
  • why is $l$ a parameter where one would actually expect something dynamic like $r\dot \theta^2$ for the angular momentum?
  • why is it $\frac k {mr^2}$, where one would expect $\frac {km} {r^2}$ for a gravitational force?

You obviously want to model celestial mechanics for earth and sun (i.e. classical 2-body problem) in 2D using polar coordinates.

Starting from a Lagrangian (see here for a nice derivation but do not consider the numerical solution) with kinetic energy $\frac {m_1} 2 (\dot r^2 + r^2 \dot \theta^2)$ and potential energy $\frac {-km_1m_2 } r$, we obtain by derivation the following equations of motion (Euler-Lagrange equations):

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

where $m_2$ is the solar mass, i.e. $m_2=1$ in solar mass units.

If I then modify your (sorry but quite ugly) code in a few places ...

#define N 4         /* number of first order equations */
    y[1]= 0.;                                       /* initial velocity */
    y[2] = 0;                                               /* initial angle */
    y[3] = 2*M_PI;                                          /* angular velocity */
        if(j % 100 == 0)
            fprintf(output, "%f\t%f\t%f\t%f\t%f\n", t, y[0], y[1], y[2], y[3]); /* only output every 100th point and output full y */
    switch(i) {
    case 0:
        return(y[1]);         /* derivative of first equation (dr/dt)*/
    case 1:
        return((-k/(y[0]*y[0])) + y[0]*y[3]*y[3]); /* derivative of second equation (d^2r/dt^2) */
    case 2:
        return(y[3]);       /* derivative of first equation (dtheta/dt)*/
    case 3:
        return(-2.0*y[1]*y[3]/y[0]);        /* derivative of second equation (d^2theta /dt^2) */

Then I get a nice time integration with $\theta$ perpetually increasing and $r$ wobbling around its initial value.

enter image description here

In other words, your RK4 code seems to be working.

One final word: If you play around with the time step (bump it up to 1e-3), you'll notice that angular velocity is lost: It does not return to its initial value of $2\pi$ but smaller values. Here is where other (symplectic) integrators work better as they preserve energy and thus are more stable for the these periodic motions, see @cpraveen's comment.


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.