Метод взрыва Рунге Кутты
Я пытался построить интегратор 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;
}
Однако вопреки тому, что должно произойти, мои результаты просто взорвались. Что тут происходит?
Решение
Ваш код 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]; }
Другие решения