Why is RK4 not updating values for my simulation of the Earth's orbit?

Question

I am trying to use a Runge Kutta method to simulate the motion of the Earth around the Sun in C, I cannot see why but my code does not update the values for position or velocity and just keeps the initial values. The code I have written is:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define dt 86400 // 1 day in seconds //

const double G = 6.67e-11;
const double au = 1.496e11;
const double M = 1.99e30;

double vx(double x, double y);
double vy(double x, double y);
double dx(double x, double y, double t);
double dy(double x, double y, double t);
double dvx(double x, double y, double t);
double dvy(double x, double y, double t);

int main(){

  double initial_x = au;
  double initial_y = 0;
  double intiial_vx = vx(initial_x, initial_y);
  double initial_vy = vy(initial_x, initial_y);
  double t = 0;

  double x = initial_x;
  double y = initial_y;
  double vx = initial_vx;
  double vy = initial_vy;

  for(int i=0;i<365;i++){
       double k1x = dt(x,y,t);  
       double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);
       double k3x = dt * dx(x + k2x/2, y + k2x/2, t + dt/2);
       double k4x = dt * dx(x + k3x, y + k3x, t + dt);
       double kx = (1/6) * (k1x + 2*k2x + 2*k3x + k);

       double k1y = dt * dy(x,y,t);
       double k2y = dt * dy(x + k1y/2, y + k1y/2, t + dt/2);
       double k3y = dt * dy(x + k2y/2, y + k2y/2, t + dt/2);
       double k4y = dt * dy(x + k3y, y + k3y, t + dt);
       double ky = (1/6) * (k1y + 2*k2y + 2*k3y + k4y);

       double k1vx = dt * dvx(x,y,t);
       double k2vx = dt * dvx(x+k1vx/2, y+k1vx/2, t + dt/2);
       double k3vx = dt * dvx(x+k2vx/2, y+k2vx/2, t + dt/2);
       double k4vx = dt * dvx(x+k3vx, y+k3vx, t+dt);
       double kvx = (1/6) * (k1vx + 2*k2vx + 2*k3vx + k4vx);

       double k1vy = dt * dvx(x,y,t);
       double k2vy = dt * dvx(x+k1vy/2, y+k1vy/2, t + dt/2);
       double k3vy = dt * dvx(x+k2vy/2, y+k2vy/2, t + dt/2);
       double k4vy = dt * dvx(x+k3vy, y+k3vy, t+dt);
       double kvy = (1/6) * (k1vy + 2*k2vy + 2*k3vy + k4vy);

       x = x + kx;
       y = y + ky;
       vx = vx + kvx;
       vy = vy + kvy;

       printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx, vy);
  }

  return 0;
}



// Function for the x velocity of a planet//
double vx(double x, double y)
{
    double theta = atan(y/x);
    double xVel = sqrt((G*M) / (sqrt(x*x + y*y))) * sin(theta);
    return xVel;
}

// Function for the y velocity of a planet //
double vy(double x, double y) 
{
    double theta = atan(y/x);
    double yVel = sqrt((G*M) / (sqrt(x*x + y*y))) * cos(theta);
    return yVel;
}

// Function for dx //
double dx(double x, double y, double t)
{
    double xVel = vx(x,y);
    double dX = xVel*t;
    return dX;
}

// Function for dy //
double dy(double x, double y, double t)
{
    double yVel = vy(x,y);
    double dY = yVel*t;
    return dY;
}

// Function for dvx //
double dvx(double x, double y, double t)
{
    double dVX = ((-G*M*x) / pow(x*x+y*y, 3/2)) * t;
    return dVX;
}

// Function for dvy //
double dvy(double x, double y, double t)
{
    double dVY = (((-G*M*x) / pow(x*x+y*y, 3/2))) * t;
    return dVY;
}

I haven't yet added functions for the Runge-Kutta simply because I can't get it to work. Thanks for any help!

Answer 1

Because you are calculating the k1vy's in terms of dvx. It should be dvy.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int dt = 86400; // 1 day in seconds //
const double G = 6.67e-11;
const double au = 1.496e11;
const double M = 1.99e30;

double vx(double x, double y);
double vy(double x, double y);
double dx(double x, double y, double t);
double dy(double x, double y, double t);
double dvx(double x, double y, double t);
double dvy(double x, double y, double t);

