Pourquoi la solution numérique d'une ODE s'éloigne-t-elle d'un équilibre instable?

15

Je souhaite simuler le comportement d'un système à double pendule. Le système est un robot manipulateur à 2 degrés de liberté qui n'est pas actionné et se comportera donc principalement comme un double pendule affecté par la gravité. La seule différence principale avec un pendule double est qu'il est composé de deux corps rigides avec des propriétés de masse et d'inertie à leurs centres de masse.

En gros, j'ai programmé ode45sous Matlab pour résoudre un système d'ODE du type suivant:

[10000M110M1200100M120M22][X˙1X˙2X˙3X˙4]=[X2-V1-g1X4-V2-g2]

x1 est l'angle du premier corps par rapport à l'horizontale, x2 est la vitesse angulaire du premier corps; x3 est l'angle du deuxième corps par rapport au premier corps, et x4 est la vitesse angulaire du deuxième corps. Tous les coefficients sont spécifiés dans le code suivant, dans les fonctions rhset que fMassj'ai créées.

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

Remarquez comment j'ai défini la condition initiale de x1 (angle du premier corps par rapport à l'horizontale) pour que le système démarre dans une position complètement verticale. De cette façon, puisque seule la gravité agit, le résultat évident est que le système ne doit pas du tout bouger de cette position.

REMARQUE: dans tous les graphiques ci-dessous, j'ai tracé les solutions x1 et x3 en fonction du temps.

ODE45

Lorsque j'exécute la simulation pendant 6 secondes avec ode45, j'obtiens la solution attendue sans aucun problème, le système reste où il est et ne bouge pas:

entrez la description de l'image ici

Cependant, lorsque j'exécute la simulation pendant 10 secondes, le système commence à se déplacer de manière déraisonnable:

entrez la description de l'image ici

ODE23

J'ai ensuite exécuté la simulation avec ode23pour voir si le problème persistait. Je me retrouve avec le même comportement, mais cette fois la divergence commence 1 seconde plus tard:

entrez la description de l'image ici

ODE15s

J'ai ensuite exécuté la simulation avec ode15spour voir si le problème persistait et non, le système semble être stable même pendant 100 secondes:

entrez la description de l'image ici

Là encore, ce ode15sn'est que du premier ordre et notez qu'il n'y a que quelques étapes d'intégration. J'ai donc exécuté une autre simulation avec ode15spendant 10 secondes mais une MaxSteptaille de 0.01 pour augmenter la précision, et malheureusement, cela conduit au même résultat qu'avec les deux ode45et ode23.

entrez la description de l'image ici

Normalement, le résultat évident de ces simulations serait que le système reste à sa position initiale car rien ne le perturbe. Pourquoi cette divergence se produit-elle? Cela a-t-il quelque chose à voir avec le fait que ces types de systèmes sont de nature chaotique? Est-ce un comportement normal pour les odefonctions de Matlab?

jrojasqu
la source
Outre les équations, je pense que le schéma a aiderait également beaucoup à comprendre la question.
nicoguaro
Si vous pensez que c'est approprié, vous pouvez accepter l'une des réponses (il y a un bouton vert).
Ertxiem - réintègre Monica le
Vous ne dites pas, mais vous semblez comploter x1et x3. (Insérez un commentaire sec sur les graphiques sans légendes ni descriptions.) Essayez de tracer les logarithmes de (les valeurs absolues de) x2et x4.
Eric Towers

Réponses:

15

Je pense que les deux points principaux ont déjà été soulevés par Brian et Ertxiem: votre valeur initiale est un équilibre instable et le fait que vos calculs numériques ne soient jamais vraiment exacts fournit la petite perturbation qui provoquera l'instabilité.

Pour donner un peu plus de détails sur la façon dont cela se passe, considérez votre problème sous la forme d'un problème général de valeur initiale

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~

Dans votre code, vous pouvez le tester en calculant

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

ce qui donne 6.191e-16 - donc presque mais pas exactement zéro. Comment cela affecte-t-il la dynamique de votre système?

fy0y~0

De plus, en très peu de temps, la solution de votre système ressemble à la solution du système linéarisé

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).

Je ne pouvais pas être dérangé pour calculer le jacobien à la main, j'ai donc utilisé la différenciation automatique pour obtenir une bonne approximation:

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

pour que votre équation devienne

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

Maintenant, nous avons besoin d'une dernière étape: nous pouvons calculer une décomposition de valeurs propres du jacobien telle que

J=QDQ1

DJQde(t):=Q1d(t)

e˙(t)=e(t),e(0)=Q-10.

e˙je(t)=λjeeje(t),eje(0)=je-e composante de Q-10

je=1,2,3,4λ1=3.2485

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

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

Daniel
la source
16

π/2π/2

