C# 中的 N 体模拟

我正在尝试使用 Runge Kutta 4 或 Velocity Verlet 集成算法在 C# 中实现 N 体模拟。

在我转向更多数量的粒子之前，我想通过模拟地球绕太阳的轨道来测试模拟，但是，由于某种原因，我得到的不是椭圆轨道，而是一个奇怪的螺旋。

我无法解决这个问题，因为我使用相同的算法对太阳系进行了更简单的模拟，其中太阳固定在位置并且一切都运行良好。积分器工作得很好，因为无论我使用哪一个，我都会得到螺旋。

任何帮助，将不胜感激。

这是代码：

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();
    }
}

这是地球轨道的输出数据图：

您的模型计算两个粒子之间的重力两次：对于第一个粒子，力基于其原始坐标，对于第二个粒子，力基于第一个粒子的更新位置。这明显违反了牛顿第三定律。在任何更新之前，您必须预先计算所有力。