Runge Kutta Verfahren 4. Ordnung

  • C#

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

    Ja genau. Das sind relativ gute Näherungen. Wenn sie zu ungenau sind, musst du den Zeitschritt nochmal unterteilen. Gerade für sehr große Geschwindigkeiten sind die numerischen Näherungen nicht so gut geeignet.
    Das * t sollte nicht nötig sein (wenn dann wäre es eh *h), da es bereits bei den k1, ..., k4 aufmultipliziert wurde. Wann hv bzw. hx welchen Wert annimmt, ist dir klar, oder?

    Viele Grüße
    ~blaze~

    Ich habe die nur so geschrieben, weil sie nicht delta entsprechen müssen. hx sollte eigentlich immer delta sein. Für hv ergibt sich der Wert ja in Abhängigkeit von den Parametern, die du an die Funktion v gibst:
    Ich glaube, dass hv im Fall
    v(t + 1/2 hx, ...)

    delta/2 sein sollte (bzw. hx/2). Aber nach meinem anstrengenden Tag bin ich mir nicht mehr ganz sicher, ob das der richtige Gedanke von vorhin war. Ich glaube aber schon... ;)

    Viele Grüße
    ~blaze~

    Ich hab es eigentlich so implementiert wie du es meintest:

    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. private double distance;
    17. public CelestialBody(double x0, double y0, double vx0, double vy0)
    18. {
    19. tmpVx = vx0;
    20. tmpVy = vy0;
    22. tmpX = x0;
    23. tmpY = y0;
    25. delta = 60 * 60 * 24 * 7 * 0.1;
    26. }
    28. public void Update()
    29. {
    30. double halfdx = .5d * delta;
    31. distance = Math.Sqrt(tmpX * tmpX + tmpY * tmpY);
    33. tmpVx = vx(tmpX);
    34. tmpVy = vy(tmpY);
    36. k1 = delta * vx(tmpX);
    37. k2 = delta * vx(tmpX + halfdx * k1);
    38. k3 = delta * vx(tmpX + halfdx * k2);
    39. k4 = delta * vx(tmpX + k3);
    40. tmpX = tmpX + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    43. k1 = delta * vy(tmpY);
    44. k2 = delta * vy(tmpY + halfdx * k1);
    45. k3 = delta * vy(tmpY + halfdx * k2);
    46. k4 = delta * vy(tmpY + k3);
    47. tmpY = tmpY + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    50. Console.WriteLine(tmpX + ";" + tmpY);
    51. }
    53. private double ax( double x)
    54. {
    55. return -SUN_G * (x / Math.Pow(distance, Exponent));
    56. }
    58. private double ay(double y)
    59. {
    60. return -SUN_G * (y / Math.Pow(distance, Exponent));
    61. }
    63. private double vx(double x)
    64. {
    65. double halfdx = .5d * delta;
    66. k1 = delta * ax(tmpX);
    67. k2 = delta * ax(tmpX + halfdx * k1);
    68. k3 = delta * ax(tmpX + halfdx * k2);
    69. k4 = delta * ax(tmpX + k3);
    70. return tmpVx + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    72. }
    74. private double vy(double y)
    75. {
    76. double halfdx = .5d * delta;
    77. k1 = delta * ay(tmpY);
    78. k2 = delta * ay(tmpY + halfdx * k1);
    79. k3 = delta * ay(tmpY + halfdx * k2);
    80. k4 = delta * ay(tmpY + k3);
    81. return tmpVy + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    82. }
    84. List<PointF> pts = new List<PointF>();
    85. public void Render(Graphics g)
    86. {
    87. PointF ptf = new PointF((float)(tmpX / ASTRONOMICUNIT) * 200 + 350f, (float)(tmpY / ASTRONOMICUNIT) * 200 + 350f);
    88. g.FillEllipse(Brushes.Blue, new RectangleF(ptf.X, ptf.Y, 15, 15));
    90. if (pts.Contains(ptf))
    91. pts.Clear();
    92. pts.Add(ptf);
    94. for (int i = 0; i < pts.Count - 1; i++)
    95. {
    96. g.DrawLine(Pens.Red, pts[i], pts[i + 1]);
    97. }
    99. }
    101. }



    Was mach ich denn nun wieder falsch?
    Und Gott alleine weiß alles am allerbesten und besser.

    Hi
    mir fallen mehrere Dinge auf:
    - du hast ein sehr großes Delta und eine sehr große Geschwindigkeit --> kleinere Schritte simulieren, d.h. kleines delta. Bei großen Werten ist die Abweichung beim Runge-Kutta-Verfahren vom tatsächlichen Wert zu groß. Du könntest höher einen Runge-Kutta höherer Ordnung anwenden, das dürfte die Performance gegenüber der schrittweisen Annäherung weit verbessern.
    - du arbeitest mit den absoluten Positionen, nicht mit den relativen (falls die Sonne sich nicht dauerhaft mit unendlicher Masse an Position (0, 0) befindet --> tmpX, tmpY müssen die Relativposition zur Sonne angeben
    - x^3 ist ebenfalls x*x*x
    - arbeite unbedingt mit Vektoren, die vereinfachen dir alles. Du kannst entweder die aus System.Numerics verwenden oder eigene definieren, die aber vmtl. weniger effizient sind, da sie nicht auf Vektoroperationen laufen dürften.
    - arbeite mehr mit Parametern statt mit Variablen, wenn diese ihre Gültigkeit eh verlieren oder nicht an die Instanz gebunden sind (k1, k2, k3, k4, tmpVx, tmpVy, tmpX, tmpY bzw. delta)

    Viele Grüße
    ~blaze~

    Zusätzlich solltest du auch noch aufpassen, dass sich die Objekte nicht zu sehr annähern... Denn für verschwindene Abstände divergiert das von dir verwendete Gravitationsgesetz. Die Konsequenz dessen ist dann das "Wegschiessen" des Objektes. Also entweder brichst du dann die Simulation ab, oder du musst das ganze so umgestallten, dass die Objekte keine Punktmassen sondern ausgedehnte Objekte mit nie verschwindenen Abständen sind.
    Aber vorher solltest du auf jeden Fall, wie von ~blaze~ angesprochen, die Schrittgrösse dramatisch verringern. Nimm dir einmal die Zeit und lass das Ganze mit einer sehr, sehr kleinen Schrittgrösse simulieren und betrachte dann mal die Trajektorie. Vielleicht löst das schon alle Probleme.

    Hallo!

    Lieben Dank Euch allen.
    Wieso denn unbedingt Vektoren nutzen?
    Geht die Akkuratesse verloren wenn es in Koordinaten zerlegt wird( im Grunde läuft es ja auf dasselbe hinaus), oder eben wegen den damit einhergehenden
    Nutzungsmöglichkeiten( exemplarisch Multiplikation mit Skalar, Länge, et cetera)?


    Zum Delta: Wie klein soll es denn sein? Hab's mit 35 versucht, das dauert doch ewig... wieso kann denn das eulersche Verfahren so gut mit großen Deltawerten umgehen?
    Ist aber die Implementierung soweit in Ordnung?
    Danke
    Und Gott alleine weiß alles am allerbesten und besser.

    Du kannst auch einen Runge-Kutta höherer Ordnung nutzen. Das sollte dann eine wesentlich bessere Näherung darstellen. Du kannst auch mal einen Testlauf durchführen mit einem dynamischen Runge-Kutta und dann abgleichen, ab welchem Wert der relative Fehler zwischen dem Runge-Kutta mit (n+1)-ter Ordnung und dem mit n-ter Ordnung insignifikant klein wird. Das wäre die übliche Vorgehensweise.

    Ich denke schon, dass die Implementierung so stimmt. Mir fällt beim schnellen Drüberschauen zumindest kein weiterer Fehler auf. Ansonsten teste mal mit langsamen Geschwindigkeiten und kleinen Zeitabschnitten, ob das passt.

    Vektoren gestalten quasi alles einfacher, was du vorhast. Wenn du außerdem auch noch die aus System.Numerics verwendest, hast du die Parallelisierung noch dazu, sofern du auf dem Framework 4.6.2 arbeitest, zumindest (davor kann ich nicht beurteilen, ob das möglich ist). Ansonsten hast du den Vorteil, dass du dir Redundanz sparen kannst. Um die Vektoren kommst du so oder so bei komplexeren 2D/3D/...-Anwendungen nicht herum und sie machen's ja quasi nicht schlechter.

    Viele Grüße
    ~blaze~

    RK4 ist nicht das höchste?
    Wie sehen denn höhere Ordnungen aus?
    Wie viele "Vorberechnungen"(also k1...kn) gibt es denn?
    Ist so etwas quasi grenzenlos? ( Scheint mir so)

    Danke für Eure Hilfe.
    Und Gott alleine weiß alles am allerbesten und besser.

    Schau dir mal das hier an:
    en.wikipedia.org/wiki/Runge%E2…nge.E2.80.93Kutta_methods
    Dort findest du die Definition.

    Du kannst dir auch den impliziten Runge-Kutta weiter unten anschauen. Allerdings empfehle ich für die Lösung von diesem die Verwendung der GPU über Compute Shader/CUDA/..., wenn es sich um Echtzeitanwendungen handelt. Die Zahl der Berechnungen steigt dann doch enorm an.
    Runge-Kutta habe ich dafür aber bisher noch nicht implementiert.

    Viele Grüße
    ~blaze~

    Hallo miteinander,

    Da öffne ich VB-Paradise nach mehreren Monaten endlich wieder und was sehe ich? Einen Thread über eines meiner Lieblingsthemen: Runge-Kutta-Verfahren. :D

    Momentan bin ich noch in den Ferien, aber übernächsten Samstag komme ich zurück, dann kann ich dir mal meine kleine Arbeit über RK-Verfahren hochladen, indem in alle gängigen RK-Verfahren explizit sowie eingebettet in Matlab implementiert habe. Bis dahin kann ich nur auf ein kleines GIF aus einem Thread von mir verlinken:


    Mit freundlichen Grüssen und bis in einer Woche+,

    Higlav

    ----------------

    EDIT: Nevermind - Ich hab's in meiner Cloud gefunden. Hier der Link zur pdf Datei. Wenn du die Rohdaten willst, einfach Bescheid geben.

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

    @φConst, @~blaze~, @Higlav Ich hab mal meinen Lösungsvorschlag in den Sourcecodeaustausch gepackt.
    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!

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

    @RodFromGermany
    Lieben Dank, interessant gelöst!

    Hab jetzt durch die Beispielaufgabe von Higlavs coolem Dokument verstanden dass die Schrittweite wirklich sehr klein gewählt sein muss:

    C#-Quellcode

    1. static double x = 1;
    3. static double k1, k2, k3, k4;
    4. static double h = 0.00000001;
    5. static double halfdx = h * .5;
    6. static double time;
    7. static void Main(string[] args)
    8. {
    9. while (true)
    10. {
    11. k1 = h* func_x(time, x);
    12. k2 = h* func_x(time + halfdx, x + halfdx * k1);
    13. k3 = h* func_x(time + halfdx, x + halfdx * k2);
    14. k4 = h* func_x(time + h, x + h * k3);
    15. x = x + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    17. Console.WriteLine("Diff: " + Decimal.Parse(Math.Sqrt(Math.Pow((exact_x(time) - x ),2)) + "", System.Globalization.NumberStyles.Float));
    19. time = time+ h;
    21. }
    22. }
    24. static double func_x(double t, double a)
    25. {
    26. return -(Math.Sin(t * t * t) + 3 * (t * t * t) * Math.Cos(t * t * t)) * a;
    27. }
    29. static double exact_x(double t)
    30. {
    31. return Math.Exp(-t * Math.Sin(t * t * t));
    32. }


    Und ich nutze die ganze Zeit Schrittweite = 60480 :D :D

    Leider ist bei meiner Implementierung trotzdem irgendein Fehler, selbst bei Schrittweite 1 scheint es in die Sonne zu stürzen.. Gucke mir das implizite Runge Kutta Verfahren an.
    Und Gott alleine weiß alles am allerbesten und besser.

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

    Hallo Leute,

    das nervt.

    Hab alles probiert, gar sowas:

    C#-Quellcode

    1. public class CelestialBody
    2. {
    3. public static readonly double ASTRONOMICUNIT = 149597870700;
    4. public static readonly double GRAVITATIONCONSTANT = 6.67408e-11;
    5. public static readonly double SUN_G = 1.989E30 * GRAVITATIONCONSTANT;
    6. public static readonly double Exponent = 3;
    8. public static readonly double Delta0 = 60 * 60 * 7 * 24 * 0.001;
    9. public static readonly double Delta1 = 60 * 60 * 7 * 24 * 0.01;
    11. public Vector2 Position, v0;
    13. public double Delta, Time;
    16. public CelestialBody(double x0, double y0, double vx0, double vy0)
    17. {
    18. Position = new Vector2(x0, y0);
    19. v0 = new Vector2(vx0, vy0);
    21. }
    22. public void Update()
    23. {
    24. Vector2 k1;
    25. Vector2 k2;
    26. Vector2 k3;
    27. Vector2 k4;
    29. k1 = Delta * v(Time, Position, v0, Delta / 2.0);
    30. k2 = Delta * v(Time + .5 * Delta, Position + .5 * k1, v0, Delta / 2.0);
    31. k3 = Delta * v(Time + .5 * Delta, Position + .5 * k2, v0, Delta / 2.0);
    32. k4 = Delta * v(Time + Delta, Position + k3, v0, Delta / 2.0);
    34. Position += (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    36. Time += Delta;
    37. }
    39. public Vector2 a(Vector2 pos)
    40. {
    41. double dist = pos.Length();
    42. return new Vector2(-SUN_G * (pos.X / (dist * dist * dist)), -SUN_G * (pos.Y / (dist * dist * dist)));
    43. }
    45. public Vector2 v(double time, Vector2 pos, Vector2 init, double delta)
    46. {
    47. Vector2 k1;
    48. Vector2 k2;
    49. Vector2 k3;
    50. Vector2 k4;
    52. Vector2 velocity = init;
    54. for (double i = 0; i <= time; i += delta)
    55. {
    56. k1 = delta * a(pos);
    57. k2 = delta * a(pos + .5 * k1);
    58. k3 = delta * a(pos + .5 * k2);
    59. k4 = delta * a(pos + k3);
    61. velocity += (1 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) ;
    62. }
    64. return velocity;
    66. }
    70. List<PointF> pts = new List<PointF>();
    71. public void Render(Graphics g, Color color)
    72. {
    74. PointF ptf = new PointF((float)(Position.X / ASTRONOMICUNIT) * 200 + 350f, (float)(Position.Y / ASTRONOMICUNIT) * 200 + 350f);
    75. g.FillEllipse(Brushes.Blue, new RectangleF(ptf.X, ptf.Y, 15, 15));
    77. if (!pts.Contains(ptf))
    78. pts.Add(ptf);
    79. for (int i = 0; i < pts.Count - 1; i++)
    80. {
    81. g.DrawLine(new Pen(color, Pens.Red.Width), pts[i], pts[i + 1]);
    82. }
    83. }
    84. }


    Ich meine, dass das nicht an am Delta liegt... immerhin gibt es beim eulerschen Verfahren mit diesem Delta kein Problem, dann sollte es noch weniger mit rk4 geben.
    Woran liegt das?

    Lieben Dank!
    Und Gott alleine weiß alles am allerbesten und besser.

    Ohne die For Schleife und mit einem Delta(= 60 * 60 * 7 * 24 * 0.01) * 0.25 für die Geschwindigkeit ergeben sich quasi stabile Bahnen:



    Problematisch sind jetzt jedoch die Daten zu Peri -und Aphel: Merkur (grau) immer 0.38 AE, Erde ebenfalls fast konstant.
    Woran liegt das?

    Funktionieren Simulationen dieser Art etwa nur mit RK5 und höher?

    Liebe Grüsse.
    Und Gott alleine weiß alles am allerbesten und besser.

    Du musst auch die Position in den Zwischenschritten updaten.
    Effizienter wäre es übrigens tatsächlich, höhere Runge-Kutta-Verfahren zu verwenden.

    Edit: Eventuell wirst du übrigens hier fündig: scholarpedia.org/article/Runge-Kutta_methods
    Ich glaube, dass die Findung von Runge-Kutta-Verfahren höherer Ordnung nicht ganz trivial ist. Ggf. muss man da auf bereits vorhandene Butcher-Tableaus zurückgreifen. Ich habe jetzt zumindest beim Überfliegen nicht herausgefunden, wie genau die
    $b_i$
    -s konstituiert sein müssen (5 Minuten Recherche ist aber auch meist nicht gerade problemlösend).

    Viele Grüße
    ~blaze~

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

Facebook
Telefonspende