物理有一些问题，代码有一些问题。
首先，物理问题。假设我们不是在模拟一个物理定律不同的替代宇宙，牛顿的万有引力定律表明 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）= dx/r，我们可以使这个 a_x = F*dx/(m*r) 几乎但不完全是你得到的在您的计算中（除数中缺少 r）。
另一种方法是使用std::complex，但我不会在假设您可能希望将此模型扩展到三个维度的情况下展示该方法。
这给我们带来了代码问题。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());
}
最后是实际的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.1s 的时间步长，这对于你的太阳系来说太短了，但对于这个快速测试来说很好，看看我们是否能得到一个已知系统的合理结果。在这种情况下，我选择了由地球和铁砧组成的二体系统。这段代码模拟了从 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。它们不是因为从数值分析的角度来看，这段代码设计得不是很好。当一个数字很大而另一个数字很小时，对浮点数进行加法和减法是一个问题。另一个是使用atan2f仅返回 a 的函数，float然后在cos()and中使用该结果sin()。如果可能，实际上最好避免使用三角函数。
最后，该程序目前仅适用于两个对象。使用这种方案添加三分之一会很痛苦，因此更好的设计是std::vector<CelestialObject>通过首先考虑所有其他对象的位置和质量来计算每个对象上的净力来处理 a。我会把它留给你，但这至少应该让你朝着正确的方向开始。