Brian Borchers
la source
4
Si vous surveillez attentivement les variables d'état (en regardant les valeurs imprimées en notation scientifique), vous devriez être en mesure de voir le mouvement très lent initial de s'éloigner de l'équilibre.
Brian Borchers
Cela a du sens et en effet, lorsque je démarre le système dans une position verticale descendante (étant un point d'équilibre stable), le système ne bouge pas du tout, du moins pour une simulation de 1000 secondes que je considère comme une très longue période de temps .
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

Regardez les composantes des forces calculées dans vos fonctions.

π

1016

a=1.0a=a+1016

alephzero
la source
4

L'hypothèse initiale était que la position initiale était à un équilibre stable (c'est-à-dire un minimum d'énergie potentielle) avec une énergie cinétique nulle et le système a commencé à s'éloigner de l'équilibre.
Puisque physiquement cela ne peut pas arriver (si l'on considère la mécanique classique), deux choses me sont venues à l'esprit:

  1. π/2π/2

  2. La seconde est qu'il y a peut-être quelque chose qui ne va pas avec les équations du mouvement (peut-être une faute de frappe quelque part?). Pouvez-vous s'il vous plaît écrire les équations explicitement? Vous pourriez peut-être tracer l'accélération angulaire en fonction de la position initiale de chaque pendule, en supposant une vitesse angulaire nulle pour vérifier s'il y a quelque chose de bizarre.

Ertxiem - réintégrer Monica
la source
1
π
2
Soit dit en passant, juste pour le plaisir, si vous vouliez garder le système à la position verticale instable, vous pourriez changer votre origine des coordonnées pour que l'angle égal à zéro pointe vers le haut.
Ertxiem - réintègre Monica le
@Ertxiem une autre option est d'introduire un petit frottement dans les broches qui mangerait des erreurs numériques.
svavil
sin(π)
Étant donné que cela ne peut pas se produire physiquement - Étant donné que nous sommes à un équilibre instable, je conteste quelque peu cela. Les systèmes physiques (sans trop de friction) ne restent pas dans des équilibres instables. Plus généralement, si vous simulez des systèmes réels, vous voudrez éviter qu'ils ne se coincent accidentellement dans un équilibre instable (quelle qu'en soit la raison) - c'est une fonctionnalité, pas un bug. (Il y a de rares exceptions à cela, comme l'état non infecté en immunologie, qui est un équilibre instable qui peut être maintenu.)
Wrzlprmft
0

Vous devriez chercher plus sur les doubles pendules: ils sont ce que nous appelons des "systèmes chaotiques". Même si elles se comportent selon des règles simples, à partir de conditions initiales légèrement différentes, les solutions divergent assez rapidement. Faire des simulations numériques pour ce type de systèmes n'est pas facile. Jetez un œil à la vidéo suivante pour mieux comprendre le problème.

Pour le pendule simple ou double, vous pouvez écrire une formule pour l'énergie totale du système. En supposant que les forces de friction sont négligées, cette énergie totale est préservée par le système analytique. Numériquement, c'est un tout autre problème.

Avant d'essayer le pendule double, essayez le pendule simple. Vous remarquerez que pour les méthodes Runge-Kutta de petit ordre, l'énergie du système augmentera dans les simulations numériques, au lieu de rester constante (c'est ce qui se passe dans vos simulations: vous obtenez le mouvement à partir de rien). Pour éviter cela, des méthodes RK d'ordre supérieur pourraient être utilisées (l'ode45 est d'ordre 4; le RK d'ordre 8 fonctionnerait mieux). Il existe d'autres méthodes appelées "méthodes symplectiques" qui sont conçues de telle sorte que les simulations numériques conservent l'énergie. En général, vous devez arrêter la simulation dès que l'énergie augmente considérablement par rapport à votre initialisation.

Et essayez de comprendre le pendule simple avant de passer au double.

Beni Bogosel
la source
2
Il ne s'agit pas pour autant d'un système chaotique. Vous pouvez également avoir un équilibre instable dans des systèmes non chaotiques, par exemple un pendule unique étant "sur sa tête", et il présentera le même comportement décrit dans la question.
Daniel
1
Il n'est pas vrai non plus que l'énergie augmente pour RKM de petit ordre: Euler implicite est de premier ordre et montre exactement le comportement opposé.
Daniel
@BeniBogosel Vous mentionnez des méthodes symplectiques qui ont retenu mon attention car évidemment, dans mon exemple, l'énergie n'est pas conservée. Cependant, pourriez-vous indiquer une méthode symplectique spécifique qui pourrait être mise en œuvre ici?
jrojasqu
@jrojasqu pourquoi dites-vous que l'énergie n'est pas conservée sur votre système?
Ertxiem - réintègre Monica le
ode45π , etc.).
jrojasqu