Правильный Runge Kutta 4-й метод подхода?
У меня есть этот код Руджи Кутты. Однако один упомянул мой подход неверен. И я не мог понять, почему от него, так что кто-нибудь здесь, кто мог бы дать подсказку, почему этот путь не так?
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);
==== ==== EDIT
Vector3d хранит координаты, x, y, z.
GetAcceleration возвращает ускорение для каждого x, y и z.
Решение
У вас есть функция ускорения
a(p,q) where p=(x,y,z) and q=(vx,vy,vz)
Ваша система заказа 1, которую можно решить с помощью RK4,
dotp = q
dotq = a(p,q)
Этапы метода РК включают смещение вектора (ов) состояния
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,
Другие решения
Вы не определили свои методы, но меня бросает в глаза то, что вы смешиваете свои результаты со своими данными. Так как Рунге-Кутта является методом для расчета y_ (n + 1) = y_n + hсумма (b_ik_i), я ожидаю, что ваше решение будет держать ваши _n условия правильно, а ваши (n + 1) члены слева. Это НЕ то, что вы делаете. Вместо этого(n + 1) зависит от r_ (n + 1), а не от r_n, t_ (n + 1) от s_ (n + 1) и так далее. Это пахнет ошибкой, когда вы пытались ограничить количество используемых переменных.
Имея это в виду, можете ли вы указать фактические промежуточные значения вычислений, которые генерирует ваша программа, и сравнить их с предполагаемыми промежуточными значениями?