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

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!

CodePudding user response：

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

CodePudding user response：

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