c++ - 太阳系模拟项目(velocity verlet 帮助)

Question

对于我的建模和仿真类(class)项目，我想模拟一个太阳系。我从一颗恒星(太阳)和一颗行星(地球)开始，但我已经遇到了一些问题。我现在花了一些时间来回顾和学习不同的公式和方法来模拟行星的轨道将如何受到恒星和周围物体的影响。我想使用 velocity verlet 并最终研究 n-body 问题。我的 velocity verlet 函数有很多问题。有时它的行为就好像它在正常地绕地球轨道运行，然后它“曲速驱动”地球进入某个随机位置。我还注意到我从来没有得到“负”加速度，所以我的 x 和 y 坐标。总是在增加，所以我看不出地球应该如何环绕太阳。任何帮助是极大的赞赏。 AGK::Prints 我刚刚看到不同的变量是如何变化的。

double velocityVerlet(float positionCalc, double position2, 
                      float &velocity, double massCalc, double mass2)
//positionCalc is the position being updated, position 2 is position of 
// other object, same with mass
{
    float force = forceFunc(positionCalc, position2, massCalc, mass2);
    agk::PrintC("Force is: ");
    agk::Print(force);
    float acceleration = accelerationFunc(force,massCalc);
    agk::PrintC("Accel is: ");
    agk::Print(acceleration);`;

    double newAccel = 0;

    positionCalc = positionCalc + velocity*dt + 
                   (.5*acceleration)*pow(dt,2); //calculates new position
    agk::PrintC("New Position is: ");
    agk::Print(positionCalc);
    force = forceFunc(positionCalc,position2,massCalc,mass2);
    newAccel = accelerationFunc(force, massCalc);

    velocity = velocity + .5*(acceleration + newAccel)*dt; //new velocity
    agk::PrintC("Velocity is: ");
    agk::Print(velocity);

    return positionCalc;
}

Answer 1

您的积分器接受标量以及您的问题是关于二维系统的事实让我认为您调用积分器两次，每个组件一次。这根本行不通，因为您的系统将在相空间中采取不切实际的移动。积分器使用 vector :

X(t+dt) = X(t) + V(t) dt + (1/2) A(t) dt²

V(t+dt) = V(t) + (1/2)(A(t) + A(t+dt)) dt

此处X(t) 是一个列 vector ，由所有 粒子的坐标组成 - 这是系统相空间的配置子空间。 V(t) 是所有粒子速度的列 vector ，技术上代表动量子空间。这同样适用于A(t)。这些必须同时更新，而不是分别更新。

整个速度 Verlet 过程在不依赖于速度的力场(例如经典重力)的代码中翻译如下:

Vector forces[num_particles];

// Compute initial forces
forces = computeForces(positions);

for (int ts = 0; ts < num_timesteps; ts++)
{
   // Update positions and half-update velocities
   for (int i = 0; i < num_particles; i++)
   {
      positions[i] += velocities[i]*dt + 0.5*(forces[i] / m[i]) * dt*dt;
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }

   // Compute new forces and half-update velocities
   forces = computeForces(positions);

   for (int i = 0; i < num_particles; i++)
   {
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }
}

请注意，在下一轮力评估之前，首先更新所有位置。此外，每次迭代只需要评估一次力，因为在第二次更新速度期间位置不会改变。在上面的示例代码中，Vector 是一个实现 n 维 vector 并包含 n 组件(例如 2d-case 中的 2)的类。它还重载了 + 和 += 运算符以实现 vector (分量)加法，以及 * 和 / 通过标量实现乘法/除法。这只是为了说明情况，可以用每个位置/速度 vector 的分量上的内部循环代替。

正确选择时间步长非常重要。太小的时间步长会显着减慢模拟速度。太大的时间步将导致不可能的物理，例如。跳跃的行星。