N body simulation in C#

Question

I'm trying to implement an N body simulation in C# using either Runge Kutta 4 or Velocity Verlet integration algorithms.

Before I move to a bigger number of particles, I wanted to test the simulation by modeling the earth's orbit around the sun, however, instead of the elliptical orbit, I get a weird spiral for some reason.

I can't figure out the problem since I made a simpler simulation of the solar system using the same algorithms where the sun was fixed in position and everything worked perfectly. The integrators work perfectly because it doesn't matter which one I use, I get the spiral with both.

Any help would be appreciated.

Here's the code:

class NBODY
{
    public static double G = 4 * Math.PI * Math.PI;

    class Particle
    {
        public double[] r;       // position vector
        public double[] v;       // velocity vector
        public double mass;    

        //constructor
        public Particle() {}
        public Particle(double x, double y, double z, double vx, double vy, double vz, double m)
        {
            this.r = new double[3];
            this.v = new double[3];

            this.r[0] = x;
            this.r[1] = y;
            this.r[2] = z;
            this.v[0] = vx;
            this.v[1] = vy;
            this.v[2] = vz;
            this.mass = m;
        }

        public void Update(Particle[] particles, double t, double h, int particleNumber)
        {
            RungeKutta4(particles, t, h, particleNumber);
        }

        private double acc(double r, Particle[] particles, int particleNumber, double[] r_temp, int l)
        {

            // dv/dt = f(x) = -G * m_i * (x - x_i) / [(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2]^(3/2)

            double sum = 0;

            switch (l)
            {
                case 0:
                    for (int i = 0; i < particles.Length; i++)
                        if (i != particleNumber)
                            sum += particles[i].mass * (r - particles[i].r[l]) / Math.Pow( Math.Pow(r - particles[i].r[l], 2)
                                + Math.Pow(r_temp[1] - particles[i].r[1], 2) + Math.Pow(r_temp[2] - particles[i].r[2], 2), 1.5);
                    break;

                case 1:
                    for (int i = 0; i < particles.Length; i++)
                        if (i != particleNumber)
                            sum += particles[i].mass * (r - particles[i].r[l]) / Math.Pow(Math.Pow(r - particles[i].r[l], 2)
                                + Math.Pow(r_temp[0] - particles[i].r[0], 2) + Math.Pow(r_temp[2] - particles[i].r[2], 2), 1.5);
                    break;

                case 2:
                    for (int i = 0; i < particles.Length; i++)
                        if (i != particleNumber)
                            sum += particles[i].mass * (r - particles[i].r[l]) / Math.Pow(Math.Pow(r - particles[i].r[l], 2)
                                + Math.Pow(r_temp[0] - particles[i].r[0], 2) + Math.Pow(r_temp[1] - particles[i].r[1], 2), 1.5);
                    break;
            }

            return -G * sum;
        }

        private void RungeKutta4(Particle[] particles, double t, double h, int particleNumber)
        {
            //current position of the particle is saved in a vector
            double[] r_temp = new double[3];

            for (int j = 0; j < 3; j++)
                r_temp[j] = this.r[j];

            //loop going over all the coordinates and updating each using RK4 algorithm
            for (int l = 0; l < 3; l++)
            {   
                double[,] k = new double[4, 2];

                k[0, 0] = this.v[l];                                                                //k1_r
                k[0, 1] = acc(this.r[l], particles, particleNumber, r_temp, l);                     //k1_v

                k[1, 0] = this.v[l] + k[0, 1] * 0.5 * h;                                            //k2_r
                k[1, 1] = acc(this.r[l] + k[0, 0] * 0.5 * h, particles, particleNumber, r_temp, l); //k2_v

                k[2, 0] = this.v[l] + k[1, 1] * 0.5 * h;                                            //k3_r
                k[2, 1] = acc(this.r[l] + k[1, 0] * 0.5 * h, particles, particleNumber, r_temp, l); //k3_v

                k[3, 0] = this.v[l] + k[2, 1] * h;                                                  //k4_r
                k[3, 1] = acc(this.r[l] + k[2, 0] * h, particles, particleNumber, r_temp, l);       //k4_v

                this.r[l] += (h / 6.0) * (k[0, 0] + 2 * k[1, 0] + 2 * k[2, 0] + k[3, 0]);
                this.v[l] += (h / 6.0) * (k[0, 1] + 2 * k[1, 1] + 2 * k[2, 1] + k[3, 1]);   
            }
        }

        /*
            Velocity Verlet algorithm:
            1. Calculate y(t+h) = y(t) + v(t)h + 0.5a(t)h*h
            2. Derive a(t+h) from dv/dt = -y using y(t+h)
            3. Calculate v(t+h) = v(t) + 0.5*(a(t) + a(t+h))*h
        */
        private void VelocityVerlet(Particle[] particles, double t, double h, int particleNumber)
        {

            double[] r_temp = new double[3];

            for (int j = 0; j < 3; j++)
                r_temp[j] = this.r[j];

            //loop going over all the coordinates and updating each using RK4 algorithm
            for (int l = 0; l < 3; l++)
            {
                //position
                this.r[l] += h * this.v[l] + 0.5 * h * h * acc(this.r[l], particles, particleNumber, r_temp, l);
                //velocity
                this.v[l] += 0.5 * h * (acc(r_temp[l], particles, particleNumber, r_temp,l)
                    + acc(this.r[l], particles, particleNumber, r_temp,l));
            }
        }     


    }

    static void Main(string[] args)
    {
        //output file
        TextWriter output = new StreamWriter("ispis.txt");

        // declarations of variables
        Particle[] particles = new Particle[2];
        particles[0] = new Particle(0, 0, 0, 0, 0, 0, 1);                  //sun
        particles[1] = new Particle(1, 0, 0, 0, 6.28, 0,  3.003467E-06);   //earth


        int N = 200;
        double h, t, tmax;
        double[,,] x = new double[particles.Length, N, 3];  //output


        //  setting initial values, step size and max time tmax
        h = 0.01;   // the step size in years
        tmax = h * N;

        // initial time        
        t = 0;

        int i = 0;

        while (t <= tmax) {

            //updates position of all particles
            for (int z = 1; z < particles.Length; z++)
                particles[z].Update(particles, t, h, z);

            //saves the position for output
            for (int j = 1; j < particles.Length ; j++)
                for (int z = 0; z < 3; z++ )
                    x[j,i,z] = particles[j].r[z];

            t += h;
            i++;
        }

        //output to file
        for (int k = 0; k < particles.Length; k++ )
        {
            for (int f = 0; f < 3; f++)

            {
                for (int l = 0; l < N; l++)
                    output.Write(string.Format("{0,-15:0.########},", x[k,l,f]));
                output.Write(string.Format("\n\n"));
            }
            output.Write(string.Format("\n\n\n\n"));
        }

        output.Close();
    }
}

And here's the plot of the output data for earth's orbit:

Answer 1

Your model calculates the gravity force between two particles twice: for the first particle the force is based on their original coordinates, and for the second particle it is based on an updated position of the first one. This is a clear violation of the Newton's 3rd law. You must precompute all the forces before any update.