Analyse numérique des équations différentielles ordinaires avec l'odeint et l'ode de Scipy

Une méthode de simulation (analyse / calcul numérique) d'un modèle physique sous la forme d'une équation différentielle normale telle qu'un simulateur de vol utilisant python. Ce que je fais, c'est utiliser odeint dans Scipy.integrate. Il semble qu'il utilise le lsode d'Odepack de Fortran, donc le calcul est rapide.

Exemple

À titre d'exemple, portons la simulation de vol de l'avion publié sous forme de code Octave (compatible Matlab) vers python. Le lien est détaillé pour le contenu du calcul. Butterfly_Effect( ) 6DOF Flight Simulation

Si vous utilisez for pour simuler quelque chose qui bouge dans le temps, même une simple méthode de calcul numérique telle que la méthode Euler est assez lente en python, vous devez donc utiliser autant que possible une bibliothèque telle que scipy.

Différence entre matlab et python

Notez que lors du portage depuis matlab, l'opérateur * de la matrice (ndarray) est différent de matlab: inner product et numpy: element product. Notez également qu'en python, si le point décimal n'est pas spécifié, une division entière se produira.

Code --odeint

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import odeint
from scipy.linalg import block_diag
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.close('all')

# import time
# start = time.time()

def Rotation_X(Psi):
	R_x = np.array([[np.cos(Psi), np.sin(Psi), 0],
					[-np.sin(Psi), np.cos(Psi), 0],
					[0, 0, 1]])
	return R_x

def Rotation_Y(Theta):
	R_y = np.array([[np.cos(Theta), 0, -np.sin(Theta)],
					[0, 1, 0],
					[np.sin(Theta), 0, np.cos(Theta)]])
	return R_y

def Rotation_Z(Phi):
	R_z = np.array([[1, 0, 0],
					[0, np.cos(Phi), np.sin(Phi)],
					[0, -np.sin(Phi), np.cos(Phi)]])
	return R_z

#Équation du mouvement
def dynamical_system(x, t, A, U0):
	# x = [u,alpha,q,theta,beta,p,r,phi,psi,x,y,z]
	dx = A.dot(x)
	u = x[0]+U0 #la vitesse
	UVW = np.array([u, u*x[4], u*x[1]]) #Vecteur de vitesse[U,V,W]
	Rotation = Rotation_X(-x[8]).dot(Rotation_Y(-x[3])).dot(Rotation_Z(-x[7]))
	dX = Rotation.dot(UVW)
	dx[9]  = dX[0]
	dx[10] = dX[1]
	dx[11] = dX[2] 
	return dx

#Coefficient fin de stabilité dimensionnelle
Xu = -0.01;   Zu = -0.1;   Mu = 0.001
Xa = 30.0;    Za = -200.0; Ma = -4.0
Xq = 0.3;     Zq = -5.0;   Mq = -1.0
Yb = -45.0;   Lb_= -2.0;   Nb_= 1.0
Yp = 0.5;     Lp_= -1.0;   Np_= -0.1
Yr = 3.0;     Lr_= 0.2;    Nr_=-0.2

#Autres paramètres
W0 = 0.0;     U0 = 100.0;  theta0 = 0.05
g  = 9.8 #Accélération gravitationnelle

#Système vertical
A_lat = np.array([[Xu, Xa, -W0, -g*np.cos(theta0)],
    			  [Zu/U0, Za/U0, (U0+Zq)/U0, -g*np.sin(theta0)/U0],
    			  [Mu, Ma, Mq, 0],
    			  [0, 0, 1, 0]])

#Système latéral / directionnel
A_lon = np.array([[Yb, (W0+Yp), -(U0-Yr), g*np.cos(theta0), 0],
        		  [Lb_, Lp_, Lr_, 0, 0],
        		  [Nb_, Np_, Nr_, 0, 0],
        		  [0, 1, np.tan(theta0), 0, 0],
        		  [0, 0, 1/np.cos(theta0), 0, 0]])

#Combinez les systèmes sous forme de blocs diagonaux
A = block_diag(A_lat, A_lon)

#De plus, un espace sécurisé pour la trajectoire de vol
A = block_diag(A, np.zeros([3,3]))

