太陽能系統仿真項目（速度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>，通過考慮所有其他人的位置和質量來處理。我會留給你的，但這至少應該給你一個正確的方向。