Warum verlässt die numerische Lösung einer ODE ein instabiles Gleichgewicht?

15

Ich möchte das Verhalten eines Doppelpendelsystems simulieren. Das System ist ein Roboter-Manipulator mit 2 Freiheitsgraden, der nicht betätigt wird und sich daher hauptsächlich wie ein durch die Schwerkraft beeinflusstes Doppelpendel verhält. Der einzige Hauptunterschied bei einem Doppelpendel besteht darin, dass es sich aus zwei starren Körpern zusammensetzt, deren Massen- und Trägheitseigenschaften in ihren Massenschwerpunkten liegen.

Grundsätzlich habe ich ode45unter Matlab programmiert , um ein System von ODEs des folgenden Typs zu lösen:

[10000M110M1200100M120M22][x˙1x˙2x˙3x˙4]=[x2V1G1x4V2G2]

wobei x1 der Winkel des ersten Körpers zur Horizontalen ist, x2 die Winkelgeschwindigkeit des ersten Körpers ist; x3 ist der Winkel des zweiten Körpers in Bezug auf den ersten Körper und x4 ist die Winkelgeschwindigkeit des zweiten Körpers. Alle Koeffizienten sind im folgenden Code in den von mir erstellten Funktionen rhsund angegeben fMass.

clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);

function F=rhs(t,x)
    m=[1 1];
    l=0.5;
    a=[0.25 0.25];
    g=9.81;
    c1=cos(x(1));
    s2=sin(x(3));
    c12=cos(x(1)+x(3));
    n1=m(2)*a(2)*l;
    V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
    V2=n1*s2*x(2)^2;
    G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
    G2=m(2)*g*a(2)*c12;

    F(1)=x(2);
    F(2)=-V1-G1;
    F(3)=x(4);
    F(4)=-V2-G2;
    F=F';     
end

function M=fMass(t,x)
    m=[1 1];
    l=0.5;
    Izz=[0.11 0.11];
    a=[0.25 0.25];
    c2=cos(x(3));
    n1=m(2)*a(2)*l;
    M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
    M12=m(2)*a(2)^2+n1*c2+Izz(2);
    M22=m(2)*a(2)^2+Izz(2);
    M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end

Beachten Sie, wie ich die Anfangsbedingung von x1 (Winkel des ersten Körpers in Bezug auf die Horizontale) so einstelle , dass das System in einer vollständig vertikalen Position startet. Auf diese Weise ist das offensichtliche Ergebnis, dass sich das System von dieser Position aus überhaupt nicht bewegen sollte, da nur die Schwerkraft wirkt.

HINWEIS: In allen folgenden Grafiken habe ich die Lösungen x1 und x3 in Bezug auf die Zeit aufgetragen .

ODE45

Wenn ich die Simulation für 6 Sekunden mit laufe ode45, erhalte ich die erwartete Lösung ohne Probleme, das System bleibt dort, wo es ist und bewegt sich nicht:

Bildbeschreibung hier eingeben

Wenn ich die Simulation jedoch 10 Sekunden lang ausführe, bewegt sich das System unangemessen:

Bildbeschreibung hier eingeben

ODE23

Ich habe dann die Simulation mit ausgeführt, ode23um festzustellen, ob das Problem weiterhin besteht. Am Ende habe ich dasselbe Verhalten, nur diesmal beginnt die Divergenz 1 Sekunde später:

Bildbeschreibung hier eingeben

ODE15s

Ich habe dann die Simulation mit ausgeführt, ode15sum festzustellen, ob das Problem weiterhin besteht, und nein, das System scheint selbst während 100 Sekunden stabil zu sein:

Bildbeschreibung hier eingeben

Andererseits ode15sist nur erster Ordnung und zu beachten, dass es nur wenige Integrationsschritte gibt. Ich habe also eine weitere Simulation mit einer ode15sZeit von 10 Sekunden und einer MaxStepGröße von 0.01 , um die Genauigkeit zu erhöhen. Leider führt dies zu dem gleichen Ergebnis wie bei beiden ode45und ode23.

Bildbeschreibung hier eingeben

Normalerweise würde das offensichtliche Ergebnis dieser Simulationen sein, dass das System an seiner ursprünglichen Position bleibt, da es durch nichts gestört wird. Warum tritt diese Divergenz auf? Hat dies damit zu tun, dass diese Art von Systemen chaotischer Natur ist? Ist dies ein normales Verhalten für odeFunktionen in Matlab?