#Définition des conditions de calcul
endurance	= 200	#Temps de vol[sec]
step		= 10	# 1.0[sec]Nombre de pas de temps par
t = np.linspace(0,endurance,endurance*step)

#Valeur initiale x0= [u,alpha,q,theta, beta,p,r,phi,psi]
#Les variables à inclure dans l'équation différentielle normale doivent être unidimensionnelles.
x0_lat		= np.array([10, 0.1, 0.4, 0.2])		#Valeur initiale verticale
x0_lon		= np.array([0.0, 0.6, 0.4, 0.2, 0.2])	#Valeur initiale horizontale / direction
x0_pos		= np.array([0, 0, -1000]) #Position initiale de l'avion
x0 = np.hstack((x0_lat, x0_lon, x0_pos))

x = odeint(dynamical_system, x0, t, args=(A,U0,))
print u"run successfully."

#Visualisation
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)

ax.plot(x[:,9], x[:,10], x[:,11])
ax.set_xlabel(u"x");ax.set_ylabel(u"y");ax.set_zlabel(u"z")
ax.set_xlim([-5000,2000])
ax.set_ylim([-2000,5000])
ax.set_zlim([-5000,2000])
plt.show()

flight1.png

Code --ode

Dans scipy.integrate, il n'y a pas seulement odeint mais aussi une interface polyvalente pour le calcul numérique des équations différentielles ordinaires appelée ode, qui est orientée objet. Contrairement à odeint, vous pouvez spécifier la méthode de calcul, donc si vous voulez décider du contenu du calcul, c'est bien. De plus, il semble possible de terminer le calcul au milieu. Cependant, comme le calcul est effectué pendant, la vitesse de calcul est assez lente, en fonction du nombre de pas. C'était environ 10 fois plus lent dans mon environnement. Notez également que l'ordre des arguments de la fonction de l'équation différentielle est inversé par rapport à celui de (temps, variable) et odeint, et que ce doit être l'argument du tapple, il est donc légèrement différent.

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import ode
from scipy.linalg import block_diag
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.close('all')

def Rotation_X(Psi):
	R_x = np.array([[np.cos(Psi), np.sin(Psi), 0],
					[-np.sin(Psi), np.cos(Psi), 0],
					[0, 0, 1]])
	return R_x

def Rotation_Y(Theta):
	R_y = np.array([[np.cos(Theta), 0, -np.sin(Theta)],
					[0, 1, 0],
					[np.sin(Theta), 0, np.cos(Theta)]])
	return R_y

def Rotation_Z(Phi):
	R_z = np.array([[1, 0, 0],
					[0, np.cos(Phi), np.sin(Phi)],
					[0, -np.sin(Phi), np.cos(Phi)]])
	return R_z

#Équation du mouvement
#Contrairement au cas d'odeint(time, x, args*)Notez que les arguments sont tapuleux
def dynamical_system(t, x, (A, U0)):
	# x = [u,alpha,q,theta,beta,p,r,phi,psi,x,y,z]
	dx = A.dot(x)
	u = x[0]+U0 #la vitesse
	UVW = np.array([u, u*x[4], u*x[1]]) #Vecteur de vitesse[U,V,W]
	Rotation = Rotation_X(-x[8]).dot(Rotation_Y(-x[3])).dot(Rotation_Z(-x[7]))
	dX = Rotation.dot(UVW)
	dx[9]  = dX[0]
	dx[10] = dX[1]
	dx[11] = dX[2] 
	return dx

#Coefficient fin de stabilité dimensionnelle
Xu = -0.01;   Zu = -0.1;   Mu = 0.001
Xa = 30.0;    Za = -200.0; Ma = -4.0
Xq = 0.3;     Zq = -5.0;   Mq = -1.0
Yb = -45.0;   Lb_= -2.0;   Nb_= 1.0
Yp = 0.5;     Lp_= -1.0;   Np_= -0.1
Yr = 3.0;     Lr_= 0.2;    Nr_=-0.2

#Autres paramètres
W0 = 0.0;     U0 = 100.0;  theta0 = 0.05
g  = 9.8 #Accélération gravitationnelle

