-1

我一直在尝试构建一个 Runge Kutta 四阶积分器来模拟简单的弹丸运动。我的代码如下

double rc4(double initState, double (*eqn)(double,double),double now,double dt)
{
        double k1 = eqn(initState,now);
        double k2 = eqn(initState + k1*dt/2.0,now + dt/2.0);
        double k3 = eqn(initState + k2*dt/2.0,now + dt/2.0);
        double k4 = eqn(initState + k3*dt, now + dt);

        return initState + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);
}

这在 while 循环中调用

while (time <= duration && yPos >=0)
                {

                        xPos = updatePosX(xPos,vx,timeStep);
                        yPos = updatePosY(yPos,vy,timeStep);


                        vx = rc4(vx,updateVelX,time,timeStep);
                        vy = rc4(vy,updateVelY,time,timeStep);

                        cout << "x Pos: " << xPos <<"\t y Pos: " << yPos << endl;

                        time+=timeStep;

                        myFile << xPos << "  " << yPos << "  " << vx << "  " << vy << endl;

                }

然而，与应该发生的事情相反，我的结果只是爆炸了。这里发生了什么？

4

1 回答 1

0

您的 rk4 代码看起来正确。但仅适用于标量微分方程。

您最肯定拥有的是一个维度大于 1 的耦合微分方程系统。在这里，您必须以向量形式应用积分方法。也就是说，x,y,vx,vy组合成一个 4 维（相位）状态向量，系统函数是向量值，k1,...k4是向量等。

作为高级说明，time <= duration对time+=timeStep;. 更好地使用time <= duration-timeStep/2在time循环结束时接近duration.

阅读关闭的上一个问题的代码，我发现您对微分方程的概念有疑问。您不应该在 RK4 实现中使用欧拉步骤的结果作为加速。没有空气摩擦的弹道运动系统是

dotx = vx
doty = vy
dotvx = 0
dotvy = -g

您必须以矢量形式实现，例如

eqn(t, [x,y,vx,vy]) // where X = array of double of dimension 4
   { return [vx,vy,0,-g]; }
于 2015-01-29T14:45:57.053 回答