0

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!

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • What is `double k1x = dt(x,y,t);`? Maybe you mean `double k1x = dt* dx(x,y,t);` – merovingian Dec 14 '22 at 10:42
  • 3
    There are few typos in your code but the main problem is the `(1 / 6)` part of statements like `double kx = (1 / 6) * (k1x + 2 * k2x + 2 * k3x + k);` (and others). That is **integer arithmetic** and will evaluate to zero. Try `1.0/6` to use floating-point arithmetic. – Adrian Mole Dec 14 '22 at 10:45
  • By the way, Runge Kutta is a very inefficient way to do this – merovingian Dec 14 '22 at 10:54
  • @AdrianMole Thanks for your help. Changing 1/6 -> 1.0/6 and 2->2.0 has helped in the sense that the code now outputs changing variables, however the values are incorrect. For example, the initial x value = 1*AU, should be the maximum possible value for x, however the code outputs forever increasing values for x and y. Any ideas? –  Dec 14 '22 at 11:01
  • @merovingian Yes I've seen a lot of people say that. Only reason I'm using it is because it's what I've been asked to do. –  Dec 14 '22 at 11:02
  • Use `atan2(y,x)` if you want to get angles on the full circle. Note also that `sin(theta)=y/r`, so the formula can be simplified. There might be missing a sign if you want the velocity to be tangent to the circle of radius `r`. – Lutz Lehmann Dec 15 '22 at 10:19
  • See https://stackoverflow.com/questions/53645649/cannot-get-rk4-to-solve-for-position-of-orbiting-body-in-python for a discussion of similar misunderstandings of the algorithm in a python code. – Lutz Lehmann Dec 15 '22 at 10:48
  • Possible typo `double kx = (1/6) * (k1x + 2*k2x + 2*k3x + k);` the last term should be `k4x` right? – John Alexiou Jan 06 '23 at 18:20
  • is `3 / 2` in `pow(x * x + y * y, 3 / 2))` integer division in C? Change it to `pow(x * x + y * y, 3 / 2.0))` – John Alexiou Jan 06 '23 at 18:30
  • `intiial_vx` is a typo. it should be `initial_vx`. Please fix the code typos and make sure the code compiles before posting. – John Alexiou Jan 06 '23 at 18:37

4 Answers4

1

  • 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
user58697
  • 7,808
  • 1
  • 14
  • 28
0

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

merovingian
  • 515
  • 2
  • 9
  • Ah yeah I've just changed the k1vy's. Although that doesn't affect the x position values. Both the x and y position values are still increasing above what the maximum value should be, and the velocities are all above the speed of light. –  Dec 14 '22 at 11:47
  • Because you are not doing enough checks – merovingian Dec 14 '22 at 12:07
  • Okay so to check whether the planet's position is within the bounds of the simulation I would just use an if statement, however I can't see what I would actually have the if statement do for those values that exceed the minimum or maximum limits, could you elaborate a little more on this? Also, for your point on the planet's distance from the star changing with time, is this not taken care of by the x and y co-ordinates of the planet changing over time? Relative to the star which is, for this purpose, fixed at the origin. –  Dec 14 '22 at 12:28
0

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
0

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.

John Alexiou
  • 28,472
  • 11
  • 77
  • 133