0

首先，我正在使用非常接近欧拉方法的方法来计算我的星球的位置。我知道这不是最准确的方法，但我一直在使用velocity verlet 近一个星期，但无法让它发挥作用。我的问题是我无法让我的星球围绕太阳循环，x 或 y 位置不断增加。任何帮助表示赞赏！谢谢！这是我的代码：

void updatePosition(CelestialObject object1, CelestialObject object2 )

{ // -----------------------------X calculations-----------------------------------
    //calc force
    float forceX = forceFuncX(object1.getX(), object2.getX(),object1.getY(), object2.getY(), object1.getMass(), object2.getMass());`

    //accel calc
    float AX = accelerationFuncX(forceX,object1.getX(), object2.getX(),object1.getMass());
    agk::PrintC("Accel X: "); 
    float AXprint = AX*dt;
    agk::Print(AXprint);

    //velocity
    float VX  = object1.getVX();
    VX = VX + AX*dt;
    agk::PrintC("Velocity X: ");
    agk::Print(VX);

    //positionCalc
    float X = object1.getX();
    X = X + VX*dt;
    agk::PrintC("Position X: ");
    agk::Print(X);

    //-------------------------Y calculations------------------------------------
    //force
    float forceY = forceFuncY(object1.getX(), object2.getX(),object1.getY(), object2.getY(), object1.getMass(), object2.getMass());

    //accel
    float AY = accelerationFuncY(forceY,object1.getY(),object2.getY(), object1.getMass()); //y
    agk::PrintC("Accel Y: ");
    float AYprint = AY*dt;
    agk::Print(AYprint);

    //velocity
    float VY = object1.getVY();
    VY = VY + AY*dt;
    agk::PrintC("Velocity Y: ");
    agk::Print(VY);

    //position
    float Y = object1.getY();
    Y = Y + VY*dt;
    agk::PrintC("Position Y: ");
    agk::Print(Y);

    object1.setPosition(X, Y);
    agk::CreateParticles(X,Y);
}

这是它调用的函数：

double forceFuncX(float object1x,float object2x,float object1y, float object2y, double mass1, double mass2)
{
    float d = object1x - object2x;
    float r = sqrt(pow(object2x - object1x,2) + pow(object2y-object1y,2));;
    //float r = sqrt(pow(object1x-object2x,2)+pow(object1y-object2y,2));
    //double F = (G*(mass1*mass2))/pow(d,2);
    float F = (G*(mass1*mass2))/(r*r);

    return F;
}


double forceFuncY(float object1x,float object2x,float object1y, float object2y, double mass1, double mass2)
{
    float d = object1y - object2y;
    float r = sqrt(pow(object2x - object1x,2) + pow(object2y-object1y,2));;
    //float r = sqrt(pow(object1x-object2x,2)+pow(object1y-object2y,2));
    //double F = (G*(mass1*mass2))/pow(d,2);
    float F = (G*(mass1*mass2))/(r*r);

    return F;
}

和

float accelerationFuncX(float force, float object1x, float object2x, double mass) //gives the acceleration of an object
{
    float accel = (force*(object2x-object1x))/mass;
    return accel;
}

float accelerationFuncY(float force, float object1y, float object2y, double mass)
{
    float accel =(force*(object2y-object1y))/mass;
    return accel;
}
4

2 回答 2

1

由于您使用的是笛卡尔坐标系，您应该会发现沿 x 的力在某些配置中为正，而在其他配置中为负，并且沿 y 相同。特别是，力 F 是一个矢量，它从行星的质心指向它绕行的太阳的质心。所以代码有两个问题：

  1. 你返回 Fy 和 Fx，它们总是正的，并且等于两个物体之间的引力：F=G * m1 * m2/r^2。但 F 是力矢量的长度。力矢量的 Fx 和 Fy 分量将类似于 Fx = - F cos alpha 和 Fy = - F sin alpha，其中 alpha 是行星相对于 x 轴的角度。例如，如果 alpha=0，并且太阳位于坐标系的中心 (x=y=0)，则 Fx = -F（因为行星位于正 x 而太阳位于 x=0，因此力指向负X）。

  2. 加速度 a 是 F/m，而你有 vec F*(delta x) / m，其中 vec F 是力矢量，delta x 是行星和太阳之间的 x 距离。a 与力矢量的方向相同。所以如果你想计算 m1 质量的加速度（比如行星），你应该有类似 ax = Fx/m1 和 ay=Fy/m1 的东西。

根据上述内容修复您的代码应该很简单，因此我将由您自行完成，但如果不清楚，请发表评论。

于 2014-04-07T17:12:26.797 回答
0

为避免使用角度并使用三角函数将其向前和向后转换，请使用矢量力

F=(Fx, Fy)=-(x,y)*G * m1 * m2/r^3。

另外，使用教程“Moving Stars Around”来参考实现各种辛积分和非辛积分方法和Hairer-etal 的论文。对于它背后的理论和历史。

于 2014-04-07T21:58:35.193 回答