物理学にいくつかの問題があり、コードにいくつかの問題があります。
まず、物理の問題。物理法則が異なる代替宇宙をモデル化していないと仮定すると、ニュートンの万有引力の法則は F=G*m1*m2/(r*r) と述べています。ただし、力はスカラーではなくベクトルであるため、大きさと方向の両方があります。
コードが計算するのforceFuncXは、単に X 軸に平行な力の成分ではなく、実際には大きさです。 forceFuncY同じ欠陥があります。
次に加速度の計算です。物理学によると、これは a=F/m です。質量はスカラーですが、加速度と力はベクトルです。したがって、a_x を計算するには、F_x/m を使用するか、F*cos( a )/m を計算することができます。cos( a ) ( aは 2D 空間での角度CelesitalObject) = dx/r であるため、これを a_x = F*dx/(m*r) にすることができますが、これはほとんどではありますが、得られたものとはまったく異なります。あなたの計算では（除数に r がありません）。
別のアプローチとして を使用するstd::complexこともできますが、このモデルを 3 次元に拡張することを想定して、そのアプローチについては説明しません。
これにより、コードの問題が発生します。CelestialObjectまず、C++ を使用して離散オブジェクトの物理システムのシミュレーションを作成しているため、クラスを定義したことは理にかなっています。あまり意味がないのは、これらのオブジェクトの個々の部分を選択してから 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());
}
最後に実際のベルレット:
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 秒のタイム ステップを使用しますが、これは太陽系には短すぎますが、既知のシステムで妥当な結果が得られるかどうかを確認するこのクイック テストには問題ありません。この場合、惑星地球と金床からなる 2 体システムを選択しました。このコードは、空気抵抗を無視した場合、370m の建物の最上部から金床を落とすことをシミュレートし、約 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 になる必要があることです。これは、このコードが数値解析の観点からあまり適切に設計されていないためではありません。ある数値が大きく、別の数値が小さい場合に浮動小数点数の加算と減算を行うことは、1 つの問題です。もう 1 つは、atan2fa のみを返す関数をfloat使用し、その結果をcos()andで使用することsin()です。三角関数は、実際には可能であれば避けるのが最善です。
最後に、このプログラムは現在 2 つのオブジェクトでしか動作しません。この種のスキームでは 3 分の 1 を追加するのは面倒なので、他のすべての位置と質量を考慮して各オブジェクトにかかる正味の力std::vector<CelestialObject>を最初に計算することにより、a を処理することをお勧めします。それはあなたに任せますが、これで少なくとも正しい方向へのスタートが切れるはずです。