#Système vertical
A_lat = np.array([[Xu, Xa, -W0, -g*np.cos(theta0)],
    			  [Zu/U0, Za/U0, (U0+Zq)/U0, -g*np.sin(theta0)/U0],
    			  [Mu, Ma, Mq, 0],
    			  [0, 0, 1, 0]])

#Système latéral / directionnel
A_lon = np.array([[Yb, (W0+Yp), -(U0-Yr), g*np.cos(theta0), 0],
        		  [Lb_, Lp_, Lr_, 0, 0],
        		  [Nb_, Np_, Nr_, 0, 0],
        		  [0, 1, np.tan(theta0), 0, 0],
        		  [0, 0, 1/np.cos(theta0), 0, 0]])

#Combinez les systèmes sous forme de blocs diagonaux
A = block_diag(A_lat, A_lon)

#De plus, un espace sécurisé pour la trajectoire de vol
A = block_diag(A, np.zeros([3,3]))

#Définition des conditions de calcul
endurance	= 200	#Temps de vol[sec]
step		= 10	# 1.0[sec]Nombre de pas de temps par
dt 			= 1.0 / step 	#Largeur du pas de temps[sec]
t0 			= 0.0	#Heure de début[sec]
t = np.linspace(0,endurance,endurance*step)

#Valeur initiale x0= [u,alpha,q,theta, beta,p,r,phi,psi,x,y,z]
#Les variables à inclure dans l'équation différentielle normale doivent être unidimensionnelles.
x0_lat		= np.array([10, 0.1, 0.4, 0.2])		#Valeur initiale verticale
x0_lon		= np.array([0.0, 0.6, 0.4, 0.2, 0.2])	#Valeur initiale horizontale / direction
x0_pos		= np.array([0, 0, -1000]) #Position initiale de l'avion
x0 = np.hstack((x0_lat, x0_lon, x0_pos))

#Paramètres ODE
solver = ode(dynamical_system).set_integrator('vode', method='bdf')
solver.set_initial_value(x0, t0)
solver.set_f_params((A, U0))
dt = 1.0 / step
x = np.zeros([int(endurance/dt)+1, 12])
t = np.zeros([int(endurance/dt)+1])

index = 0
while solver.successful() and solver.t < endurance:
	solver.integrate(solver.t+dt)
	x[index] = solver.y
	t[index] = solver.t
	index += 1

print u"run successfully."

#Visualisation
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)

ax.plot(x[:,9], x[:,10], x[:,11])
ax.set_xlim([-5000,2000])
ax.set_ylim([-2000,5000])
ax.set_zlim([-5000,2000])
plt.show()

Recommended Posts

Analyse numérique des équations différentielles ordinaires avec l'odeint et l'ode de Scipy
L'histoire du calcul numérique des équations différentielles avec TensorFlow 2.0
Résolvez des équations différentielles normales simultanées avec Python et SymPy.
Trouvez la solution numérique de l'équation différentielle ordinaire du second ordre avec scipy
Résolution du problème de la valeur initiale des équations différentielles ordinaires avec JModelica
Calcul numérique d'équations aux dérivées partielles avec singularité (en utilisant les équations thermiques de type Hardy-Hénon à titre d'exemple)
Résolution d'équations différentielles normales avec Python ~ Gravitation universelle
[Calcul scientifique / technique par Python] Solution numérique d'équations différentielles ordinaires du premier ordre, problème de valeur initiale, calcul numérique
[Calcul scientifique / technique par Python] Solution numérique d'une équation différentielle ordinaire du second ordre, problème de valeur initiale, calcul numérique
[Méthode de calcul numérique, python] Résolution d'équations différentielles ordinaires par la méthode Eular
Effectuer une analyse isocurrent des canaux en eau libre avec Python et matplotlib
Analyse de correspondance des phrases avec l'API COTOHA et sauvegarde dans un fichier
Résoudre des équations différentielles normales en Python
Coexistence de Python2 et 3 avec CircleCI (1.0)
J'ai remplacé le calcul numérique de Python par Rust et comparé la vitesse
[Calcul scientifique / technique par Python] Résolution du problème de la valeur aux limites des équations différentielles ordinaires au format matriciel, calcul numérique