int main(){

  double initial_x = au;
  double initial_y = 0;
  double initial_vx = vx(initial_x, initial_y);
  double initial_vy = vy(initial_x, initial_y);
  double t = 0;

  double x = initial_x;
  double y = initial_y;
  double vx = initial_vx;
  double vy = initial_vy;

  for(int i=0;i<365;i++){
       double k1x = dt* dx(x,y,t);  
       double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);
       double k3x = dt * dx(x + k2x/2, y + k2x/2, t + dt/2);
       double k4x = dt * dx(x + k3x, y + k3x, t + dt);
       double kx = (1.0/6) * (k1x + 2*k2x + 2*k3x + k4x);

       double k1y = dt * dy(x,y,t);
       double k2y = dt * dy(x + k1y/2, y + k1y/2, t + dt/2);
       double k3y = dt * dy(x + k2y/2, y + k2y/2, t + dt/2);
       double k4y = dt * dy(x + k3y, y + k3y, t + dt);
       double ky = (1.0/6) * (k1y + 2*k2y + 2*k3y + k4y);

       double k1vx = dt * dvx(x,y,t);
       double k2vx = dt * dvx(x+k1vx/2, y+k1vx/2, t + dt/2);
       double k3vx = dt * dvx(x+k2vx/2, y+k2vx/2, t + dt/2);
       double k4vx = dt * dvx(x+k3vx, y+k3vx, t+dt);
       double kvx = (1.0/6) * (k1vx + 2*k2vx + 2*k3vx + k4vx);

       double k1vy = dt * dvy(x,y,t);
       double k2vy = dt * dvy(x+k1vy/2, y+k1vy/2, t + dt/2);
       double k3vy = dt * dvy(x+k2vy/2, y+k2vy/2, t + dt/2);
       double k4vy = dt * dvy(x+k3vy, y+k3vy, t+dt);
       double kvy = (1.0/6) * (k1vy + 2*k2vy + 2*k3vy + k4vy);

       x = x + kx;
       y = y + ky;
       vx = vx + kvx;
       vy = vy + kvy;

       printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx, vy);
  }

  return 0;
}



// Function for the x velocity of a planet//
double vx(double x, double y)
{
    double theta = atan(y/x);
    double xVel = sqrt((G*M) / (sqrt(x*x + y*y))) * sin(theta);
    return xVel;
}

// Function for the y velocity of a planet //
double vy(double x, double y) 
{
    double theta = atan(y/x);
    double yVel = sqrt((G*M) / (sqrt(x*x + y*y))) * cos(theta);
    return yVel;
}

// Function for dx //
double dx(double x, double y, double t)
{
    double xVel = vx(x,y);
    double dX = xVel*t;
    return dX;
}

// Function for dy //
double dy(double x, double y, double t)
{
    double yVel = vy(x,y);
    double dY = yVel*t;
    return dY;
}

// Function for dvx //
double dvx(double x, double y, double t)
{
    double dVX = ((-G*M*x) / pow(x*x+y*y, 3/2)) * t;
    return dVX;
}

// Function for dvy //
double dvy(double x, double y, double t)
{
    double dVY = (((-G*M*x) / pow(x*x+y*y, 3/2))) * t;
    return dVY;
} 

Some other issues that you should check:

Take into account the fact that the planet's distance from the star will change over time

Your code does not check whether the planet's position is within the bounds of the simulation

Answer 2

• Your physics is wrong. The equations of motion do not explicitly depend on time. E.g, why dvx

  double dVX = ((-G*M*x) / pow(x*x+y*y, 3/2)) * t;
  

  should be different at different time? It shall only depend on x, y.

  Do not multiply by t.

• The Runge-Kutta implementation is wrong too. In your

   double k1x = dt(x,y,t);  
   double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);
   double k3x = dt * dx(x + k2x/2, y + k2x/2, t + dt/2);
   double k4x = dt * dx(x + k3x, y + k3x, t + dt);
   double kx = (1/6) * (k1x + 2*k2x + 2*k3x + k);
  

  you multiply by dt far too much. Keep in mind the dt is small. Every multiplication by a small value diminish the corrective power of subsequent k*x. You should do

   double k1x = dx(x, y, t);
   double k2x = dx(x + dt * k1x/x, y + dt * k1x/2, t + dt/2);
   ....
   double kx = (dt/6) * (k1x + 2*k2x + 2*k3x + k4x);
  
   // etc
  

Answer 3

Let's also look a little deeper at

 double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);

It should be obvious from the variable names that y + k1x/2 is wrong. It should be

y + k1y/2

Which then of course requires computing k1y before computing k2x.

You can only start computing the next stage of a Runge-Kutta method when you have completed the current stage in all components of the state vector.

---

By the way, the correct system is

x' = vx
y' = vy
vx' = ax(x,y)
vy' = ay(x,y)

so that the correct line for k2x would read

 double k2x = dt * (vx + k1vx/2);

which of course again requires that k1vx be computed beforehand.

Answer 4

