右龙格-库塔第四种方法

Question

我有这个runge-kutta代码。然而，有人提到我的做法是错误的。我真的不明白为什么他，所以这里的任何人，谁能暗示为什么这种方式是错误的？

    Vector3d r = P.GetAcceleration();
    Vector3d s = P.GetAcceleration() + 0.5*m_dDeltaT*r;
    Vector3d t = P.GetAcceleration() + 0.5*m_dDeltaT*s;
    Vector3d u = P.GetAcceleration() + m_dDeltaT*t;
    P.Velocity += m_dDeltaT * (r + 2.0 * (s + t) + u) / 6.0);

===编辑===

Vector3d存储坐标x，y，z。

GetAcceleration返回每个x、y和z的加速度。

Answer 1

您有一些加速功能

a(p,q) where p=(x,y,z) and q=(vx,vy,vz)

可以通过RK4解决的订单1系统是

dotp = q
dotq = a(p,q)

RK方法的各阶段涉及状态向量的偏移

k1p = q
k1q = a(p,q)
p1 = p + 0.5*dt*k1p
q1 = q + 0.5*dt*k1q
k2p = q1
k2q = a(p1,q1)
p2 = p + 0.5*dt*k2p
q2 = p + 0.5*dt*k2q
k3p = q2
k3q = a(p2,q2)

等等。可以为每个步骤调整点P的状态向量，保存原始坐标，也可以使用P的临时副本来计算k2, k3, k4。

Answer 2

您还没有定义方法，但让我感到惊讶的是，您将结果与输入混合在一起。由于Runge-Kutta是一种计算y_（n+1）=y_n+hsum（b_ik_i）的方法，我希望您的解决方案将_n项保持在右侧，将（n+1）项保持在左侧。这不是你要做的。相反，s（n+1）依赖于r_。

考虑到这一点，您能否指出程序生成的计算的实际中间值，并将其与预期中间值进行比较？