爆炸龙格库塔法

Question

我一直在尝试构建一个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;
                }

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

Answer 1

您的 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]; }