Exploding Runge Kutta Method

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 (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]; }