[Pendel-Simulation] Was mache ich falsch?

Es gibt 3 Antworten in diesem Thema. Der letzte Beitrag () ist von RodFromGermany.

    [Pendel-Simulation] Was mache ich falsch?

    Hallo,

    ich versuche die Dynamiken eines simplen Pendels mittels RK4/Euler-Verfahren zu simulieren (wie es analytisch geht, ist mir relativ bekannt).

    Das Vektor-Diagramm sollte doch eigentlich so aussehen:



    Erklärung: Auf den "Bob"(= Kugel) wirken sowohl die Gravitations -als auch "Zug"kraft vom "Bob" zum Seil.
    Der Winkel "theta" ist nichts anderes als das Skalarprodukt der normierten Vektoren g und d.

    Man stelle sich vor der Bob befinde sich in den Extremstellen (genau dann wenn exemplarisch die Y-Koordinate des Bobs und des Zentrums übereinstimmen).
    Das die Formel (a = g + (-d * |g| * theta)) Sinn macht ist offensichtlich: Wenn der Bob sich in jener Extremstelle befindet, liegen die Vektoren g und d orthogonal zueinander:
    Das Skalarprodukt ist also 0.

    Setzt man das in die Formel ein
    a = g + (-d * 0 * theta) = g dann wirkt in jener Extremstelle ausschließlich die Gravitationskraft, was ja auch richtig sein sollte.

    Dennoch "reißt" der Bob in der Simulation einfach ab... Habe es mit verschiedenen Delta-Werten versucht: Stets ist die Geschwindigkeit in Y Richtung einfach zu hoch, sodass der Bob gar nicht hochgezogen wird.
    Ich mache also irgendwo was falsch.

    Wo und vor allem was?

    Hier der Source-Code:

    Die Form:
    Spoiler anzeigen

    C#-Quellcode

    1. public partial class Form1 : Form
    2. {
    3. public Form1()
    4. {
    5. InitializeComponent();
    6. }
    8. Vector2 center = new Vector2(500, 250);
    10. Vector2 bobPosition = new Vector2(550, 250);
    11. Vector2 bobVelocity = new Vector2(0, 0);
    13. IApproximation positionIntegrator;
    14. IApproximation velocityIntegrator;
    16. Vector2 gravitationalForce = new Vector2(0, 10);
    18. private void Form1_Load(object sender, EventArgs e)
    19. {
    20. this.Size = new Size(1000, 1000);
    22. this.DoubleBuffered = true;
    23. this.Paint += new PaintEventHandler(draw);
    24. Application.Idle += new EventHandler((a, b) =>
    25. {
    27. bobPosition = positionIntegrator.Run();
    29. Console.WriteLine(velocityIntegrator.Output);
    31. this.Invalidate();
    32. });
    34. positionIntegrator = new RK4(bobPosition, 0.005f);
    35. velocityIntegrator = new RK4(bobVelocity, 0.005f);
    37. positionIntegrator.f = pPrime;
    38. velocityIntegrator.f = pPrimePrime;
    40. }
    42. Vector2 pPrime(float dt, Vector2 y)
    43. {
    44. return velocityIntegrator.Run();
    45. }
    47. Vector2 pPrimePrime(float dt, Vector2 y)
    48. {
    49. Vector2 dir = bobPosition - center;
    50. dir.Normalize();
    51. Vector2 force = new Vector2(gravitationalForce.X, gravitationalForce.Y);
    52. force.Normalize();
    54. float forceMagnitude = gravitationalForce.Length();
    55. float theta = Vector2.Dot(dir, force);
    57. dir = -dir;
    59. Vector2 acceleration = gravitationalForce + dir * theta * (forceMagnitude);
    62. return acceleration;
    63. }
    65. private void draw( object sender, PaintEventArgs e)
    66. {
    68. Graphics g = e.Graphics;
    70. g.FillEllipse(Brushes.Red, center.X - 25f, center.Y - 25f, 50, 50);
    71. g.FillEllipse(Brushes.Blue, bobPosition.X - 12.5f, bobPosition.Y - 12.5f, 25, 25);
    73. g.DrawLine(Pens.Black, center.X, center.Y, bobPosition.X, bobPosition.Y);
    75. }
    76. }



    Hier die numerische Integrator-Klasse:

    IApproximation:
    Spoiler anzeigen

    C#-Quellcode

    1. public abstract class IApproximation
    2. {
    3. public delegate Vector2 f_delegate(float x, Vector2 y);
    5. public f_delegate f { get; set; }
    7. public float Delta { get; set; }
    8. public float Time { get; set; }
    10. public Vector2 Output { get; set; }
    12. public abstract Vector2 Run();
    13. }


    RK4:
    Spoiler anzeigen

    C#-Quellcode

    1. public class RK4 : IApproximation
    2. {
    4. public RK4(Vector2 _0, float delta)
    5. {
    6. Output = _0;
    7. Delta = delta;
    8. }
    10. public override Vector2 Run()
    11. {
    12. Vector2 k1 = Delta * f(Time, Output);
    13. Vector2 k2 = Delta * f(Time + .5f * Delta, Output + .5f * k1);
    14. Vector2 k3 = Delta * f(Time + .5f * Delta, Output + .5f * k2);
    15. Vector2 k4 = Delta * f(Time + Delta, Output + k3);
    17. Output += 1.0f / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4);
    19. Time += Delta;
    21. return Output;
    22. }
    23. }



    Liebe Grüße!

    EDIT: Okay, wenigstens eine kleine Besserung.
    Benutzt man stattdessen folgende Formel:

    a = (1 - |theta|) * g + (-d * |g| * |theta|) reißt er wenigstens nicht ab.
    Eine Pendel-Simulation ist das dennoch nicht, eher eine "Sprungfeder" -Simulation...

    _
    Und Gott alleine weiß alles am allerbesten und besser.

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

    @φConst Ich hab Deinen Code in ein neues Projekt gepackt und Vector2 durch meine eigene Vektor-Klasse ersetzt.
    Es pendelt permanent, aber der "Faden" wird immer länger.
    Spoiler anzeigen

    C#-Quellcode

    1. using System;
    3. namespace WindowsFormsApplication1
    4. {
    5. public class Vector2
    6. {
    7. public float X { get; set; }
    8. public float Y { get; set; }
    9. public Vector2(Vector2 vector)
    10. {
    11. this.X = vector.X;
    12. this.Y = vector.Y;
    13. }
    15. public Vector2(float x, float y)
    16. {
    17. this.X = x;
    18. this.Y = y;
    19. }
    21. public Vector2()
    22. {
    23. this.X = 0;
    24. this.Y = 0;
    25. }
    27. public static Vector2 operator *(float factor, Vector2 vector)
    28. {
    29. Vector2 v = new Vector2(vector);
    30. v.X *= factor;
    31. v.Y *= factor;
    32. return v;
    33. }
    35. public static Vector2 operator *(Vector2 vector, float factor)
    36. {
    37. Vector2 v = new Vector2(vector);
    38. v.X *= factor;
    39. v.Y *= factor;
    40. return v;
    41. }
    42. public static Vector2 operator +(Vector2 v1, Vector2 v2)
    43. {
    44. Vector2 v = new Vector2(v1);
    45. v.X += v2.X;
    46. v.Y += v2.Y;
    47. return v;
    48. }
    49. public static Vector2 operator -(Vector2 v1, Vector2 v2)
    50. {
    51. Vector2 v = new Vector2(v1);
    52. v.X -= v2.X;
    53. v.Y -= v2.Y;
    54. return v;
    55. }
    57. public static Vector2 operator -(Vector2 v1)
    58. {
    59. return new Vector2(-v1.X, -v1.Y);
    60. }
    62. /// <summary>
    63. /// Vektor dividiert durch Skalar
    64. /// </summary>
    65. /// <param name="d"></param>
    66. /// <param name="v"></param>
    67. /// <returns></returns>
    68. public static Vector2 operator /(float d, Vector2 v)
    69. {
    70. return new Vector2(v.X / d, v.Y / d);
    71. }
    73. /// <summary>
    74. /// Vektor dividiert durch Skalar
    75. /// </summary>
    76. /// <param name="v"></param>
    77. /// <param name="d"></param>
    78. /// <returns></returns>
    79. public static Vector2 operator /(Vector2 v, float d)
    80. {
    81. return new Vector2(v.X / d, v.Y / d);
    82. }
    84. /// <summary>
    85. /// Berechnung der Länge des Vektors
    86. /// </summary>
    87. /// <returns>die Länge</returns>
    88. public float Length()
    89. {
    90. return (float)(Math.Sqrt(this.X * this.X + this.Y * this.Y));
    91. }
    93. internal void Normalize()
    94. {
    95. float bet = this.Length();
    96. this.X /= this.Length();
    97. this.Y /= this.Length();
    98. }
    100. /// <summary>
    101. /// Skalarprodukt
    102. /// </summary>
    103. /// <param name="v1"></param>
    104. /// <param name="v2"></param>
    105. /// <returns></returns>
    106. internal static float Dot(Vector2 v1, Vector2 v2)
    107. {
    108. return v1.X * v2.X + v1.Y * v2.Y;
    109. }
    110. }
    111. }
    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!

    RodFromGermany schrieb:

    Es pendelt permanent, aber der "Faden" wird immer länger.


    Eben, genau das ist das Problem.. ^^
    Ich habe mir mal die Werte des Beschleunigungsvektors ausgegeben: Sind - ohne den "Fix" im Edit - stets positiv: Genau das ist ja falsch.
    Aber wie sonst würde man denn den Beschleunigungsvektor konstruieren?
    Und Gott alleine weiß alles am allerbesten und besser.

    @φConst Ohne jetzt alles nachvollziehen zu wollen, schau mal hier rein:
    Runge Kutta Verfahren 4. Ordnung
    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!

Facebook
Telefonspende