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

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:

$$\left[ \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & M_{11} & 0 & M_{12}\\ 0 & 0 & 1 & 0\\ 0 & M_{12} & 0 & M_{22} \end{array} \right] \left[ \begin{array}{c} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3\\ \dot{x}_4 \end{array} \right]= \left[ \begin{array}{c} x_2\\ -V_1-G_1\\ x_4\\ -V_2-G_2 \end{array} \right]$$

où $x_1$ est l'angle du premier corps par rapport à l'horizontale, $x_2$ est la vitesse angulaire du premier corps; $x_3$ est l'angle du deuxième corps par rapport au premier corps, et $x_4$ 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 $x_1$ (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 $x_1$ et $x_3$ 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:

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

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:

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:

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.

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?

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

$$\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{M}^{-1} \mathbf{f}(t, \mathbf{y}(t)) \end{equation}$$ $\mathbf{y}(t) = (x_1(t), x_2(t), x_3(t), x_4(t))$

$$\begin{equation} \mathbf{f}(t, \mathbf{y}(t)) = \begin{bmatrix} x_2 \\ -V_1 - G_1 \\ x_4 \\ -V_2 - G_2 \end{bmatrix} \end{equation}$$

$\mathbf{f}(0, \mathbf{y}_0) = 0$$\dot{\mathbf{y}}(0) = 0$$\tilde{\mathbf{f}}$

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?

$\mathbf{f}$$\mathbf{y}_0$$\tilde{\mathbf{y}}_0$

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

$$\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{f}(0, \mathbf{y}_0) + \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) = \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) \end{equation}$$

$\mathbf{f}'$$\mathbf{f}$rhs$\mathbf{y}_0$$\mathbf{d}(t) := \mathbf{y}(t) - \mathbf{y}_0$$\mathbf{d}$

$$\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{f}'(0, \mathbf{y}_0) \mathbf{d}(t). \end{equation}$$

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:

$$\begin{equation} \mathbf{J} := \mathbf{f}'(0, \mathbf{y}_0) = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 9.81 & 0 & 2.4525 & 0 \\ 0 & 0 & 0 & 1 \\ 2.4525 & 0 & 2.4525 & 0 \end{bmatrix} \end{equation}$$

pour que votre équation devienne

$$\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{J} \mathbf{d}(t), \mathbf{d}(0) = \mathbf{\tilde{y}}_0 - \mathbf{y}_0 \end{equation}$$

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

$$\begin{equation} \mathbf{J} = \mathbf{Q} \mathbf{D} \mathbf{Q}^{-1} \end{equation}$$

$\mathbf{D}$$\mathbf{J}$$\mathbf{Q}$$\mathbf{d}$$\mathbf{e}(t) := \mathbf{Q}^{-1} \mathbf{d}(t)$

$$\begin{equation} \dot{\mathbf{e}}(t) = \mathbf{D} \mathbf{e}(t), \mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}_0. \end{equation}$$

$\mathbf{D}$

$$\begin{equation} \dot{e}_i(t) = \lambda_i e_i(t), e_i(0) = i-\textrm{th component of}~\mathbf{Q}^{-1} \mathbf{d}_0 \end{equation}$$

$i=1,2,3,4$$\lambda_1 = 3.2485$

$$\begin{equation} e_1(t) = e_1(0) e^{3.2485 t}. \end{equation}$$

$\mathbf{d}(0) = 0$$\mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}(0) = 0$$e_1(0) = 0$$e_1(0)$

$\pi/2$$\pi/2$

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

$\pi$

$10^{-16}$

$a = 1.0$$a = a + 10^{-16}$

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. $\pi/2$$-\pi/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.

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.