Résolution d'équations de mouvement en Python (odeint)

Cet article

Résolvez des équations cinétiques à l'aide d'une bibliothèque appelée odeint en python (scipy pour être exact). J'utilise python3.6.3 que j'ai mis dans anaconda.

Il résout les équations cinétiques (équations différentielles du second ordre), mais il peut également être appliqué à d'autres équations différentielles telles que les équations différentielles du premier ordre.

Je veux résoudre une équation de mouvement avec Python

Qu'est ce que la vie

Il y a des moments dans la vie où vous voulez résoudre une simple équation de mouvement. En particulier, il y a des moments où vous souhaitez résoudre et tracer l'état de l'exercice. Que faire à ce moment

Cette fois, le but est de profiter du traçage en utilisant la bibliothèque python "odeint".

odeint

On dirait Oden. Je pense que c'est intégral ODE. ODE = équation différentielle ordinaire

Référence officielle https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

Il existe également une bibliothèque appelée ode https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html

Après examen, il semble qu'il existe un nouveau paquetage "résoudre_ivp", qui est recommandé. Je veux écrire un article à ce sujet un jour. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

À propos, odeint et les deux autres (ode, solve_ivp) ont un ordre d'arguments différent. Veillez à ne pas vous confondre.

L'équation de mouvement que vous souhaitez résoudre cette fois

Toute équation cinétique fera l'affaire, mais cette fois je résoudrai un simple lancer.

m\ddot y=-rv+mg

Il est nécessaire de convertir l'équation différentielle du second ordre avec odeint en deux différentiels du premier ordre en particulier. Par exemple, dans ce cas

\dot y=v \\
m\dot v=-rv+mg

Est en fait résolu par odeint.

code

Flux global

Le flux de code est le suivant

1. Chargez la bibliothèque requise
2. Décrivez l'équation différentielle que vous souhaitez résoudre
3. Décidez des conditions initiales et demandez à odeint de résoudre l'équation.
4. Dessinez le résultat

Préparation

Tout d'abord, mettez dans diverses bibliothèques.

import

import numpy as np #numpy
from scipy.integrate import odeint #odeint
import matplotlib.pyplot as plt #to draw graphs

Ecrire une équation différentielle

Définissez ensuite l'équation différentielle. La résistance de l'air r et la masse m du point de qualité sont définies à l'étape suivante. Puisque la constante de gravité g n'est pas modifiée, elle est définie ici cette fois.

function

def func(s, t, r, m):
    
    y, v = s #s est une paire de variables y et v
    g=9.80665 #m/s^2
    dsdt = [v, (-r*v-m*g)/m]
    return dsdt 

Décidez des conditions initiales et demandez-leur de résoudre l'équation différentielle

Déterminez la résistance de l'air r et la masse m, déterminez les conditions initiales et résolvez l'équation différentielle. Notez un point sur le dernier argument de sol = odeint (func, y0, t, args = (r, m)). Vous passez une liste comme args = (r, m), donc si vous avez un argument, vous devez faire quelque chose comme args = (r,). Cette fois, ce n'est pas un problème, mais soyez prudent lorsque vous résolvez d'autres équations cinétiques.

solve

r=12
m=100
y0 = [0,10]#Lancer de la position 0 à la vitesse initiale 10
t = np.linspace(0, 7, 201)#Calculer par incréments de 201 étapes entre 0 et 7

sol = odeint(func, y0, t, args=(r,m))

Dessinez le résultat dans un diagramme

Visualisons le résultat.

visualize

plt.plot(t, sol[:, 0], 'b', label='y')#à propos de y plot
plt.plot(t, sol[:, 1], 'g', label='v')#à propos de v plot
plt.legend(loc='best')#Ajouter une légende
plt.xlabel('t')
plt.grid()#Ajouter une grille
plt.show()

Ensuite, vous devriez obtenir ce chiffre.

Vient ensuite

Dans un premier temps, changer t de 0 à 100 pour vérifier la vitesse de terminaison, ou résoudre le problème en le jetant en diagonale est une bonne pratique pour comprendre le code.