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 ode45
unter Matlab programmiert , um ein System von ODEs des folgenden Typs zu lösen:
wobei der Winkel des ersten Körpers zur Horizontalen ist, die Winkelgeschwindigkeit des ersten Körpers ist; ist der Winkel des zweiten Körpers in Bezug auf den ersten Körper und ist die Winkelgeschwindigkeit des zweiten Körpers. Alle Koeffizienten sind im folgenden Code in den von mir erstellten Funktionen rhs
und 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 (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 und 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:
Wenn ich die Simulation jedoch 10 Sekunden lang ausführe, bewegt sich das System unangemessen:
ODE23
Ich habe dann die Simulation mit ausgeführt, ode23
um festzustellen, ob das Problem weiterhin besteht. Am Ende habe ich dasselbe Verhalten, nur diesmal beginnt die Divergenz 1 Sekunde später:
ODE15s
Ich habe dann die Simulation mit ausgeführt, ode15s
um festzustellen, ob das Problem weiterhin besteht, und nein, das System scheint selbst während 100 Sekunden stabil zu sein:
Andererseits ode15s
ist nur erster Ordnung und zu beachten, dass es nur wenige Integrationsschritte gibt. Ich habe also eine weitere Simulation mit einer ode15s
Zeit von 10 Sekunden und einer MaxStep
Größe von , um die Genauigkeit zu erhöhen. Leider führt dies zu dem gleichen Ergebnis wie bei beiden ode45
und ode23
.
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 ode
Funktionen in Matlab?
quelle
x1
undx3
. (Fügen Sie einen trockenen Kommentar zu Diagrammen ohne Legenden oder Beschreibungen ein.) Zeichnen Sie die Logarithmen von (die absoluten Werte von)x2
undx4
.Antworten:
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
In Ihrem Code können Sie dies durch Berechnen testen
was 6.191e-16 ergibt - also fast aber nicht genau null. Wie wirkt sich das auf die Dynamik Ihres Systems aus?
Darüber hinaus sieht die Lösung Ihres Systems in sehr kurzer Zeit aus wie die Lösung des linearisierten Systems
rhs
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:
so dass Ihre Gleichung wird
Jetzt brauchen wir noch einen letzten Schritt: Wir können eine Eigenwertzerlegung des Jacobi so berechnen, dass
quelle
quelle
sin(0)
cos(0)
sin(pi/2)
cos(pi/2)
rhs(t,[0,0,0 0] == [0,0,0,0]
Schauen Sie sich die Komponenten der in Ihren Funktionen berechneten Kräfte an.
quelle
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:
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.
quelle
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.
quelle
ode45
berechne, 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