Hej!
Prosty przykład, przyjmijmy, że masz układ 2 rzędu, który sprowadzisz sobie (lub masz sprowadzony) do 2 równań różniczkowych liniowych zwyczajnych rzędu 1:
x1'(t) = 2x1(t) + x2(t)
x2'(t) = 3x2(t)
W metodzie Rungego Kutty 4 rzędu ze stałym krokiem całkowania (zakładam, że to o nią chodzi, gdyż jest chyba najpopularniejsza) rozwiązujemy równanie różniczkowe obliczając w każdym kroku wartość pochodnych (lewa strona równania) bazując na wartości funkcji (prawa strona) w poprzednim kroku.
Dla ułatwienia napiszmy metodę rhs (ang.right hand side), która będzie przyjmować jako argumenty stan systemu oraz przeładujmy operatory + i *.
struct state
{
static const int SYSTEM_DIMENSION = 2;
double x[SYSTEM_DIMENSION];
state& operator+(const double number)
{
for(int i = 0; i < SYSTEM_DIMENSION; i++)
x[i] += number;
return *this;
}
state& operator*(const double number)
{
for(int i = 0; i < SYSTEM_DIMENSION; i++)
x[i] *= number;
return *this;
}
state& operator+(const state& system)
{
for(int i = 0; i < SYSTEM_DIMENSION; i++)
x[i] += system.x[i];
return *this;
}
};
void printSystemState(const state& systemState)
{
for(int i = 0; i < systemState.SYSTEM_DIMENSION; i++)
std::cout << systemState.x[i] << std::endl;
}
state rhs(const state& systemState)
{
state newState;
newState.x[0] = 2*systemState.x[0] -1*systemState.x[1];
newState.x[1] = 3*systemState.x[1];
return newState;
}
Teraz możemy zaimplementować bardzo prostą wersję metody Rungego-Kutty (opieram się na wikipedii: https://pl.wikipedia.org/wiki/Algorytm_Rungego-Kutty).
W algorytmie ważne jest, aby zachować stały krok metody, dlatego warto poświęcić chwilę na opracowanie lepszego mechanizmu dzielenia czasu na przedziały.
state& rungeKutta(state& systemState, double finalTime)
{
float h = 0.1; //simulation step, s
int numberOfIntervals = finalTime/h;
for(int i = 0; i < numberOfIntervals; i++)
{
state k1, k2, k3, k4;
k1 = rhs(systemState) * h;
k2 = rhs(systemState + k1*h*0.5)*h;
k3 = rhs(systemState + k2*h*0.5)*h;
k4 = rhs(systemState) * h;
state deltax = (k1 + k2 * 2 + k3 * 2 + k4) * (1/6.0);
systemState = systemState + deltax;
}
return systemState;
}
int main() {
state system;
//initial conditions
system.x[0] = 1.5;
system.x[1] = 0.1;
std::cout << "Initial system state" << std::endl;
printSystemState(system);
rungeKutta(system, 3.0);
std::cout << "System state after 3 seconds" << std::endl;
printSystemState(system);
return 0;
}
U mnie daje to wynik:
Initial system state
1.5
0.1
System state after 3 seconds
35.0392
16.4262
Przedstawione rozwiązanie wymaga jeszcze trochę nakładu pracy w celu poprawy jakości działania i uproszczenia kodu, ale mam nadzieję, że udało mi się zasygnalizować możliwy sposób poradzenia sobie z problemem.
Pozdrawiam