太阳能系统仿真项目（速度verlet的帮助）

Question

对于我的建模和模拟类的项目，我想模拟太阳系。我从一个明星（太阳）和一个行星（地球）开始，但我已经遇到了一些问题。我花了一些时间来回顾并学习不同的公式，以及如何模拟星球轨道受星星和周围物体影响的方式。我想使用速度verlet，并最终研究n体问题。我有很多问题与我的速度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

，你的积分接受标量和你的问题是关于2维系统让我觉得，你是为每个组件调用积分两次，一次是事实。这完全不起作用，因为你的系统将在阶段空间中采取不切实际的动作。积分器可与载体量：

X（T + DT）= X（T）+ V（T）dt的+（1/2）甲（T）dt的

V（T + DT）= V（T）+（1/2）（甲（T）+ 甲（T + DT））dt的

这里X（t）是一个列向量，它由所有颗粒的坐标的 - 这是该系统的相空间的配置的子空间。 V（t）是所有粒子的速度的列向量，在技术上表示动量子空间。这同样适用于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维向量并包含n组件（例如，在您的2d案例中为2）的类。它还重载+和+=运营商实施矢量（逐个分量）另外，以及*和/由一个标量来实现乘法/除法。这只是为了说明这种情况，可以用每个位置/速度矢量分量的内部循环代替。

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

Answer 2

物理上存在一些问题，代码存在一些问题。

首先，物理问题。假设我们不是建模物理定律不同的另一个宇宙，牛顿的万有引力定律说F = G * m1 * m2 /（r * r）。然而，力量是一个矢量，而不是一个标量，所以它既有量级又有方向。

代码在forceFuncX中计算的实际值是幅度，而不仅仅是平行于X轴的力的分量。 forceFuncY也有同样的缺陷。

接下来是计算加速度。物理学说这是a = F/m。质量是一个标量，但加速度和力量都是矢量。因此，要计算a_x，我们可以使用F_x/m或者我们可以计算F * cos（a）/ m。因为cos（a）（其中a是从一个CelesitalObject到2D空间中另一个的角度）= dx/r，我们可以使这个a_x = F * dx /（m * r）你的计算中有什么（它在除数中缺少r）。

另一种方法是使用std::complex，但我不会在假设您希望将此模型扩展到三维的情况下显示该方法。

这给我们带来了代码的问题。首先，由于您使用C++并编写了离散对象物理系统的模拟，所以您定义了一个CelestialObject类是有道理的。更不重要的是，你的函数是通过挑选这些对象的各个部分然后调用C风格的函数来调用的。我们可以更好地使用这些对象来改进代码。首先，因为你还没有公布一个，这里是基于我从你的代码推断在接口上CelestialObject类：

class CelestialObject 
{ 
public: 
    CelestialObject(std::string name, float mass, float X, float Y, 
     float VX, float VY) 
      : myname(name), m(mass), x(X), y(Y), vx(VX), vy(VY) {} 
    void setPosition(float X, float Y) { x=X; y=Y; } 
    void setVelocity(float VX, float VY) { vx=VX; vy=VY; } 
    float getMass() const { return m; } 
    float getX() const { return x; } 
    float getY() const { return y; } 
    float getVX() const { return vx; } 
    float getVY() const { return vy; } 
    friend std::ostream& operator<<(std::ostream& out, 
            const CelestialObject& obj) { 
     return out << obj.myname << '\t' << obj.x << '\t' << obj.y 
        << '\t' << obj.vx << '\t' << obj.vy << std::endl; 
    } 
private: 
    std::string myname; 
    float m, x, y; 
    float vx, vy; 
}; 

接下来，一些辅助功能：

// returns square of distance between objects 
float distance_sq(const CelestialObject &a, const CelestialObject &b) 
{ 
    // distance squared is (dy^2 + dx^2) 
    return pow(a.getY()-b.getY(),2) + pow(a.getX()-b.getX(),2); 
} 

// returns magnitude of the force between the objects 
float force(const CelestialObject &a, const CelestialObject &b) 
{ 
    // F=(G * m1 * m1)/(r^2) in the direction a->b and b->a 
    return G*a.getMass()*b.getMass()/distance_sq(a, b); 
} 

// returns the angle from a to b 
float angle(const CelestialObject &a, const CelestialObject &b) 
{ 
    return atan2f(b.getY()-a.getY(),b.getX()-a.getX()); 
} 

最后实际verlet的：

