Wie kann ich die Euler-Methode durch Runge-Kutta 4. Ordnung ersetzen, um die Bewegung des freien Falls in nicht konstanter Gravitationsgröße (z. B. freier Fall aus 10 000 km Höhe) zu bestimmen?
Bisher habe ich eine einfache Integration nach der Euler-Methode geschrieben:
while()
{
v += getMagnitude(x) * dt;
x += v * dt;
time += dt;
}
x Variable bedeutet aktuelle Position, v bedeutet Geschwindigkeit, getMagnitude (x) gibt die Beschleunigung auf x Position zurück.
Ich habe versucht, RK4 zu implementieren:
while()
{
v += rk4(x, dt) * dt; // rk4() instead of getMagintude()
x += v * dt;
time += dt;
}
Dabei ist der Funktionskörper von rk4 ():
inline double rk4(double tx, double tdt)
{
double k1 = getMagnitude(tx);
double k2 = getMagnitude(tx + 0.5 * tdt * k1);
double k3 = getMagnitude(tx + 0.5 * tdt * k2);
double k4 = getMagnitude(tx + tdt * k3);
return (k1 + 2*k2 + 2*k3 + k4)/6.0;
}
Aber etwas stimmt nicht, weil ich nur einmal mit RK4 (Beschleunigung) integriere. Die Integration der Geschwindigkeit mit RK4 ist nicht sinnvoll, da sie mit v * dt identisch ist.
Können Sie mir sagen, wie man Differentialgleichungen zweiter Ordnung mithilfe der Runge-Kutta-Integration löst? Sollte ich RK4 implementieren, indem ich die Koeffizienten k1, l1, k2, l2 ... l4 berechne? Wie kann ich das machen?
quelle
Antworten:
Es scheint ziemlich verwirrend zu sein, wie mehrstufige (z. B. Runge-Kutta) Methoden auf ODEs 2. oder höherer Ordnung oder ODE-Systeme angewendet werden sollen. Der Prozess ist sehr einfach, wenn Sie ihn erst einmal verstanden haben, aber ohne eine gute Erklärung vielleicht nicht offensichtlich. Die folgende Methode finde ich am einfachsten.
k1
k4
X
RHS( t, X )
Leider unterstützt C ++ solche Vektoroperationen nicht nativ, sodass Sie entweder eine Vektorbibliothek verwenden, Schleifen verwenden oder die einzelnen Teile manuell ausschreiben müssen.In C ++ können Siestd::valarray
den gleichen Effekt erzielen. Hier ist ein einfaches Arbeitsbeispiel mit konstanter Beschleunigung.quelle
typedef std::valarray<double> Vector
für häufig verwendete Typen. 3) Verwenden Sieconst int NDIM = 2
anstelle von#define
für Typensicherheit und Korrektheit. 4) Seit C ++ 11 können Sie den Körper von RHS einfach durch ersetzenreturn {X[1], 1}
. 5) In C ++ ist es wirklich ungewöhnlich (im Gegensatz zu C), zuerst Variablen zu deklarieren und später zu initialisieren. Variablen lieber an derselben Stelle zu deklarieren, an der Sie sie initialisieren (double t = 0.
usw.)RHS()
berechnet die rechte Seite der Differentialgleichung. Der Zustandsvektor X ist (x, v), also ist dX / dt = (dx / dt, dv / dt) = (v, a). Für Ihr Problem (wenn a = G * M / x ^ 2) sollte RHS zurückkehren{ X[1], G*M/(X[0]*X[0]) }
.