There is a big mistake here by updating x first, then y, then vx and finally vy because you have to make sure that values on the same stage are used for the evaluation of dvx and dvy.

I recommend re-arranging the terms to update x, y, vx and vy on each stage side by side.

There are numerous other mistakes and problems with your code. I'd rather show the correct way of doing this, rather than try to deal with all the bugs.

const double year = 365.0; // 1 year in days //
const double day = 86400.0; // 1 day in seconds //
const double G = 6.67e-11;
const double au = 1.496e11;
const double M = 1.99e30;

int main()
{
    // Initial Conditions
    double t = 0;
    double x = au;
    double y = 0;
    double vx = init_vx(x, y);
    double vy = init_vy(x, y);

    const int steps = 365;
    double dt = day * (year / steps);

    // Initial Values
    printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx, vy);

    // Simulation Loop
    for (int i = 0; i < steps; i++)
    {
        // stage 1
        double k1x = dx(vx, vy);
        double k1y = dy(vx, vy);
        double k1vx = dvx(x, y);
        double k1vy = dvy(x, y);

        //stage 2
        double k2x = dx(vx + (dt / 2) * k1vx, vy + (dt / 2) * k1vy);
        double k2y = dy(vx + (dt / 2) * k1vx, vy + (dt / 2) * k1vy);
        double k2vx = dvx(x + (dt / 2) * k1x, y + (dt / 2) * k1y);
        double k2vy = dvy(x + (dt / 2) * k1x, y + (dt / 2) * k1y);

        //stage 3
        double k3x = dx(vx + (dt / 2) * k2vx, vy + (dt / 2) * k2vy);
        double k3y = dy(vx + (dt / 2) * k2vx, vy + (dt / 2) * k2vy);
        double k3vx = dvx(x + (dt / 2) * k2x, y + (dt / 2) * k2y);
        double k3vy = dvy(x + (dt / 2) * k2x, y + (dt / 2) * k2y);

        //stage 4
        double k4x = dx(vx + (dt) * k3vx, vy + (dt) * k3vy);
        double k4y = dy(vx + (dt) * k3vx, vy + (dt) * k3vy);
        double k4vx = dvx(x + (dt) * k3x, y + (dt) * k3y);
        double k4vy = dvy(x + (dt) * k3x, y + (dt) * k3y);

        // slope calculation
        double kx = (k1x + (2.0 * k2x) + (2.0 * k3x) + k4x) / 6.0;
        double ky = (k1y + (2.0 * k2y) + (2.0 * k3y) + k4y) / 6.0;
        double kvx = (k1vx + (2.0 * k2vx) + (2.0 * k3vx) + k4vx) / 6.0;
        double kvy = (k1vy + (2.0 * k2vy) + (2.0 * k3vy) + k4vy) / 6.0;

        // integration step
        t = t + dt;
        x = x + dt * kx;
        y = y + dt * ky;
        vx = vx + dt * kvx;
        vy = vy + dt * kvy;

        // Simulation results
        printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx, vy);
    }

     return 0;
}

// Function for the initial x velocity of a planet//
double init_vx(double x, double y)
{
    double theta = atan2(y, x);
    double xVel = -sqrt((G * M) / (sqrt(x * x + y * y))) * sin(theta);
    return xVel;
}

// Function for the initial y velocity of a planet //
double init_vy(double x, double y)
{
    double theta = atan2(y, x);
    double yVel = +sqrt((G * M) / (sqrt(x * x + y * y))) * cos(theta);
    return yVel;
}

// Function for dx //
double dx(double vx, double vy)
{
    double xVel = vx;
    return xVel;
}

// Function for dy //
double dy(double vx, double vy)
{
    double yVel = vy;
    return yVel;
}

// Function for dvx //
double dvx(double x, double y)
{
    double dVX = ((-G * M * x) / pow(x * x + y * y, 3 / 2.0));
    return dVX;
}

// Function for dvy //
double dvy(double x, double y)
{
    double dVY = ((-G * M * y / pow(x * x + y * y, 3 / 2.0)));
    return dVY;
}

The code above has the flexibility to simulate 1 year in steps=365 or some other value. I would report out t/day values which is the time in days.

The big issues where

• Incorrect handling of time in the d() functions.
• Incorrect handling of intermediate values for each sub-step in RK4.
• Careless use of integer division, 3/2 = 1 and not 1.5
• Several typos all over the place, for example in dvy() it should use the y coordinate in the numerator.
• Bug in initial conditions when (x,y) not int he first quadrant. Should use the atan2() function for finding the angle, and the x component should have a negative sign in front of it.

_{There are some optimizations that can be done to my code because some values are evaluated two times, but I think it won't make much of a difference as the floating point math is rather fast nowadays. I left them out for clarity really.}