void updatePosition(CelestialObject &a, CelestialObject &b) 
{ 
    float F = force(a,b); 
    float theta = angle(a,b); 
    float accela = F/a.getMass(); 
    float accelb = -F/b.getMass(); 

    // now that we have the acceleration of both objects, update positions 
    // x = x +v *dt + a*dt*dt/2 
    // = x + dt * (v + a*dt/2) 
    a.setPosition(
    a.getX() + dt * (a.getVX() + accela*cos(theta)*dt/2), 
    a.getY() + dt * (a.getVY() + accela*sin(theta)*dt/2) 
    ); 
    b.setPosition(
    b.getX() + dt * (b.getVX() + accelb*cos(theta)*dt/2), 
    b.getY() + dt * (b.getVY() + accelb*sin(theta)*dt/2) 
    ); 
    // get new acceleration a' 
    F = force(a,b); 
    float thetap = angle(a,b); 
    float accelap = F/a.getMass(); 
    float accelbp = -F/b.getMass(); 
    // and update velocities 
    // v = v + (a + a')*dt/2 
    a.setVelocity(
    a.getVX() + (accela*cos(theta) + accelap*cos(thetap))*dt/2, 
    a.getVY() + (accela*sin(theta) + accelap*sin(thetap))*dt/2 
    ); 
    b.setVelocity(
    b.getVX() + (accelb*cos(theta) + accelbp*cos(thetap))*dt/2, 
    b.getVY() + (accelb*sin(theta) + accelbp*sin(thetap))*dt/2 
    ); 
} 

最后是一些简单的测试代码。

#include <string> 
#include <iostream> 
#include <vector> 
#include <cmath> 

const float G(6.67e-11); // N*(m/kg)^2 
const float dt(0.1);  // s 
// all of the other code goes here... 
int main() 
{ 
    CelestialObject anvil("anvil", 70, 370, 0, 0, 0); 
    CelestialObject earth("earth", 5.97e+24, -6.378e6, 0, 0, 0); 
    std::cout << "Initial values:\n" << earth << anvil; 
    std::cout << "Dropping an anvil from the top of a 370m building...\n" 
       "It should hit the ground in about 8.7 seconds.\n"; 
    int t; 
    for (t=0; anvil.getX() > 0; ++t) { 
    std::cout << dt*t << '\t' << anvil; 
    updatePosition(anvil, earth); 
    } 
    std::cout << "Final values at t = " << dt*t << " seconds:\n" 
       << earth << anvil; 
    return 0; 
} 

测试代码使用的0.1秒时间步远远太短，你们的太阳系，但罚款这种快速测试，就是看如果我们得到了一个已知的系统一个合理的结果。在这种情况下，我选择了由地球和铁砧组成的双体系统。该代码模拟从370米建筑物的顶部放下一个砧座，如果我们忽略空气阻力，我们可以轻松计算出大约8.7秒的距离。为了保持简单的坐标，我选择将原点（0,0）放置在地球表面，并考虑建筑物顶部（370,0）。当代码被编译和运行，它会产生如下：

Initial values: 
earth -6.378e+06 0 0 0 
anvil 370 0 0 0 
Dropping an anvil from the top of a 370m building... 
It should hit the ground in about 8.7 seconds. 
0 anvil 370 0 0 0 
0.1 anvil 369.951 -4.27834e-09 -0.97877 -8.55668e-08 
0.2 anvil 369.804 -1.71134e-08 -1.95754 -1.71134e-07 
0.3 anvil 369.56 -3.85051e-08 -2.93631 -2.567e-07 
    ... 
8.3 anvil 32.8567 -2.9474e-05 -81.2408 -7.1023e-06 
8.4 anvil 24.6837 -3.01885e-05 -82.2197 -7.18787e-06 
8.5 anvil 16.4127 -3.09116e-05 -83.1985 -7.27345e-06 
8.6 anvil 8.04394 -3.16432e-05 -84.1774 -7.35902e-06 
Final values at t = 8.7 seconds: 
earth -6.378e+06 3.79705e-28 9.98483e-22 8.72901e-29 
anvil -0.422744 -3.23834e-05 -85.1563 -7.4446e-06 

正如你所看到的，这似乎工作，但也存在问题。第一个问题是由于对象只能沿着X轴移动，所以所有的Y分量都应该是0.它们并不是因为从数值分析的角度来看，这些代码设计得不是很好。当一个数字很大而另一个很小时，对浮点数进行加法和减法是一个问题。另一种是使用atan2f函数，该函数仅返回float，然后使用cos()和sin()中的结果。如果可能，实际上最好避免三角函数。

最后，这个程序目前只适用于两个对象。添加三分之一将是痛苦的这种方案，所以更好的设计将通过首先计算每个对象上的净力std::vector<CelestialObject>，通过考虑所有其他人的位置和质量来处理。我会留给你的，但这至少应该给你一个正确的方向。