Metoda Runge Kutty

Cześć, mam do napisania jako projekt z informatyki program, który implementuje algorytm Runge Kutty. Ma on wczytywać współczynniki równania różniczkowego zwyczajnego n-tego rzędu z pliku z tym, że nie miałem jeszcze w ogóle do czynienia z rónaniami różniczkowymi (na studiach mieliśmy tylko całkowanie funkcji jednej zmiennej) i mam problem z rozpocząciem samego programu. Jak program ten powinien wczytywać te współczynniki (w jakiej formie itp. ?). Najlepiej jeśli któś mógłby mi to pokazać na przykładzie jakiegoś równania.

Z góry dzięki za pomoc!

Współczynniki to zwykłe liczby zmiennoprzecinkowe, czyli double.

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