Runge Kutta Verfahren 4. Ordnung

  • C#

Es gibt 36 Antworten in diesem Thema. Der letzte Beitrag () ist von ~blaze~.

    Runge Kutta Verfahren 4. Ordnung

    Hallo,

    ich versuche gerade das klassische Runge Kutta Verfahren nachzuvollziehen und will es auf Bewegung von Satelliten anwenden.
    Folgende Seite habe ich gefunden: dieter-heidorn.de/Physik/SS/K0…3_Gravitationsgesetz.html

    Erklärt gut das Gravitationsgesetz und zeigt seine Applikation.
    Zur Simulation wird das explizite eulersche Verfahren genutzt.

    Das eben dadurch die Akkuratesse der Simulation zu wünschen übrig lässt, habe ich mir überlegt es mit RK4 zu realisieren.
    Folgende Frage habe ich:

    Auf der ehemals vorhandenen Seite letsthinkabout.us fand ich folgende Implementierung des RK4 in C#:
    web.archive.org/web/2016033113….us/post/runge-kutta-in-c

    Müsste ich jetzt quasi 2 Equations definieren, einmal für Velocity, Position und Acceleration?
    Das ist etwas, was ich nicht richtig nachvollzogen haben.

    Was wäre denn der Initialwert?
    Etwa die Beschleunigung zum Zeitpunkt 0?

    Lieben Dank für eure Hilfen.

    ~blaze~: Thema verschoben
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „~blaze~“ ()

    @φConst Jou.
    Wir haben einst im Strudium das allgemeine Pendel (große Ausschläge) mit RK programmiert. Die RK-Bibliothek war gegeben, den Aufruf und die Darstellung war der Plan.
    Die Gleichungen gehorchen dieser Form: f(i) = h * {f[y(j), t]} mit h - Schrittweite, t - Koordinate (Zeit).
    Mit den DGLn
    • y(0)' = y(1) #Geschwindigkeit#
    • y(1)' = -sin(y(0)) #Beschleunigung#
    und den Startbedingungen (Einheiten weiß ich jetzt nicht)
    • y(0) = 0 #Geschwindigkeit#
    • y(1) = 1 #Beschleunigung#
    kamen da die richtigen Kurven.
    Kannst dies ja mal ausprobieren, um das Herangehen an Dein Problem zu tunen.
    Falls Du diesen Code kopierst, achte auf die C&P-Bremse.
    Jede einzelne Zeile Deines Programms, die Du nicht explizit getestet hast, ist falsch :!:
    Ein guter .NET-Snippetkonverter (der ist verfügbar).
    Programmierfragen über PN / Konversation werden ignoriert!

    Hallo,
    also soweit die Implementierung mit Euler:

    C#-Quellcode

    1. public class CelestialBody
    2. {
    4. public static readonly double ASTRONOMICUNIT = 149597870700;
    5. public static readonly double GRAVITATIONCONSTANT = 6.67408e-11;
    6. public static readonly double SUN_G = 1.989E30 * GRAVITATIONCONSTANT;
    7. public static readonly double Exponent = 3;
    9. public double X { get; set; }
    10. public double Y { get; set; }
    12. public double Ax { get; private set; }
    13. public double Ay { get; private set; }
    15. public double Vx { get; set; }
    16. public double Vy { get; set; }
    18. public double Time { get; private set; }
    19. public double Delta { get; set; }
    21. public CelestialBody(PointF initPos, PointF initVel, double d)
    22. {
    23. X = initPos.X;
    24. Y = initPos.Y;
    26. Vx = initVel.X;
    27. Vy = initVel.Y;
    29. Delta = d;
    31. }
    33. public void Update()
    34. {
    35. double distance = Math.Sqrt(Math.Pow(X, 2) + Math.Pow(Y, 2));
    37. Ax = -SUN_G * (X / Math.Pow(distance, Exponent));
    38. Ay = -SUN_G * (Y / Math.Pow(distance, Exponent));
    40. Vx += Ax * Delta;
    41. Vy += Ay * Delta;
    43. X += Vx * Delta;
    44. Y += Vy * Delta;
    45. }
    47. public void Render(Graphics g)
    48. {
    49. g.FillEllipse(Brushes.Blue, new RectangleF((float)(X / ASTRONOMICUNIT) * 200 + 350f, (float)(Y / ASTRONOMICUNIT) * 200 + 350f, 15, 15));
    50. }
    52. }


    Gibt man nun die Initialwerte 149597870700 Meter und 29780 m/s simuliert man die Bewegung der Erde um die Sonne, klappt wunderbar mit dem Delta = 60 * 60 * 24 * 7 * 0.1

    Wären nun, wenn ich RK4 implementiere, lediglich Acceleration betroffen?
    Weil im Grunde ist ja Geschwindigkeit und Position eigentlich abhängig von Beschleunigung.
    Das Ziel ist also y(0) = V0 , was wäre denn y'(x)=f(x,y) ? Etwa Velocity t+delta = Velocity t + delta * a?
    Und Gott alleine weiß alles am allerbesten und besser.

    Prinzipiell musst du immer, egal welches Lösungsverfahren du nutzt den gleichen Ablauf durchgehen:
    1) Problem spezifieren
    2) DGLs/PDEs aufstellen
    3) Lösungsverfahren auswählen
    4) System lösen

    1) hast du bereits durch die von dir verlinkte Seite gelöst.
    Nun hast du (wie grundsätzlich erstmal bei allen rein mechanischen Problemen) das Gleichungssystem . Das möchtest Du jetzt mit einem expliziten Runge-Kutta-Verfahren lösen. Gut, okay. Dabei ist das vorgehen (zumindest für die explizten Lösungsverfahren und autonome Systeme) eigentlich immer identisch: du berechnest für jede der beiden Komponenten des Gleichungsystems den nächsten Schritt. Das, was auf der Seite x genannt wird, kannst du bei dir als die Zeit betrachten. Da das System autonom ist, ist für dich f(x,y) = f(y).

    Falls das deine Fragen noch nicht ganz geklärt hat, frag nochmal nach.

    Lieben Dank!
    Ich meine jetzt mein Problem erfasst haben zu können:

    Das Problem ist f(t)!
    Dem Anschein nach ist ja a(x) nicht von der Zeit abhängig, ebenfalls Geschwindigkeit nicht, sondern von Delta.
    Wie müsste ich die Gleichung umstellen um ein Wert zum Zeitpunkt t zu erhalten?

    Ich hätte jetzt an sowas gedacht:




    Ist das denn so?
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 2 mal editiert, zuletzt von „φConst“ ()

    So, hab das jetzt besser illustriert wie ich das meine und komm zur Frage: Wird das kein Stackoverflow?
    Hier:


    Anderenfalls hab ich ja Initialwerte, am besten ich probiere das mal _ irgendwann aus.
    Falls du was findest, sag jetzt schon ^^


    Addendum: Kleiner Fehler in den Summenzeichen: Es muss immer im Summenzeichen der Initialwert subtrahiert werden, sonst hat man am Ende Initialwert * (t + 1).
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 2 mal editiert, zuletzt von „φConst“ ()

    @φConst Setz mal auf Deinen RK-Algo mein Pendel an, probiere das einfach mal aus.
    An so ner Skizze lässt sich kein Nicht-Stackoverflow absichern.
    Und:
    Klar brauchst Du auch die Orte, hab ich oben vergessen. Sorry.
    Also:
    Ort, Geschwindigkeit, Beschleunigung und los.
    Falls Du diesen Code kopierst, achte auf die C&P-Bremse.
    Jede einzelne Zeile Deines Programms, die Du nicht explizit getestet hast, ist falsch :!:
    Ein guter .NET-Snippetkonverter (der ist verfügbar).
    Programmierfragen über PN / Konversation werden ignoriert!

    Hi
    wenn
    $x_0$
    und
    $v_0$
    gegeben sind und du alle Konstanten hast, um die Beschleunigung a(t, x) zu berechnen, kannst du ja eigentlich alles bestimmen, was du brauchst.

    Deine Beschleunigung ist abhängig vom Ort und dieser ist wiederum abhängig von der Beschleunigung:
    x'' = a(t, x)

    Deine Beschleunigung erhältst du ja aus den Gleichungen auf der von dir geposteten Webseite.
    Du musst jetzt zwei mal integrieren - einmal für Geschwindigkeit und einmal für Ort:

    Deine Funktion a(t, x) ist gegeben, das ist deine Funktion für die Beschleunigung.
    Jetzt berechnest du für
    $v_{n+1} = v(t + h, x)$
    den neuen Wert aus dem alten
    $v_n = v(t, x)$
    . Hierbei berechnest du k1, ..., k4 nach dem gegebenen Prinzip. In deinem Link zum Runge-Kutta-Verfahren ist jeweils xn deine alte Geschwindigkeit, h die Schrittweite, yn dein alter Ort:
    $ k_1 = h_v \cdot a(t, x_n)\\ k_2 = h_v \cdot a(t + \frac{1}{2}h_v, x_n + \frac{1}{2}k_1)\\ k_3 = h_v \cdot a(t + \frac{1}{2}h_v, x_n + \frac{1}{2}k_2)\\ k_4 = h_v \cdot a(t + h_v, x_n + k_3)\\ v_{n+1} = v_n + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) $

    Der ki-te Werte basiert immer auf dem k(i-1)-ten Wert. Daher berechnest du die einfach der Reihenfolge nach.

    Für den Ort ist die zu integrierende Funktion dann eben v(t, x). Wenn du an dieser Stelle erneut den Runge Kutta anwendest, erhältst du
    $ k_1 = h_x \cdot v(t, x_n)\\ k_2 = h_x \cdot v(t + \frac{1}{2}h_x, x_n + \frac{1}{2}k_1)\\ k_3 = h_x \cdot v(t + \frac{1}{2}h_x, x_n + \frac{1}{2}k_2)\\ k_4 = h_x \cdot v(t + h_x, x_n + k_3)\\ x_{n+1} = x_n + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) $


    Was jetzt der Unterschied ist, ist, dass v(...) jeweils ebenfalls über den Runge-Kutta definiert ist (weshalb die Indices von v btw. auch nicht ganz so toll gewählt sind. Du hast quasi auch ein
    $v_{n + \frac{1}{2}}$
    ), d.h. du musst jeweils die Schrittweite
    $h_v$
    für die Berechnung von v anpassen.

    Und in deinem Link wird der Runge-Kutta halt über einen Delegaten ausgeführt, was das ganze halbwegs generisch macht. Es sollte auf jeden Fall auch für Vektoren funktionieren.

    Edit: Hab's mal in LaTeX neu geschrieben, damit es übersichtlicher ist...

    Viele Grüße
    ~blaze~

    Dieser Beitrag wurde bereits 5 mal editiert, zuletzt von „~blaze~“ ()

    Lieben Dank euch allen,
    stimmt, an Integrieren hab ich gar nicht gedacht.

    Im Grunde habe ich folgenden Funktionsterm

    ax(t) = - const. * M * x(t)/sqrt(x(t)^2+y(t)^2)
    ay(t) = - const. * M * y(t)/sqrt(x(t)^2+y(t)^2)

    Was ich aber bei dir nicht ganz nachvollzogen habe ist, wieso du a(t,x) schreibst, wo doch die Position x (oder y) selber abhänging von t ist.
    Das irritiert mich ja die ganze Zeit: Es ergibt sich eine Art Ring-Formel die am Ende auf sich selbst verweist.
    Und die Formel aus der Website gibt ja nicht explizit an, dass a von x abhängig ist, sondern eher x(t), oder y(t). Deswegen meine Frage: Wie kann ich die Formel umschreiben, sodass ich zu jedem Zeitpunkt t, einen Wert für a kriege.

    Falls ich gerade Mist rede, verzeiht mir, will mal wieder einsteigen ^^

    @Rod werde ich mal anwenden_
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „φConst“ ()

    Wie gesagt, dein x ist das x aus dem letzten Schritt. Schau dir nochmal meine Formeln an, da steht drin, was zu tun ist.
    Dein a(t, x) ist durch die Parameter t und x festgelegt, du setzt nur noch x0, xn, k1, usw. ein.
    Runge-Kutta ist eine Näherung des Integrals. Du hast ja im Prinzip eine Differentialgleichung.
    Dein r ist von x abhängig.

    Viele Grüße
    ~blaze~

    Ich meine jetzt verstanden zu haben:

    k1 für v wäre dann : hv * a(0, 149597870700)
    und vn ist 29780 m/s

    der Wert wird gespeichert, dann wird inkrementiert, x kriegt nun v, neues x berechnen und v erhält dann den Wert von neuem x, alles gut.
    Mich irritierte nur halt a(t,x) denn



    gibt keinen Wert zum Zeitpunkt t zurück, sondern zu x, oder y.

    Ich glaube es ist einfach nur die Formalie, ich meine aber nachvollzogen haben zu können.
    Werde es _ mal nachher implementieren.

    Schade das es kein hilfreich Button gibt, danke euch.
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „φConst“ ()

    φConst schrieb:

    gibt keinen Wert zum Zeitpunkt t zurück, sondern zu x, oder y.
    Da wirst Du Schwierigkeiten bekommen, denn das Verfahren geht davon aus, aus Startwerten und den Differentialgleichungen für Geschwindigkeit und Beschleunigung die Position als Funktion der Zeit zu berechnen.
    Die Geschwindigkeit ist ja die erste, die Beschleunigung die zweite Ableitung des Ortes nach der Zeit.
    Bei Dir käme dann wohl die Bahnkurve als solche raus, zumal Dein Darstellungscode eine ellipsoide Form impliziert.
    Bei meinem Projekt damals kam so was raus wie ein Sinus, und bei großen Auslenkungen eine sehr rechteckige abgerundete Kurve.
    ====
    Kann es sein, dass in Deinem Bildchen iwas falsch ist?
    x0 kommt zwei Mal vor, y0 und y je einmal.
    Falls Du diesen Code kopierst, achte auf die C&P-Bremse.
    Jede einzelne Zeile Deines Programms, die Du nicht explizit getestet hast, ist falsch :!:
    Ein guter .NET-Snippetkonverter (der ist verfügbar).
    Programmierfragen über PN / Konversation werden ignoriert!

    Ja, eben das frage ich mich ja auch.
    Deswegen irritiert es mich ja die ganze Zeit wieso ihr a(t,x) schreibt, als sei a von t abhängig ( was es ja ist, aber die Seite es nicht so formuliert).

    Hier mein Versuch das zu implementieren:


    C#-Quellcode

    1. public class CelestialBody
    2. {
    4. public static readonly double ASTRONOMICUNIT = 149597870700;
    5. public static readonly double GRAVITATIONCONSTANT = 6.67408e-11;
    6. public static readonly double SUN_G = 1.989E30 * GRAVITATIONCONSTANT;
    7. public static readonly double Exponent = 3;
    9. private double tmpVx, tmpVy, tmpX, tmpY;
    11. private double deltaVx, deltaVy, deltaXx, deltaYy;
    12. private double k1, k2, k3, k4;
    14. public CelestialBody(double x0, double y0, double vx0, double vy0)
    15. {
    16. tmpVx = vx0;
    17. tmpVy = vy0;
    19. tmpX = x0;
    20. tmpY = y0;
    22. deltaVx = deltaVy= deltaXx = deltaYy = .1f;
    23. }
    25. public void Update()
    26. {
    27. // VELOCITY
    29. k1 = deltaVx * ax(tmpX);
    30. k2 = deltaVx * ax(tmpX + (1.0 / 2.0) * k1);
    31. k3 = deltaVx * ax(tmpX + (1.0 / 2.0) * k2);
    32. k4 = deltaVx * ax(tmpX + k3);
    33. tmpVx = tmpVx + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    35. k1 = deltaVy * ay(tmpY);
    36. k2 = deltaVy * ay(tmpY + 1.0 / 2.0 * k1);
    37. k3 = deltaVy * ay(tmpY + 1.0 / 2.0 * k2);
    38. k4 = deltaVy * ay(tmpY + k3);
    39. tmpVy = tmpVy + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    41. // POSITION
    43. k1 = deltaXx * vx(tmpVx);
    44. k2 = deltaXx * vx(tmpVx + (1.0 / 2.0) * k1);
    45. k3 = deltaXx * vx(tmpVx + (1.0 / 2.0) * k2);
    46. k4 = deltaXx * vx(tmpVx + k3);
    47. tmpX = tmpX + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    49. k1 = deltaYy * vy(tmpVy);
    50. k2 = deltaYy * vy(tmpVy + (1.0 / 2.0) * k1);
    51. k3 = deltaYy * vy(tmpVy + (1.0 / 2.0) * k2);
    52. k4 = deltaYy * vy(tmpVy + k3);
    53. tmpY = tmpY + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    55. Console.WriteLine(tmpX + " ; " + tmpY);
    57. }
    61. private double ax(double x)
    62. {
    63. double distance = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(tmpY, 2));
    64. return -SUN_G * (x / Math.Pow(distance, Exponent));
    65. }
    67. private double ay(double y)
    68. {
    69. double distance = Math.Sqrt(Math.Pow(tmpX, 2) + Math.Pow(y, 2));
    70. return -SUN_G * (y / Math.Pow(distance, Exponent));
    71. }
    73. private double vx(double x)
    74. {
    75. double distance = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(tmpY, 2));
    76. return -(1.58500795e27 * Math.Pow(x, 2)) / (23880051 * Math.Pow(distance, 3));
    77. }
    79. private double vy(double y)
    80. {
    81. double distance = Math.Sqrt(Math.Pow(tmpX, 2) + Math.Pow(y, 2));
    82. return -(1.58500795e27 * Math.Pow(y, 2)) / (23880051 * Math.Pow(distance, 3));
    83. }


    Initialisiert wird mit:

    C#-Quellcode

    1. CelestialBody test = new CelestialBody(149597870700, 0, 0, 29780);


    Ergebnisse sind arg falsch:
    Und Gott alleine weiß alles am allerbesten und besser.

    @φConst Was sagt denn die Beschreibung dazu, was ax, ay und vx, vy physikalisch genau verkörpert, wenn nicht Geschwindigkeit und Beschleunigung?
    Du kannst alle Deltas zu einer einzigen Variable zusammenfassen, gekoppelte DGLn können nur mit einer gemeinsamen Schrittweite gelöst werden.
    Welche physikalische Einheit hat Deine Schrittweite?
    Falls Du diesen Code kopierst, achte auf die C&P-Bremse.
    Jede einzelne Zeile Deines Programms, die Du nicht explizit getestet hast, ist falsch :!:
    Ein guter .NET-Snippetkonverter (der ist verfügbar).
    Programmierfragen über PN / Konversation werden ignoriert!

    UPDATE: Also es funktioniert jetzt!

    C#-Quellcode

    1. public class CelestialBody
    2. {
    4. public static readonly double ASTRONOMICUNIT = 149597870700;
    5. public static readonly double GRAVITATIONCONSTANT = 6.67408e-11;
    6. public static readonly double SUN_G = 1.989E30 * GRAVITATIONCONSTANT;
    7. public static readonly double Exponent = 3;
    9. private double X, Y, Vx, Vy;
    10. private double tmpVx, tmpVy, tmpX, tmpY;
    12. private double delta;
    13. private double k1, k2, k3, k4;
    15. public CelestialBody(double x0, double y0, double vx0, double vy0)
    16. {
    17. tmpVx = vx0;
    18. tmpVy = vy0;
    20. tmpX = x0;
    21. tmpY = y0;
    23. delta = 60 * 60 * 24 * 7 * 0.1;
    24. }
    26. public void Update()
    27. {
    28. // VELOCITY
    30. double halfdx = .5d * delta;
    32. k1 = delta * ax(tmpX);
    33. k2 = delta * ax(tmpX + halfdx * k1);
    34. k3 = delta * ax(tmpX + halfdx * k2);
    35. k4 = delta * ax(tmpX + k3);
    36. Vx = tmpVx + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    38. tmpVx = Vx;
    40. k1 = delta * ay(tmpY);
    41. k2 = delta * ay(tmpY + halfdx * k1);
    42. k3 = delta * ay(tmpY + halfdx * k2);
    43. k4 = delta * ay(tmpY + k3);
    44. Vy = tmpVy + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    46. tmpVy = Vy;
    48. // POSITION
    50. k1 = delta * vx(time, tmpX);
    51. k2 = delta * vx(time + halfdx, tmpX + halfdx * k1);
    52. k3 = delta * vx(time + halfdx, tmpX + halfdx * k2);
    53. k4 = delta * vx(time + delta , tmpX + k3);
    54. X = tmpX + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    56. tmpX = X;
    58. k1 = delta * vy(time, tmpY);
    59. k2 = delta * vy(time + halfdx, tmpY + halfdx * k1);
    60. k3 = delta * vy(time + halfdx, tmpY + halfdx * k2);
    61. k4 = delta * vy(time + delta, tmpY + k3);
    62. Y = tmpY + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    64. tmpY = Y;
    66. Console.WriteLine(tmpX + " ; " + tmpY);
    68. time += delta;
    69. }
    71. double time;
    74. private double ax(double x)
    75. {
    76. double distance = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(tmpY, 2));
    77. return -SUN_G * (x / Math.Pow(distance, Exponent));
    78. }
    80. private double ay(double y)
    81. {
    82. double distance = Math.Sqrt(Math.Pow(tmpX, 2) + Math.Pow(y, 2));
    83. return -SUN_G * (y / Math.Pow(distance, Exponent));
    84. }
    86. private double vx(double t, double x)
    87. {
    88. return tmpVx + ax(x) * t; ;
    89. }
    91. private double vy(double t, double y)
    92. {
    93. return tmpVy + ay(y) * t;
    94. }
    96. public void Render(Graphics g)
    97. {
    98. g.FillEllipse(Brushes.Blue, new RectangleF((float)(tmpX / ASTRONOMICUNIT) * 200 + 350f, (float)(tmpY / ASTRONOMICUNIT) * 200 + 350f, 15, 15));
    99. }
    101. }


    Hab einfach für v(t,x) = a(x)*t+v(t-1) genommen...
    Aber wir müssten alle tot sein, die Erde fliegt in die Sonne, why?

    Also hab jetzt time = time + delta rausgenommen und stattdessen eben
    f(0,x)
    f(halfdx,x)
    f(halfdx,x)
    f(dx,x) berechnet:


    Woran liegt das?
    Überlauffehler?
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „φConst“ ()

    Mit v(t, x) = a(x) * t + v(t-1) hast du jetzt doch wieder einen Euler-Schritt gemacht, oder?

    Ersetze übrigens Pow(x, 2) durch x * x. Das ist weit effizienter, als der ineffiziente Power-Algorithmus.

    Die Rechnung in ax und ay macht nur Sinn, wenn du Vektoren verwendest. Formuliere eine Funktion a(t, x), wobei x ein Vektor ist. Die Formel lautet
    $-G \cdot \frac{r}{{\lvert r \rvert}^3}$

    du hast quasi
    $-G \cdot \frac{r}{{\lvert r \rvert}} \cdot \frac{1}{{\lvert r \rvert}^2}$


    hierbei ist
    $\frac{r}{{\lvert r \rvert}}$
    dann eben der normalisierte relative Vektor r zwischen den beiden Objekten, zwischen denen die Beschleunigung vorherrscht und 1/r^2 die Länge zum Quadrat. Aus Ersterem folgt die Richtung der Beschleunigung.

    Bei mehreren Körpern würde sich die Beschleunigung nach dem Superpositionsprinzip aufaddieren, d.h. einfach alle Beschleunigungen aufaddieren. Es ist dann allerdings die Frage, ob du für einen Schritt die k1, ..., k4 für alle Objekte benötigst und im Schritt t+h bzw. t+h/2 abfrägst, da bin ich gerade unschlüssig.

    Btw.: Meine Gleichungen sollten schon passen, sofern kein Flüchtigkeitsfehler vorherrscht. Es ergibt sich die neue Zeit zu
    $t_{n+1} = t_n + h$
    mit
    $t_0 = 0$
    .

    Viele Grüße
    ~blaze~

    @~blaze~
    Hab schon geahnt das es wieder eulersches Verfahren ist, wusste eben nicht was wie der Funktionsterm für v(t) lautet... a(t,x) integriert kommen sehr komische Sachen raus.

    Formuliere eine Funktion a(t, x), wobei x ein Vektor ist.


    Wo in der Formel -G*r/|r|^3 hast du ein t? '
    Ich hab echt keine Ahnung wie du auf dieses t kommst, der ist doch nicht ausschlaggebend für die Funktion. Bin einfach nur bekloppt glaub ich :D

    Also salopp gesagt ist die Formel die gesucht es ja dann
    f(x)= -G * x/|x|^3 wobei eben x ein Vektor ist (zB: (149 000 000 000; 0)).
    Aber wieso ist jetzt diese Funktion nochmal abhängig von t?
    Wenn hier irgendetwas abhängig von t sein soll haben wir doch automatisch einen Eulerschritt irgendwo.
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 3 mal editiert, zuletzt von „φConst“ ()

    Das r ist entsprechend (x2 - x1), wobei x2 und x1 die Punkte sind, an denen sich die Objekte aufhalten.
    Du hast an dieser Stelle kein t. Die Beschleunigung hängt nur über den Ort von der Zeit ab und der ist ja offensichtlich von der Zeit abhängig.
    Das Runge-Kutta-Verfahren ist sehr ähnlich zum Euler-Verfahren. Es macht quasi mehrere verschiedene Eulerschritte und vereinigt diese dann zu einer Lösung. Die k2, ... k4 sind jeweils Näherungen auf Basis des jeweiligen Vorgängers und zum Schluss erfolgt die Gewichtung aller, wobei k1 mit 1 gewichtet wird, k2 mit 2, k3 mit 2 und k4 mit 1.

    Du besitzt zu jedem Durchlauf x und v des vorherigen Zeitschritts. Ausgehend von diesem berechnest du die Beschleunigung. Während du beim normalen Euler dann mit diesem Ergebnis weiterrechnest (du integrierst ja quasi anhand der Tangente, die durch die gegebene Funktion gelegt wird), verwendet das Runge-Kutta-Verfahren das neu ermittelte Ergebnis weiter, um eine präzisere Schätzung durchzuführen. D.h. deine k1, ..., k4 sind ebenfalls Vektoren, die eine Schätzung für die Geschwindigkeit angeben - eben auf Basis der vorliegenden Geschwindigkeit v bzw. des vorliegenden Orts x oder dem vorhergehenden ki.

    Vergiss übrigens nicht: Das x ist das im Parameter x übergebene x, nicht das aus der Variablen.

    Viele Grüße
    ~blaze~

    @~blaze~
    Ach, ... ja,... du meinst also :

    sind meine beiden "Integrale"?

    Dieses f(t, x) irritiert mich die ganze Zeit.

    Also mein Funktionsterm für v(t,x..) ist doch dann:



    ?
    Falls ja, meine ich das jetzt ganz verstanden zu haben.

    Addendum: Kann es sein das du bei

    am Ende "* t "vergessen hast?
    Also Vn+1 = Vn + 1/6(...)*t
    Und Gott alleine weiß alles am allerbesten und besser.

    Dieser Beitrag wurde bereits 2 mal editiert, zuletzt von „φConst“ ()

Facebook
Telefonspende