对 Runge Kutta 第 4 种方法的方法？

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 returns 每个 x、y 和 z 的加速度。

Answer 1

您还没有定义您的方法，但让我感到震惊的是您将结果与输入混合在一起。由于 Runge-Kutta 是一种计算 y_(n+1) = y_n + hsum(b_ik_i) 的方法，我希望您的解决方案是让您的 _n 条款在右边，而您的 (n+1) 条款在左边。这不是你在做什么。相反，s(n+1) 依赖于 r_(n+1) 而不是 r_n，t_(n+1) 依赖于 s_(n+1)，等等。这有点像您试图限制正在使用的变量数量的错误。

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

Answer 2

你有一些加速功能

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

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

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.

对 Runge Kutta 第 4 种方法的方法？

Right Runge Kutta 4th method approach?

c++

vector

runge-kutta