[PYTHON] Numerische Analyse gewöhnlicher Differentialgleichungen mit Scipys Odeint und Ode
Eine Methode zur Simulation (numerische Analyse / numerische Berechnung) eines physikalischen Modells in Form einer normalen Differentialgleichung wie eines Flugsimulators unter Verwendung von Python. Ich verwende odeint in Scipy.integrate. Es scheint, dass es den lsode von Odepack of Fortran verwendet, daher ist die Berechnung schnell.
Beispiel
Lassen Sie uns als Beispiel die Flugsimulation des Flugzeugs, das als Octave-Code (Matlab-kompatibel) veröffentlicht wurde, auf Python portieren. Der Link ist für den Inhalt der Berechnung detailliert. Butterfly_Effect( ) 6DOF Flight Simulation
Wenn Sie für verwenden, um etwas zu simulieren, das sich im Laufe der Zeit bewegt, ist selbst eine einfache numerische Berechnungsmethode wie die Euler-Methode in Python ziemlich langsam. Verwenden Sie daher so viel wie möglich eine Bibliothek wie scipy.
Unterschied zwischen Matlab und Python
Beachten Sie, dass sich der * -Operator der Matrix (ndarray) beim Portieren von matlab von matlab: inner product und numpy: element product unterscheidet. Beachten Sie außerdem, dass in Python, wenn Sie den Dezimalpunkt nicht angeben, es sich um eine Ganzzahldivision handelt.
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
#Bewegungsgleichung
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 #Geschwindigkeit
UVW = np.array([u, u*x[4], u*x[1]]) #Geschwindigkeitsvektor[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
#Dimensionsstabilität Feinkoeffizient
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
#Andere Parameter
W0 = 0.0; U0 = 100.0; theta0 = 0.05
g = 9.8 #Schwerkraftbeschleunigung
#Vertikales System
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]])
#Seiten- / Richtungssystem
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]])
#Kombinieren Sie Systeme als diagonale Blöcke
A = block_diag(A_lat, A_lon)
#Sichern Sie außerdem Platz für die Flugbahn
A = block_diag(A, np.zeros([3,3]))
#Berechnungsbedingungen einstellen
endurance = 200 #Flugzeit[sec]
step = 10 # 1.0[sec]Anzahl der Zeitschritte pro
t = np.linspace(0,endurance,endurance*step)
#Anfangswert x0= [u,alpha,q,theta, beta,p,r,phi,psi]
#Die Variablen, die in die normale Differentialgleichung aufgenommen werden sollen, müssen eindimensional sein.
x0_lat = np.array([10, 0.1, 0.4, 0.2]) #Vertikaler Anfangswert
x0_lon = np.array([0.0, 0.6, 0.4, 0.2, 0.2]) #Anfangswert Horizontal / Richtung
x0_pos = np.array([0, 0, -1000]) #Ausgangsposition des Flugzeugs
x0 = np.hstack((x0_lat, x0_lon, x0_pos))
x = odeint(dynamical_system, x0, t, args=(A,U0,))
print u"run successfully."
#Visualisierung
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()
Code --ode
In scipy.integrate gibt es nicht nur odeint, sondern auch eine universelle Schnittstelle zur numerischen Berechnung gewöhnlicher Differentialgleichungen namens ode, die objektorientiert ist. Im Gegensatz zu odeint können Sie die Berechnungsmethode angeben. Wenn Sie also den Inhalt der Berechnung festlegen möchten, ist dies gut. Es scheint auch möglich zu sein, die Berechnung in der Mitte zu beenden. Da die Berechnung jedoch währenddessen durchgeführt wird, ist die Berechnungsgeschwindigkeit abhängig von der Anzahl der Schritte ziemlich langsam. In meiner Umgebung war es ungefähr zehnmal langsamer. Beachten Sie auch, dass die Reihenfolge der Argumente der Funktion der Differentialgleichung von der von (Zeit, Variable) und Odeint umgekehrt ist und dass es sich um das Argument des Tapples handeln muss, sodass es sich geringfügig unterscheidet.
# -*- 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
#Bewegungsgleichung
#Im Gegensatz zu Odeint(time, x, args*)Beachten Sie, dass Argumente tapulous sind
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 #Geschwindigkeit
UVW = np.array([u, u*x[4], u*x[1]]) #Geschwindigkeitsvektor[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
#Dimensionsstabilität Feinkoeffizient
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
#Andere Parameter
W0 = 0.0; U0 = 100.0; theta0 = 0.05
g = 9.8 #Schwerkraftbeschleunigung
#Vertikales System
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]])
#Seiten- / Richtungssystem
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]])
#Kombinieren Sie Systeme als diagonale Blöcke
A = block_diag(A_lat, A_lon)
#Sichern Sie außerdem Platz für die Flugbahn
A = block_diag(A, np.zeros([3,3]))
#Berechnungsbedingungen einstellen
endurance = 200 #Flugzeit[sec]
step = 10 # 1.0[sec]Anzahl der Zeitschritte pro
dt = 1.0 / step #Breite des Zeitschritts[sec]
t0 = 0.0 #Startzeit[sec]
t = np.linspace(0,endurance,endurance*step)
#Anfangswert x0= [u,alpha,q,theta, beta,p,r,phi,psi,x,y,z]
#Die Variablen, die in die normale Differentialgleichung aufgenommen werden sollen, müssen eindimensional sein.
x0_lat = np.array([10, 0.1, 0.4, 0.2]) #Vertikaler Anfangswert
x0_lon = np.array([0.0, 0.6, 0.4, 0.2, 0.2]) #Anfangswert Horizontal / Richtung
x0_pos = np.array([0, 0, -1000]) #Ausgangsposition des Flugzeugs
x0 = np.hstack((x0_lat, x0_lon, x0_pos))
#ODE-Einstellungen
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."
#Visualisierung
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