jrojasqu
quelle
Neben den Gleichungen, denke ich, würde das Schema auch viel helfen, die Frage zu verstehen.
Nicoguaro
Wenn Sie es für angemessen halten, akzeptieren Sie möglicherweise eine der Antworten (es gibt einen grünen Knopf).
Ertxiem - wieder Monica
Sie sagen nicht, aber Sie scheinen zu planen x1und x3. (Fügen Sie einen trockenen Kommentar zu Diagrammen ohne Legenden oder Beschreibungen ein.) Zeichnen Sie die Logarithmen von (die absoluten Werte von) x2und x4.
Eric Towers

Antworten:

15

Ich denke, die beiden Hauptpunkte wurden bereits von Brian und Ertxiem angesprochen: Ihr Anfangswert ist ein instabiles Gleichgewicht, und die Tatsache, dass Ihre numerischen Berechnungen niemals wirklich genau sind, sorgt für die kleine Störung, die die Instabilität auslöst.

Betrachten Sie Ihr Problem in Form eines allgemeinen Anfangswertproblems, um etwas detaillierter zu beschreiben, wie sich dies entwickelt

y˙(t)=M1f(t,y(t))
y(t)=(x1(t),x2(t),x3(t),x4(t))

f(t,y(t))=[x2V1G1x4V2G2]

f(0,y0)=0y˙(0)=0f~

In Ihrem Code können Sie dies durch Berechnen testen

norm(rhs(0,[pi/2 0 0 0]))

was 6.191e-16 ergibt - also fast aber nicht genau null. Wie wirkt sich das auf die Dynamik Ihres Systems aus?

fy0y~0

Darüber hinaus sieht die Lösung Ihres Systems in sehr kurzer Zeit aus wie die Lösung des linearisierten Systems

y˙(t)=f(0,y0)+f(0,y0)(y(t)y0)=f(0,y0)(y(t)y0)

ffrhsy0d(t):=y(t)y0d

d˙(t)=f(0,y0)d(t).

Ich konnte mich nicht die Mühe machen, den Jacobian von Hand zu berechnen, und benutzte die automatische Differenzierung , um eine gute Annäherung zu erhalten:

J:=f(0,y0)=[01009.8102.4525000012.452502.45250]

so dass Ihre Gleichung wird

d˙(t)=Jd(t),d(0)=y~0y0

Jetzt brauchen wir noch einen letzten Schritt: Wir können eine Eigenwertzerlegung des Jacobi so berechnen, dass

J=QDQ1

DJQ.de(t): =Q.-1d(t)

e˙(t)=De(t),e(0)=Q.-1d0.

D

e˙ich(t)=λicheich(t),eich(0)=ich-Komponente von Q.-1d0

ich=1,2,3,4λ1=3,2485

e1(t)=e1(0)e3.2485t.

d(0)=0e(0)=Q1d(0)=0e1(0)=0e1(0)

Daniel
quelle
16

π/2π/2

