4th Order Runge Kutta: Integration of Differential Equations for Planetary Orbit

Question

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 */
{
    t=j*dist;
    runge4(t, y, dist);

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

fclose(output);
}

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.

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)*/
        break;;
    case 1:
        return((-k/(y[0]*y[0])) + y[0]*y[3]*y[3]); /* derivative of second equation (d^2r/dt^2) */
        break;;
    case 2:
        return(y[3]);       /* derivative of first equation (dtheta/dt)*/
        break;;
    case 3:
        return(-2.0*y[1]*y[3]/y[0]);        /* derivative of second equation (d^2theta /dt^2) */
        break;;
    }

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

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.