Système d'EDOs: Runge-Kutta 4
Bonsoir,
A titre personnel, je souhaiterais implémenter l'algorithme de Runge-Kutta d'ordre $4$ afin de résoudre un système d'équations différentielles couplées.
Pour arriver à cela, je commence par résoudre le problème du pendule circulaire (équation différentielle ordinaire d'ordre $2$):
$$\frac{d^2 \theta}{dt^2} = -g \sin(\theta)$$
Afin de résoudre numériquement ce problème, je pose
$$\begin{cases} \frac{d \theta}{dt} = f(t, \theta, \dot{\theta}) \\ \frac{df}{dt} = \frac{- g}{l} \sin(\theta) = g(t, \theta, \dot{\theta}) \end{cases}$$
Les conditions initiales sont: $$\theta(0) = \theta_0 \textrm{ et } \dot{\theta}_0 = 0$$
J'aimerais connaître les solutions dans l'intervalle de temps (en secondes) $[0; 10]$ par pas de $\Delta t = 0.1$ secondes.
Pour ce faire, j'initie
\begin{align*}
t = 0 \\
\theta = \theta_0 \\
\dot{\theta} = 0
\end{align*}
A présent, j'applique la méthode de Runge-Kutta d'ordre $4$ en respectant bien l'ordre. J'obtiens donc ceci:
\begin{align*}
k_1 &= \Delta t \ f(t, \theta, \dot{\theta}) \\
l_1 &= \Delta t \ g(t, \theta, \dot{\theta}) \\
\\
k_2 &= \Delta t f(t + \frac{\Delta t}{2}, \theta + \frac{k_1}{2}, \dot{\theta} + \frac{l1}{2}) \\
l_2 &= \Delta t \ g(t + \frac{\Delta t}{2}, \theta + \frac{k_1}{2}, \dot{\theta} + \frac{l1}{2}) \\
\\
k_3 &= \Delta t \ f \left(t + \frac{\Delta t}{2}, x + \frac{k_2}{2}, v + \frac{l_2}{2}\right) \\
l_3 &= \Delta t \ g \left(t + \frac{\Delta t}{2}, x + \frac{k_2}{2}, v + \frac{l_2}{2} \right) \\
\\
k_4 &= \Delta t \ f(t + \Delta t, \theta + k_3, \dot{\theta} + l_3 ) \\
l_4 &= \Delta t \ g(t + \Delta t, \theta + k_3, \dot{\theta} + l_3 )
\end{align*}
Une fois ces approximations intermédiaires calculées, je peux injecter ces valeurs dans les anciennes valeurs de $t$, $\theta$ et $\dot{\theta}$. Il vient que,
\begin{align*}
\theta &= \theta + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} \\
\dot{\theta} &= \dot{\theta} + \frac{l_1}{6} + \frac{l_2}{3} + \frac{l_3}{3} + \frac{l_4}{6} \\
t &= t + \Delta t
\end{align*}
J'ai pu coder en python l'algorithme précédent afin de représenter le mouvement du pendule circulaire (et de le comparer avec l'équation dans l'approximation des petites oscillations calculée analytiquement).
Maintenant, j'aimerais pouvoir résoudre un système de $2$ équations différentielles d'ordre $2$ éventuellement couplées afin d'étudier le mouvement du pendule inversé attaché à un chariot.
Je n'ai malheureusement pas la moindre idée de la façon dont je pourrais entamer le problème. J'ai pensé que l'algorithme devait fortement ressembler au précédent. Voici ma façon de procéder dans le cas très simple de la chute libre d'un corps.
$$\begin{cases} \frac{d^2 x}{dt^2} = 0, \quad \textrm{avec, $x(0) = 0$ et $\dot{x}(0) = 1$ } \\ \frac{d^2 y}{dt^2} = - g, \quad \textrm{avec, $y(0) = 0$ et $\dot{y}(0) = 1$ } \end{cases}$$
En écrivant ces $2$ équations d'ordre $2$ en $4$ équations d'ordre $1$, j'obtiens,
$$\begin{cases} \frac{d x}{dt} = f_1(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{dy}{dt} = g_1(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{d f_1}{dt} = f_2(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{d g_1}{dt} = g_2(t, x, \dot{x}, y, \dot{y}) \end{cases}$$
Il ne me reste à présent, plus qu'à implémenter l'algorithme. Voici mon code:
borne_inf (dans le problème, borne_inf = 0) et borne_sup (ici borne_sup = 1) correspondent à l'intervalle sur lequel j'aimerais résoudre le système d'équations, les paramètres suivants sont mes conditions initiales, x0_1 = $x(0) = 0$, v0_1 = $\dot{x}(0) = 1$, x0_2 = $y(0) = 0$, v0_2 = $\dot{y}(0) = 1$ et Delta_t représente la largeur des intervalles découpés dans l'intervalle [borne_inf; borne_sup].
Je remarque que les valeurs de la liste_x sont toutes égales à 0. Lors du calcul avec le pendule inversé sur le chariot, j'obtiens un plot qui ne ressemble pas vraiment à la physique du problème.
Merci de votre aide.
A titre personnel, je souhaiterais implémenter l'algorithme de Runge-Kutta d'ordre $4$ afin de résoudre un système d'équations différentielles couplées.
Pour arriver à cela, je commence par résoudre le problème du pendule circulaire (équation différentielle ordinaire d'ordre $2$):
$$\frac{d^2 \theta}{dt^2} = -g \sin(\theta)$$
Afin de résoudre numériquement ce problème, je pose
$$\begin{cases} \frac{d \theta}{dt} = f(t, \theta, \dot{\theta}) \\ \frac{df}{dt} = \frac{- g}{l} \sin(\theta) = g(t, \theta, \dot{\theta}) \end{cases}$$
Les conditions initiales sont: $$\theta(0) = \theta_0 \textrm{ et } \dot{\theta}_0 = 0$$
J'aimerais connaître les solutions dans l'intervalle de temps (en secondes) $[0; 10]$ par pas de $\Delta t = 0.1$ secondes.
Pour ce faire, j'initie
\begin{align*}
t = 0 \\
\theta = \theta_0 \\
\dot{\theta} = 0
\end{align*}
A présent, j'applique la méthode de Runge-Kutta d'ordre $4$ en respectant bien l'ordre. J'obtiens donc ceci:
\begin{align*}
k_1 &= \Delta t \ f(t, \theta, \dot{\theta}) \\
l_1 &= \Delta t \ g(t, \theta, \dot{\theta}) \\
\\
k_2 &= \Delta t f(t + \frac{\Delta t}{2}, \theta + \frac{k_1}{2}, \dot{\theta} + \frac{l1}{2}) \\
l_2 &= \Delta t \ g(t + \frac{\Delta t}{2}, \theta + \frac{k_1}{2}, \dot{\theta} + \frac{l1}{2}) \\
\\
k_3 &= \Delta t \ f \left(t + \frac{\Delta t}{2}, x + \frac{k_2}{2}, v + \frac{l_2}{2}\right) \\
l_3 &= \Delta t \ g \left(t + \frac{\Delta t}{2}, x + \frac{k_2}{2}, v + \frac{l_2}{2} \right) \\
\\
k_4 &= \Delta t \ f(t + \Delta t, \theta + k_3, \dot{\theta} + l_3 ) \\
l_4 &= \Delta t \ g(t + \Delta t, \theta + k_3, \dot{\theta} + l_3 )
\end{align*}
Une fois ces approximations intermédiaires calculées, je peux injecter ces valeurs dans les anciennes valeurs de $t$, $\theta$ et $\dot{\theta}$. Il vient que,
\begin{align*}
\theta &= \theta + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} \\
\dot{\theta} &= \dot{\theta} + \frac{l_1}{6} + \frac{l_2}{3} + \frac{l_3}{3} + \frac{l_4}{6} \\
t &= t + \Delta t
\end{align*}
J'ai pu coder en python l'algorithme précédent afin de représenter le mouvement du pendule circulaire (et de le comparer avec l'équation dans l'approximation des petites oscillations calculée analytiquement).
Maintenant, j'aimerais pouvoir résoudre un système de $2$ équations différentielles d'ordre $2$ éventuellement couplées afin d'étudier le mouvement du pendule inversé attaché à un chariot.
Je n'ai malheureusement pas la moindre idée de la façon dont je pourrais entamer le problème. J'ai pensé que l'algorithme devait fortement ressembler au précédent. Voici ma façon de procéder dans le cas très simple de la chute libre d'un corps.
$$\begin{cases} \frac{d^2 x}{dt^2} = 0, \quad \textrm{avec, $x(0) = 0$ et $\dot{x}(0) = 1$ } \\ \frac{d^2 y}{dt^2} = - g, \quad \textrm{avec, $y(0) = 0$ et $\dot{y}(0) = 1$ } \end{cases}$$
En écrivant ces $2$ équations d'ordre $2$ en $4$ équations d'ordre $1$, j'obtiens,
$$\begin{cases} \frac{d x}{dt} = f_1(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{dy}{dt} = g_1(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{d f_1}{dt} = f_2(t, x, \dot{x}, y, \dot{y}) \\[3pt] \frac{d g_1}{dt} = g_2(t, x, \dot{x}, y, \dot{y}) \end{cases}$$
Il ne me reste à présent, plus qu'à implémenter l'algorithme. Voici mon code:
borne_inf (dans le problème, borne_inf = 0) et borne_sup (ici borne_sup = 1) correspondent à l'intervalle sur lequel j'aimerais résoudre le système d'équations, les paramètres suivants sont mes conditions initiales, x0_1 = $x(0) = 0$, v0_1 = $\dot{x}(0) = 1$, x0_2 = $y(0) = 0$, v0_2 = $\dot{y}(0) = 1$ et Delta_t représente la largeur des intervalles découpés dans l'intervalle [borne_inf; borne_sup].
def RK_4(f1, g1, f2, g2, borne_inf, borne_sup, x0_1, v0_1, x0_2, v0_2, Delta_t): liste_t = [] liste_x_1 = [] liste_v_1 = [] liste_x_2 = [] liste_v_2 = [] t = borne_inf x_1 = x0_1 v_1 = v0_1 x_2 = x0_2 v_2 = v0_2 while t < borne_sup: liste_t.append(t) liste_x_1.append(x_1) liste_v_1.append(v_1) liste_x_2.append(x_2) liste_v_2.append(v_2) k1 = Delta_t * f1(t, x_1, v_1, x_2, v_2) l1 = Delta_t * g1(t, x_1, v_1, x_2, v_2) m1 = Delta_t * f2(t, x_1, v_1, x_2, v_2) n1 = Delta_t * g2(t, x_1, v_1, x_2, v_2) k2 = Delta_t * f1(t + Delta_t/2, x_1 + k1/2, v_1 + l1/2, x_2 + m1/2, v_2 + n1/2) l2 = Delta_t * g1(t + Delta_t/2, x_1 + k1/2, v_1 + l1/2, x_2 + m1/2, v_2 + n1/2) m2 = Delta_t * f2(t + Delta_t/2, x_1 + k1/2, v_1 + l1/2, x_2 + m1/2, v_2 + n1/2) n2 = Delta_t * g2(t + Delta_t/2, x_1 + k1/2, v_1 + l1/2, x_2 + m1/2, v_2 + n1/2) k3 = Delta_t * f1(t + Delta_t/2, x_1 + k2/2, v_2 + l2/2, x_2 + m2/2, v_2 + n2/2) l3 = Delta_t * g1(t + Delta_t/2, x_1 + k2/2, v_2 + l2/2, x_2 + m2/2, v_2 + n2/2) m3 = Delta_t * f2(t + Delta_t/2, x_1 + k2/2, v_2 + l2/2, x_2 + m2/2, v_2 + n2/2) n3 = Delta_t * g2(t + Delta_t/2, x_1 + k2/2, v_2 + l2/2, x_2 + m2/2, v_2 + n2/2) k4 = Delta_t * f1(t + Delta_t, x_1 + k3, v_1 + l3, x_2 + m3, v_2 + n3) l4 = Delta_t * g1(t + Delta_t, x_1 + k3, v_1 + l3, x_2 + m3, v_2 + n3) m4 = Delta_t * f2(t + Delta_t, x_1 + k3, v_1 + l3, x_2 + m3, v_2 + n3) n4 = Delta_t * g2(t + Delta_t, x_1 + k3, v_1 + l3, x_2 + m3, v_2 + n3) x_1 += k1/6 + k2/3 + k3/3 + k4/6 v_1 += l1/6 + l2/3 + l3/3 + l4/6 x_2 += m1/6 + m2/3 + m3/3 + m4/6 v_2 += n1/6 + n2/3 + n3/3 + n4/6 t += Delta_t return liste_t, liste_x_1, liste_v_1, liste_x_2, liste_v_2 def f1(t, x, deriv_x, y, deriv_y): return deriv_x def g1(t, x, deriv_x, y, deriv_y): return deriv_y def f2(t, x, deriv_x, y, deriv_y): return 0 def g2(t, x, deriv_x, y, deriv_y): return - 9.81 liste_t, liste_x, liste_deriv_x, liste_y, liste_deriv_y = RK_4(f1, g1, f2, g2, 0, 5, 0, 1, 0, 1, 0.05) print(liste_x)
Je remarque que les valeurs de la liste_x sont toutes égales à 0. Lors du calcul avec le pendule inversé sur le chariot, j'obtiens un plot qui ne ressemble pas vraiment à la physique du problème.
Merci de votre aide.
Réponses
-
Chez moi, c'est liste_y qui est pleine de zéros. Vu que f2 renvoie toujours zéro (au lieu de $-9{,}81$ pour g2) et que f2 est associé à m1,...,m4 qui contribuent à x_2, ce n'est pas si étonnant.
Dans les définitions, k, l, m, n correspondent dans cet ordre à f1, g1, f2, g2, qui devraient correspondre à x1, x2, v1, v2 alors que tu les ajoutes à x1, v1, x2, v2 (enlève ces _ par homogénéité !). Ça pourrait être un problème. -
Effectivement, il s'agit bien de liste_y qui ne contient que des $0$ (la modification a été effectuée dans le message principal).
En enlevant les "_" et en tenant compte de vos remarques, j'associe les bons coefficients $k$, $l$, $m$ et $n$ aux bonnes valeurs $x1$, $x2$, $v1 $et $v2$
La toute fin de ma fonction devient alors,x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat
Je remarque que le mouvement obtenu ne correspond pas vraiment au résultat recherché (cnfr. image jointe à ce message).def RK4(f1, g1, f2, g2, borneInf, borneSup, x01, v01, x02, y02, Deltat): listet = [] listex1 = [] listex2 = [] listev1 = [] listev2 = [] t = borneInf x1 = x01 v1 = v01 x2 = x02 v2 = x02 while t < borneSup: listet.append(t) listex1.append(x1) listex2.append(x2) listev1.append(v1) listev2.append(v2) k1 = Deltat * f1(t, x1, v1, x2, v2) l1 = Deltat * g1(t, x1, v1, x2, v2) m1 = Deltat * f2(t, x1, v1, x2, v2) n1 = Deltat * g2(t, x1, v1, x2, v2) k2 = Deltat * f1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) l2 = Deltat * g1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) m2 = Deltat * f2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) n2 = Deltat * g2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) k3 = Deltat * f1(t + Deltat/2, x1 + k2/2, v2 + l2/2, x2 + m2/2, v2 + n2/2) l3 = Deltat * g1(t + Deltat/2, x1 + k2/2, v2 + l2/2, x2 + m2/2, v2 + n2/2) m3 = Deltat * f2(t + Deltat/2, x1 + k2/2, v2 + l2/2, x2 + m2/2, v2 + n2/2) n3 = Deltat * g2(t + Deltat/2, x1 + k2/2, v2 + l2/2, x2 + m2/2, v2 + n2/2) k4 = Deltat * f1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) l4 = Deltat * g1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) m4 = Deltat * f2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) n4 = Deltat * g2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat return listet, listex1, listev1, listex2, listev2 def f1(t, x, derivx, y, derivy): return derivx def g1(t, x, derivx, y, derivy): return derivy def f2(t, x, derivx, y, derivy): return 0 def g2(t, x, derivx, y, derivy): return - 9.81 listet, listex, listederivx, listey, listederivy = RK4(f1, g1, f2, g2, 0, 0.2, 0, 1, 0, 1, 0.01) def representation_graphique(t_i, x_i): plt.scatter(t_i, x_i) plt.xlabel('t') plt.legend() representation_graphique(listet, listex) representation_graphique(listet, listey) plt.show() -
Il y a d'autres incohérences du même genre, dues principalement à un ordre fluctuant sur les éléments centraux de la liste (x1,v1,x2,v2) :
- dans le calcul de k2 etc., tu fais intervenir v1+l1/2 et x2+m1/2. Or dans la somme de la fin de la boucle, les l* correspondent à x2 et les m* correspondent à v1, ce qui n'est pas cohérent ;
- dans le calcul de k3 etc., c'est pire : il apparaît v2+l2/2 au lieu de v1+l2/2 pour faire la même erreur qu'au point précédent ;
- dans le calcul de k4 etc., il y a la même incohérence qu'au premier point.
- dans la définition de RK, on croit que g1 correspond à la vitesse du point 1 et f2 à la position du point 2 mais dans la définition de ces fonctions qui suit, c'est l'inverse.
- dans les lignes k3, l3, m3, n3 : remplacer les v2 « entre x1 et x2 » par v1 ;
- dans la définition de x1, v1, x2, v2 : permuter v1 et x2 ;
- permuter les définitions de f2 et g1.
- y02 au lieu de v02
- "v2 = v02" au lieu de "v2 = x02"
Après, comment vas-tu faire quand tu auras trois points dans l'espace ? Vas-tu écrire une procédure par problème physique ? Il vaudrait mieux, je crois, écrire une bonne fois les formules de Runge-Kutta avec une seule inconnue vectorielle (ou de type list ou tuple, si tu n'as pas de vecteurs). Tu aurais ainsi une seule ligne pour chaque calcul intermédiaire de k1, k2, k3, k4, ce serait quand même plus facile d'être cohérent. -
Modifications apportées:
- Permutation des paramètres de ma fonction RK4: x01, v01, x02, v02 devient x01, x02, v01, v02
- Modification des coefficients $k_3$, $l_3$, $m_3$, $n_3$. Les v_2 ont été remplacés par les v_1
- Corrections des fautes de frappe
- Permutation de $f_2$ et $g_1$ dans la fonction et à l'appel de la fonction (pour rester cohérent avec ce qui a été écrit)
J'obtiens enfin le résultat voulu. Grâce à votre aide, j'ai pu corriger mon code et j'ai pu comprendre les erreurs que j'avais faites. Je vous en remercie. Cela me sera très utile !
Je ne code en python que depuis peu de temps.
De plus, l'algorithme de résolution d'EDOs de Runge-Kutta est un algorithme assez compliqué à comprendre au premier abord. J'ai d'abord écrit une fonction RK4 pour résoudre une équation différentielle du premier ordre.
Ensuite, j'en ai écrit une qui résout des EDOs du second ordre et me voilà enfin à écrire une fonction de résolution d'un système de deux EDOs du second ordre.
Cette dernière procédure m'a donné beaucoup de fil à retordre. Je vais donc envisager l'écriture d'une procédure avec une seule inconnue vectorielle.
Code modifié:def RK4(f1, f2, g1, g2, borneInf, borneSup, x01, x02, v01, v02, Deltat): listet = [] listex1 = [] listex2 = [] listev1 = [] listev2 = [] t = borneInf x1 = x01 v1 = v01 x2 = x02 v2 = v02 while t < borneSup: listet.append(t) listex1.append(x1) listex2.append(x2) listev1.append(v1) listev2.append(v2) k1 = Deltat * f1(t, x1, v1, x2, v2) l1 = Deltat * g1(t, x1, v1, x2, v2) m1 = Deltat * f2(t, x1, v1, x2, v2) n1 = Deltat * g2(t, x1, v1, x2, v2) k2 = Deltat * f1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) l2 = Deltat * g1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) m2 = Deltat * f2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) n2 = Deltat * g2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) k3 = Deltat * f1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) l3 = Deltat * g1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) m3 = Deltat * f2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) n3 = Deltat * g2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) k4 = Deltat * f1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) l4 = Deltat * g1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) m4 = Deltat * f2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) n4 = Deltat * g2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat return listet, listex1, listev1, listex2, listev2 def f1(t, x, derivx, y, derivy): return derivx def g1(t, x, derivx, y, derivy): return derivy def f2(t, x, derivx, y, derivy): return 0 def g2(t, x, derivx, y, derivy): return - 9.81 listet, listex, listederivx, listey, listederivy = RK4(f1, f2, g1, g2, 0, 0.220, 0, 0, 1, 1, 0.01) plt.scatter(listet, listex, label = "x(t)") plt.scatter(listet, listey, label = "y(t)") plt.xlabel('t') plt.legend() plt.show() -
Cette fois, on y croit !
Et si on mettait un ressort entre ces deux points ? -
Le problème étudié ici est la chute libre d'un seul corps dans le champ de la pesanteur (en l'absence des frottements de l'air). De quels points parlez-vous ?
-
Oui. En regardant le système un peu différemment, on peut aussi l'interpréter comme deux points matériels qui se déplacent sur une droite verticale, un de masse nulle et un autre de masse $1$, que l'on envoie initialement avec une vitesse verticale $1$ (dirigée vers le haut). Dans ce sens, on peut ajouter un ressort entre ces deux points.
À tort ou à raison, je me suis mis à voir la situation comme ça à cause des notations x1 et x2 au lieu de x et y (alors que c'était introduit comme « première » et « deuxième coordonnées ») et du dessin des deux coordonnées séparément en fonction du temps, au lieu d'une « courbe paramétrée » $(x(t),y(t))$. -
Si je saisis bien, on demande d'approximer numériquement l'évolution de la position du système constitué d'une particule $y_1$ de masse $m_1$ connectée au moyen d'un ressort de constante $k$ à une autre particule $y_2$ de masse $m_2$.
Les deux particules sont initialement en $y = 0$ et on lance $y_1$ avec une vitesse initiale d'intensité $v_0 = 1$.
Le système d'équations différentielles couplées qui régit le mouvement des particules est alors donné par,
\begin{align*}
&\begin{cases} m \frac{d^2 y_1}{dt^2} = -ky_1 + ky_2 - mg \\[3pt] m \frac{d^2 y_2}{dt} = ky_1 - ky_2 -mg \end{cases} \\[5pt]
\iff & \begin{cases} \frac{d^2 y_1}{dt^2} = -\frac{k}{m}y_1 + \frac{k}{m}y_2 - g \\[3pt] \frac{d^2 y_2}{dt} = \frac{k}{m}y_1 - \frac{k}{m}y_2 -g \end{cases} \\[5pt]
\iff & \begin{cases} \frac{d y_1}{dt} = f_1(t, y_1, \dot{y_1}, y_2, \dot{y_2})\\[3pt] \frac{df_1}{dt} = f_2 (t, y_1, \dot{y_1}, y_2, \dot{y_2}) = -\frac{k}{m}y_1 + \frac{k}{m}y_2 - g \\[4pt] \frac{d y_2}{dt} = g_1(t, y_1, \dot{y_1}, y_2, \dot{y_2})\\[3pt] \frac{dg_1}{dt} = g_2 (t, y_1, \dot{y_1}, y_2, \dot{y_2}) = \frac{k}{m}y_1 - \frac{k}{m}y_2 - g \end{cases}
\end{align*}
avec,
$y_1(0) = 0$
$y_1(0) = 0$
$\dot{y}_1(0) = 1$
$\dot{y}_2(0) = 0$ -
Oui, quelque chose comme ça. Mon idée, c'était de profiter du code en place pour résoudre un problème différent avec un code presque identique (il suffit de changer f2 et g2 en ajoutant un terme $\pm K(y_1-y_2)$, éventuellement avec des constantes différentes pour indiquer des masses différentes). Ce n'est pas palpitant palpitant non plus, hein...
-
Je n'ai pas du modifier grand chose à mon code initial. Le résultat à l'air de bien rendre compte de la physique du problème.
Note: j'ai du changer un peu les conditions initiales pour pouvoir avoir deux plots bien distincts. En effet, lorsque la vitesse de la première particule est de $1 \ ms^{-1}$, les mouvements sont ratatinés.from math import * from pylab import * def RK4(f1, f2, g1, g2, borneInf, borneSup, x01, x02, v01, v02, Deltat): listet = [] listex1 = [] listex2 = [] listev1 = [] listev2 = [] t = borneInf x1 = x01 v1 = v01 x2 = x02 v2 = v02 while t < borneSup: listet.append(t) listex1.append(x1) listex2.append(x2) listev1.append(v1) listev2.append(v2) k1 = Deltat * f1(t, x1, v1, x2, v2) l1 = Deltat * g1(t, x1, v1, x2, v2) m1 = Deltat * f2(t, x1, v1, x2, v2) n1 = Deltat * g2(t, x1, v1, x2, v2) k2 = Deltat * f1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) l2 = Deltat * g1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) m2 = Deltat * f2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) n2 = Deltat * g2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) k3 = Deltat * f1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) l3 = Deltat * g1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) m3 = Deltat * f2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) n3 = Deltat * g2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) k4 = Deltat * f1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) l4 = Deltat * g1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) m4 = Deltat * f2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) n4 = Deltat * g2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat return listet, listex1, listev1, listex2, listev2 k = 100 m = 1 def f1(t, y1, derivy1, y2, derivy2): return derivy1 def f2(t, y1, derivy1, y2, derivy2): return - (k/m)*y1 + (k/m)*y2 - 9.81 def g1(t, y1, derivy1, y2, derivy2): return derivy2 def g2(t, y1, derivy1, y2, derivy2): return (k/m)*y1 - (k/m)*y2 - 9.81 y01 = 0 y02 = 0 v01 = 10 v02 = 0 listet = [] listey1 = [] listey2 = [] listev1 = [] listev2 = [] listet, listey1, listev1, listey2, listev2 = RK4(f1, f2, g1, g2, 0, 1, y01, y02, v01, v02, 0.01) plt.scatter(listet, listey1, label = "y1(t)") plt.scatter(listet, listey2, label = "y2(t)") plt.xlabel('t') plt.legend() plt.show() -
Ça swingue !
Et maintenant, un codage qui prend une variable vectorielle en argument, pour résoudre tous les problème de physique à la fois ? -
Bien sûr ! Mais avant cela, j'aimerais discuter du problème du pendule inversé attaché à un mobile (cnfr. image).
Les équations qui régissent le mouvement du mobile sont données par,
\begin{align*}
&\begin{cases}
\frac{d^2 x}{dt^2} = \frac{1}{M - m \cos^2(\theta)} \left(F + mg \cos(\theta)\sin(\theta) - c \dot{x} - m l \sin(\theta) \dot{\theta}^2 \right) \\[10pt] \frac{d^2 \theta}{dt^2} = \frac{1}{l} \left\{g \sin(\theta) + (F + mg \cos(\theta) \sin(\theta) - c \dot{x} - m l \sin(\theta) \dot{\theta}^2) \right\} \frac{\cos(\theta)}{M - m \cos^2(\theta)}
\end{cases} \\[20 pt]
\iff & \begin{cases}
\frac{d x}{dt} = f_1(t, x, \dot{x}, \theta, \dot{\theta}) \\[5pt] \frac{d f_1}{dt} = f_2(t, x, \dot{x}, \theta, \dot{\theta}) = \frac{1}{M - m \cos^2(\theta)} \left(F + mg \cos(\theta)\sin(\theta) - c \dot{x} - m l \sin(\theta) \dot{\theta}^2 \right) \\[5pt] \frac{d \theta}{dt} = g_1(t, x, \dot{x}, \theta, \dot{\theta}) \\[5pt] \frac{d g_1}{dt} = g_2(t, x, \dot{x}, \theta, \dot{\theta}) = \frac{1}{l} \left\{g \sin(\theta) + (F + mg \cos(\theta) \sin(\theta) - c \dot{x} - m l \sin(\theta) \dot{\theta}^2) \right\}
\end{cases} \\
\end{align*}
$l = 0.5$ longueur du fil du pendule
g = 9.81
$M = 2$ masse du mobile
$m = 0.1$ masse de la particule attachée au pendule
$c = 1$ constante associée à la force de frottement
$F = 2$ intensité de $F$
Je souhaiterais étudier le mouvement de ce système aux conditions initiales $x(0) = \dot{x}(0) = \theta(0) = \dot{\theta}(0) = 0$.
Compte-tenu de $\theta$, il est clair que tant que $F = 0$, le système est en équilibre. Mais dans mon cas, j'applique une perturbation $F = 2$ et je représente l'évolution de $\theta$ et de $x$ au cours du temps.
Je suis peu satisfait de la courbe obtenue car étant donné la force de départ, le mobile devrait se déplacer avant "d'osciller" vers une nouvelle position d'équilibre. De plus, les oscillations ont l'air de diverger lorsque $t \to \infty$from math import * from pylab import * def RK4(f1, f2, g1, g2, borneInf, borneSup, x01, x02, v01, v02, Deltat): listet = [] listex1 = [] listex2 = [] listev1 = [] listev2 = [] t = borneInf x1 = x01 v1 = v01 x2 = x02 v2 = v02 while t < borneSup: listet.append(t) listex1.append(x1) listex2.append(x2) listev1.append(v1) listev2.append(v2) k1 = Deltat * f1(t, x1, v1, x2, v2) l1 = Deltat * g1(t, x1, v1, x2, v2) m1 = Deltat * f2(t, x1, v1, x2, v2) n1 = Deltat * g2(t, x1, v1, x2, v2) k2 = Deltat * f1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) l2 = Deltat * g1(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) m2 = Deltat * f2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) n2 = Deltat * g2(t + Deltat/2, x1 + k1/2, v1 + l1/2, x2 + m1/2, v2 + n1/2) k3 = Deltat * f1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) l3 = Deltat * g1(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) m3 = Deltat * f2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) n3 = Deltat * g2(t + Deltat/2, x1 + k2/2, v1 + l2/2, x2 + m2/2, v2 + n2/2) k4 = Deltat * f1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) l4 = Deltat * g1(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) m4 = Deltat * f2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) n4 = Deltat * g2(t + Deltat, x1 + k3, v1 + l3, x2 + m3, v2 + n3) x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat return listet, listex1, listev1, listex2, listev2 l = 0.5 g = 9.81 M = 2 m = 0.1 c = 1 F = 5 def f1(t, x, derivx, theta, derivtheta): if t <= 0.010: F = 2 else: F = 0 return derivx def f2(t, x, derivx, theta, derivtheta): if t <= 0.010: F = 2 else: F = 0 return (1/(M - (m*(cos(theta))*(cos(theta))))) * (F + (m*g*(cos(theta))*(sin(theta))) - (c*derivx) - (m*l*derivtheta*derivtheta*(sin(theta)))) def g1(t, x, derivx, theta, derivtheta): if t <= 0.010: F = 2 else: F = 0 return derivtheta def g2(t, x, derivx, theta, derivtheta): if t <= 0.010: F = 2 else: F = 0 return (1/l) * ((g*(sin(theta))) + (F + (m*g*(cos(theta))*(sin(theta))) - (c*derivx) - (m*l*derivtheta*derivtheta*(sin(theta))))) * ((cos(theta))/(M - (m*(cos(theta))*(cos(theta))))) theta0 = 0 theta0rad = np.radians(theta0) derivtheta0 = 0 derivtheta0rad = np.radians(derivtheta0) x0 = 0 derivx0 = 0 listet, listex, listederivx, listetheta, listederivtheta = RK4(f1, f2, g1, g2, 0, 5, x0, theta0, derivx0, derivtheta0rad, 0.01) plt.scatter(listet, listex, label = "x(t)") plt.scatter(listet, listetheta, label = "thêta(t)") plt.xlabel('t') plt.legend() plt.show() -
La situation physique m'échappe : pourquoi est-ce que ça se mettrait à osciller ? Si on part de l'équilibre et qu'on tire le chariot vers la droite pendant 1/10 s, le pendule ne va-t-il pas se mettre à tomber vers la gaucheet faire des tours complets ?
À part ça, tu es reparti dans les mêmes travers que plus haut : tu calcules k2, k3, etc. en faisant intervenir "v1+l*/*" et tu poses à la fin "x2 += l1/6 + l2/3 + l3/3 + l4/6". En permutant x2 et v1 dans la définition, voici la figure que j'obtiens (l'abscisse est graduée en nombre de deltat et pas en secondes). Comme je ne comprends pas la situation, je ne sais pas dire si c'est convaincant.
PS : animation – crédible ? -
Effectivement, il y avait encore une coquille au niveau de la fonction. J'ai switché $x2$ et $v1$.
Enfin,x1 += k1/6 + k2/3 + k3/3 + k4/6 v1 += l1/6 + l2/3 + l3/3 + l4/6 x2 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat
EDIT 13h00: les $k_i$ sont associées à $f_1$, les $l_i$ à $g_1$. Donc, je dois switcher les $m$ et les $l$ lors de l'appel des fonctions $f_i$ et $g_i$ pour déterminer les coefficients. Voici mon nouveau codedef RK4(f1, f2, g1, g2, borneInf, borneSup, x01, x02, v01, v02, Deltat): listet = [] listex1 = [] listex2 = [] listev1 = [] listev2 = [] t = borneInf x1 = x01 v1 = v01 x2 = x02 v2 = v02 while t < borneSup: listet.append(t) listex1.append(x1) listex2.append(x2) listev1.append(v1) listev2.append(v2) k1 = Deltat * f1(t, x1, v1, x2, v2) l1 = Deltat * g1(t, x1, v1, x2, v2) m1 = Deltat * f2(t, x1, v1, x2, v2) n1 = Deltat * g2(t, x1, v1, x2, v2) k2 = Deltat * f1(t + Deltat/2, x1 + k1/2, v1 + m1/2, x2 + l1/2, v2 + n1/2) l2 = Deltat * g1(t + Deltat/2, x1 + k1/2, v1 + m1/2, x2 + l1/2, v2 + n1/2) m2 = Deltat * f2(t + Deltat/2, x1 + k1/2, v1 + m1/2, x2 + l1/2, v2 + n1/2) n2 = Deltat * g2(t + Deltat/2, x1 + k1/2, v1 + m1/2, x2 + l1/2, v2 + n1/2) k3 = Deltat * f1(t + Deltat/2, x1 + k2/2, v1 + m2/2, x2 + l2/2, v2 + n2/2) l3 = Deltat * g1(t + Deltat/2, x1 + k2/2, v1 + m2/2, x2 + l2/2, v2 + n2/2) m3 = Deltat * f2(t + Deltat/2, x1 + k2/2, v1 + m2/2, x2 + l2/2, v2 + n2/2) n3 = Deltat * g2(t + Deltat/2, x1 + k2/2, v1 + m2/2, x2 + l2/2, v2 + n2/2) k4 = Deltat * f1(t + Deltat, x1 + k3, v1 + m3, x2 + l3, v2 + n3) l4 = Deltat * g1(t + Deltat, x1 + k3, v1 + m3, x2 + l3, v2 + n3) m4 = Deltat * f2(t + Deltat, x1 + k3, v1 + m3, x2 + l3, v2 + n3) n4 = Deltat * g2(t + Deltat, x1 + k3, v1 + m3, x2 + l3, v2 + n3) x1 += k1/6 + k2/3 + k3/3 + k4/6 x2 += l1/6 + l2/3 + l3/3 + l4/6 v1 += m1/6 + m2/3 + m3/3 + m4/6 v2 += n1/6 + n2/3 + n3/3 + n4/6 t += Deltat return listet, listex1, listev1, listex2, listev2
Dès que j'ai tiré le mobile vers la droite, le pendule devrait tomber vers la gauche et effectuer des oscillations de forte amplitude autour de la position d'équilibre $\theta = \pi$. Or ici, lorsque le mobile part vers la droite, le pendule oscille un peu en direction de la droite avant de revenir (comme si un force de rappel était appliquée sur le pendule). En plus de cela, avec les frottements entre le mobile et le sol, le mobile devrait s'arrêter à un moment donné, ce qui n'est pas le cas ici.
EDIT 14h00: Avec ma nouvelle fonction, j'obtiens un résultat qui est bien plus satisfaisant. Dès que je perturbe le système vers la droite, le pendule se met à osciller (vers les $\theta$ positifs). Le pendule fait quelques tours sur lui même avant que le mouvement du mobile n'empêche le pendule d'effectuer des tours complets à nouveau. Le pendule se met donc à osciller.
Comment arrivez-vous à créer une telle animation ? Je suis curieux.
(à 9:19, il y a une petite animation montrant le comportement d'un tel dispositif constitué d'un pendule inversé et d'un mobile. On remarque que le mouvement ne diverge pas. Or, si on fait la simulation sur notre programme, on obtient un mouvement qui "explose" lorsque $t$ devient grand) -
Bizarre quand même...
Pour faire l'animation, rien de sorcier : on dessine un chariot et un pendule, on sauve un png par instant et on convertit ça en un gif animé. Voici les instructions en Sage, elles doivent être assez transparentes pour être transcrites en Python "pur".def char(x): G = polygon([(x-1,.5),(x+1,.5),(x+1,1.5),(x-1,1.5)],fill=True,color="lightgray") G+= circle((x-.5,.25),.25,fill=True,color="lightgray") G+= circle((x+.5,.25),.25,fill=True,color="lightgray") return G def pendule(x,th,l=2): u,v = x-l*sin(th),1.5+l*cos(th) G = line([(x,1.5),(u,v)],color="black",thickness=3) G+= circle((u,v),.1,fill=True,color="black") return G def allonge(k,l=2): ch = str(k) if len(ch)<l: ch = (l-len(ch))*"0"+ch return ch N = 10 for k in range(3000//N): G = char(5*listex[N*k])+pendule(5*listex[N*k],listetheta[N*k]) G.save('/home/jerome/Bureau/cyclo/chariot/pendule_%s.png'%allonge(k,3),xmin=-5,xmax=8,ymax=3.8) -
C'est vraiment très bizarre...
Les $10$ premières secondes se passent plutôt bien (même si le pendule semble ralentir à l'approche de $\theta = \pi$) mais après ces $10$ secondes, le mouvement ne rend plus du tout compte de la réalité.
Je ne vois pas du tout d'où pourrait bien provenir les erreurs. Je suis pratiquement sûr de la fonction $RK4$, à présent et les équations d'évolution de $x$ et de $\theta$ sont correctes... -
Pour l'heure, je vais supposer que la fonction est correcte mais qu'il y a une erreur autre part (mauvaises équations, équations mal écrites en python avec toutes les parenthèses, ...)
Le graphique de $y_1$ et $y_2$ (dans le problème de la chute libre de deux particules ponctuelles attachées à un ressort) n'est pas correct car je n'avais pas encore switché les $m_i$ et $l_i$ dans ma fonction. C'est chose faite!
J'ai écrit deux nouvelles fonctions RK3 et RK5, j'ai tout représenté graphiquement avec le graphe de la solution calculée analytiquement. On peut s'apercevoir que RK3 est moins précis que RK4 qui est lui-même moins précis que RK5.
Je vais passer à la partie avec les vecteurs. Existe-t-il un objet vecteur en python ? -
Il y a les tableaux de numpy, cf. par exemple https://stackoverflow.com/questions/12049154/python-numpy-vector-math.
Sinon, ça ne coûte pas cher de le créer. Il suffit de créer une classe "vecteur" qui étend la classe "tuple" et permet la multiplication par un scalaire, surtout la somme de deux vecteurs, le produit scalaire, le produit vectoriel de deux vecteurs à trois composantes, etc.
Je n'ai pas de Python sous la main mais ça pourrait partir comme ça.class vecteur(tuple): def __init__(self,L): if type(L)=='list': self.v = tuple(L) else: self.v = L def ajoute(self,W): return vecteur([self.v[k]+W.v[k] for k in range(len(self.v))]) def multiplie(self,x): return vecteur([x*u for u in self.v]) V = vecteur([1,2,3]) print V.ajoute(V) print V+V -
Je vois. Je vous remercie beaucoup !
Connectez-vous ou Inscrivez-vous pour répondre.
Bonjour!
Catégories
- 163.2K Toutes les catégories
- 9 Collège/Lycée
- 21.9K Algèbre
- 37.1K Analyse
- 6.2K Arithmétique
- 53 Catégories et structures
- 1K Combinatoire et Graphes
- 11 Sciences des données
- 5K Concours et Examens
- 11 CultureMath
- 47 Enseignement à distance
- 2.9K Fondements et Logique
- 10.3K Géométrie
- 65 Géométrie différentielle
- 1.1K Histoire des Mathématiques
- 69 Informatique théorique
- 3.8K LaTeX
- 39K Les-mathématiques
- 3.5K Livres, articles, revues, (...)
- 2.7K Logiciels pour les mathématiques
- 24 Mathématiques et finance
- 314 Mathématiques et Physique
- 4.9K Mathématiques et Société
- 3.3K Pédagogie, enseignement, orientation
- 10K Probabilités, théorie de la mesure
- 773 Shtam
- 4.2K Statistiques
- 3.7K Topologie
- 1.4K Vie du Forum et de ses membres
Qui est en ligne 21
+11 Invités