Brian Borchers
quelle
4
Wenn Sie die Zustandsvariablen sorgfältig überwachen (indem Sie die in wissenschaftlicher Notation ausgedruckten Werte betrachten), sollten Sie in der Lage sein, die anfänglich sehr langsame Bewegung vom Gleichgewicht weg zu erkennen.
Brian Borchers
Dies ist sinnvoll, und wenn ich das System in einer senkrechten Position nach unten starte (da dies ein Punkt des stabilen Gleichgewichts ist), bewegt sich das System zumindest für eine Simulation von 1000 Sekunden, die ich für einen sehr langen Zeitraum halte, überhaupt nicht .
Jrojasqu
6
x1sin(0)cos(0)sin(pi/2)cos(pi/2)rhs(t,[0,0,0 0] == [0,0,0,0]
π/2
1
θ=0 1016
4

Schauen Sie sich die Komponenten der in Ihren Funktionen berechneten Kräfte an.

π

1016

a=1.0a=a+1016

Alephzero
quelle
4

Die anfängliche Annahme war, dass sich die Ausgangsposition in einem stabilen Gleichgewicht (dh einem Minimum der potentiellen Energie) mit null kinetischer Energie befand und das System begann, sich aus dem Gleichgewicht zu entfernen.
Da es physikalisch nicht passieren kann (wenn wir die klassische Mechanik betrachten), kamen mir zwei Dinge in den Sinn:

  1. π/2π/2

  2. Die zweite ist, dass vielleicht etwas mit den Bewegungsgleichungen nicht stimmt (vielleicht irgendwo ein Tippfehler?). Können Sie die Gleichungen bitte explizit schreiben? Vielleicht können Sie die Winkelbeschleunigung als Funktion der Ausgangsposition jedes Pendels darstellen, wobei Sie eine Winkelgeschwindigkeit von Null voraussetzen, um zu prüfen, ob etwas Unheimliches vorliegt.

Ertxiem - Setzen Sie Monica wieder ein
quelle
1
π
2
Übrigens, nur zum Spaß, wenn Sie das System in einer instabilen vertikalen Position halten möchten, können Sie den Koordinatenursprung so ändern, dass der Winkel Null nach oben zeigt.
Ertxiem - wieder Monica
@Ertxiem eine andere Möglichkeit ist, eine kleine Reibung in die Stifte einzuführen, die numerische Fehler verursachen würde.
Svavil
sin(π)
Da es physisch nicht passieren kann - Angesichts der Erkenntnis, dass wir uns in einem instabilen Gleichgewicht befinden, fordere ich dies etwas heraus. Physikalische Systeme (ohne zu große Reibung) bleiben nicht in einem instabilen Gleichgewicht. Wenn Sie reale Systeme simulieren, möchten Sie im Allgemeinen vermeiden, dass sie versehentlich in einem instabilen Gleichgewicht stecken bleiben (wie auch immer) - es ist eine Funktion, kein Fehler. (Es gibt einige seltene Ausnahmen, wie den nicht infizierten Zustand in der Immunologie, der ein instabiles Gleichgewicht darstellt, das aufrechterhalten werden kann.)
Wrzlprmft
0

Sie sollten mehr über Doppelpendel suchen: Sie sind das, was wir "chaotische Systeme" nennen. Obwohl sie sich nach einfachen Regeln verhalten, ausgehend von leicht unterschiedlichen Anfangsbedingungen, weichen die Lösungen recht schnell voneinander ab. Numerische Simulationen für diese Art von Systemen durchzuführen ist nicht einfach. Schauen Sie sich das folgende Video an, um mehr über das Problem zu erfahren.

Für das einfache oder doppelte Pendel können Sie eine Formel für die Gesamtenergie des Systems schreiben. Unter der Annahme, dass Reibungskräfte vernachlässigt werden, wird diese Gesamtenergie vom Analysesystem erhalten. Numerisch ist dies ein ganz anderes Problem.

Bevor Sie das Doppelpendel ausprobieren, probieren Sie das einfache Pendel aus. Sie werden feststellen, dass bei Runge-Kutta-Methoden kleiner Ordnung die Energie des Systems in den numerischen Simulationen zunimmt, anstatt konstant zu bleiben (dies passiert in Ihren Simulationen: Sie bewegen sich aus dem Nichts heraus). Um dies zu verhindern, könnten RK-Methoden höherer Ordnung verwendet werden (ode45 hat die Ordnung 4; RK der Ordnung 8 würde besser funktionieren). Es gibt andere als "symplektische Methoden" bezeichnete Methoden, die so ausgelegt sind, dass die numerischen Simulationen die Energie sparen. Im Allgemeinen sollten Sie die Simulation stoppen, sobald die Energie im Vergleich zu Ihrer Initialisierung erheblich ansteigt.

Und versuchen Sie, das einfache Pendel zu verstehen, bevor Sie zum doppelten Pendel übergehen.

Beni Bogosel
quelle
2
Dies ist jedoch keine Frage des chaotischen Systems. Sie können auch in nicht-chaotischen Systemen ein instabiles Gleichgewicht haben, z. B. wenn ein einzelnes Pendel "auf dem Kopf" steht und dasselbe Verhalten aufweist, wie in der Frage beschrieben.
Daniel
1
Es ist auch nicht wahr, dass die Energie für RKM kleiner Ordnung zunimmt: impliziter Euler ist erster Ordnung und zeigt genau das Gegenteil.
Daniel
@BeniBogosel Sie erwähnen symplektische Methoden, die meine Aufmerksamkeit auf sich gezogen haben, weil in meinem Beispiel die Energie offensichtlich nicht erhalten bleibt. Könnten Sie jedoch eine bestimmte symplektische Methode angeben, die hier implementiert werden könnte?
Jrojasqu
@jrojasqu warum sagst du, dass Energie auf deinem System nicht konserviert wird?
Ertxiem - wieder Monica
@Ertxiem Wenn ich die gesamte mechanische Energie des Systems (kinetische + potenzielle Energien) mit den von bereitgestellten Ausgaben ode45berechne, erhalte ich einen Wert, der mit Null beginnt und dann mit der Zeit wächst. Der Wert ist in den ersten Sekunden sehr sehr klein, wächst aber immer noch von Null weg. Ich glaube, es liegt an den Problemen, die in den obigen Antworten angesprochen wurden (Annäherung anπ, etc.).
Jrojasqu