[PYTHON] Doppelpendelsimulation
Was ist ein Doppelpendel?
https://www.youtube.com/watch?v=25feOUNQB2Y Dies ist ein berühmtes Experiment als Chaos-Phänomen. Das Chaos-Phänomen ist ein Phänomen, bei dem der Ausgangszustand einen großen Einfluss auf den nachfolgenden Zustand hat. Die Geschichte des Schmetterlingseffekts, bei dem die Flügel eines Schmetterlings aufgrund des Luftstroms einen Taifun erzeugen, wird auch in der Erklärung des gleichen Chaosphänomens erwähnt. (Kann nicht korrekt sein)
Ableitung und Simulation der kinetischen Gleichung mit der von [hier] erstellten Euler-Lagrange-Methode (https://qiita.com/hiRina/items/519d5168c2c2fdf012cd).
Variablendefinition
Stellen Sie das Problem im Qualitätspunktesystem / zwei Dimensionen ein.
Die Geschwindigkeit von m1 und die Geschwindigkeit von m2
\vec{R_{m_1}} =\left(
\begin{array}{ccc}
r_1 sin(\theta_1) \\
-r_1 cos(\theta_1)
\end{array}
\right) ,
\vec{R_{m_2}}=\vec{R_{m_1}} +\left(
\begin{array}{ccc}
r_2 sin(\theta_1) \\
-r_2 cos(\theta_1)
\end{array}
\right) \\
\vec{V_{m_1}} =\left(
\begin{array}{ccc}
r_1 cos(\theta_1) \\
r_1 sin(\theta_1)
\end{array}
\right) \dot{\theta_1},
\vec{V_{m_2}}=\vec{V_{m_1}} +\left(
\begin{array}{ccc}
r_2 cos(\theta_1) \\
r_2 sin(\theta_1)
\end{array}
\right) \dot{\theta_2}\\
Übungsenergie T.
T=\frac{1}{2} m_1 |\vec{V_{m_1}}|^2+\frac{1}{2} m_2 |\vec{V_{m_2}}|^2
Position Energie U.
U=m_1gr_1(1-cos(\theta_1)) + m_2g(r_1(1-cos(\theta_1)) + r_2(1-cos(\theta_2))
Wenn Sie Lagrange berechnen
L= m_1*r_1^2/2 + m2*r_2^2 /2 \dot{ x_1}^{2.0} +
m_2*r_2^2/2 \dot{ x_2}^{2.0} +
m_2*r_1*r_2 \dot{ x_1}^{1.0} sin( x_1)^{1.0} \dot{ x_2}^{1.0} sin( x_2)^{1.0} +
m_2*r_1*r_2 \dot{ x_1}^{1.0} cos( x_1)^{1.0} \dot{ x_2}^{1.0} cos( x_2)^{1.0} +
-1*(g*(m_1*r_1+m_2*(r_1+r_2))) +
-1*(-1*g*(m_1*r_1+m_2*r_1)) cos( x_1)^{1.0} +
-1*(-1*g*r_2*m_2) cos( x_2)^{1.0}
\\
\theta_1 = x_1\\
\theta_2 = x_2
Bei der Ableitung der Bewegungsgleichung
(m_2*r_1*r_2*1.0) + -1*( ( (m_2*r_1*r_2*1.0) *1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0} \dot{ x_2}^{1.0} sin( x_2)^{1.0} +
(m_2*r_1*r_2*(-1)*1.0) + -1*( ( (m_2*r_1*r_2*1.0) *(-1)*1.0) ) \dot{ x_1}^{1.0} sin( x_1)^{1.0} \dot{ x_2}^{1.0} cos( x_2)^{1.0} +
(-1*(-1*g*(m_1*r_1+m_2*r_1))*(-1)*1.0) sin( x_1)^{1.0} +
-1*( ( (m_1*r_1^2/2 + m2*r_2^2 /2*2.0) *1.0) ) \ddot { x_1}^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) sin( x_1)^{1.0} \dot{ x_2}^{2.0} cos( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *(-1)*1.0) ) cos( x_1)^{1.0} \dot{ x_2}^{2.0} sin( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) sin( x_1)^{1.0} \ddot { x_2}^{1.0} sin( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) cos( x_1)^{1.0} \ddot { x_2}^{1.0} cos( x_2)^{1.0} = 0 \\
(m_2*r_1*r_2*1.0) + -1*( ( (m_2*r_1*r_2*1.0) *1.0) ) \dot{ x_1}^{1.0} sin( x_1)^{1.0} \dot{ x_2}^{1.0} cos( x_2)^{1.0} +
(m_2*r_1*r_2*(-1)*1.0) + -1*( ( (m_2*r_1*r_2*1.0) *(-1)*1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0} \dot{ x_2}^{1.0} sin( x_2)^{1.0} +
(-1*(-1*g*r_2*m_2)*(-1)*1.0) sin( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) \dot{ x_1}^{2.0} cos( x_1)^{1.0} sin( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *(-1)*1.0) ) \dot{ x_1}^{2.0} sin( x_1)^{1.0} cos( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) \ddot { x_1}^{1.0} sin( x_1)^{1.0} sin( x_2)^{1.0} +
-1*( ( (m_2*r_1*r_2*1.0) *1.0) ) \ddot { x_1}^{1.0} cos( x_1)^{1.0} cos( x_2)^{1.0} +
-1*( ( (m_2*r_2^2/2*2.0) *1.0) ) \ddot { x_2}^{1.0} =0
Der Code, den ich verwendet habe, ist schwer zu erkennen, daher werde ich ihn am Ende anhängen.
Simulation (Euler-Methode)
Bei Simulation nach der Euler-Methode unter Verwendung der erhaltenen Formel Aufbau
m1=1.0\\
m2=1.0\\
r1=1.0\\
r2=1.0\\
g=9.81\\
Ursprünglicher Wert:
\theta_1 = 180deg\\
\theta_2 = 180deg\\
\dot{\theta_1 }= 2rad/s\\
\dot{\theta_2} = 2rad/s\\
Ursprünglicher Wert:
\theta_1 = 180deg\\
\theta_2 = 180deg\\
\dot{\theta_1 }= 0.1rad/s\\
\dot{\theta_2} = 0.1rad/s
Dieser ist chaotischer (mittleres 2 durchschnittliches Gefühl)
\theta_1 = 45deg\\
\theta_2 = 45deg\\
\dot{\theta_1 }= 0.1rad/s\\
\dot{\theta_2} = 0.1rad/s
Ohne Geschwindigkeit kann man nicht chaotisch sein
Bewegungsgleichungs-Ableitungscode
Ab hier ist es das Programm, das Sie erstellt haben. Es sollte mit Kopieren und Einfügen funktionieren. Ich laufe auf Anaconda Python2. gif Klicken Sie hier, um zu animieren https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2 Ich bezog mich auf.
Analyseklasse
Formula2
import numpy as np #Numpy Bibliothek
import copy
import math
'''
Formelverarbeitungsprogramm
Spalte: Polygon
Zeile: θ θ θ θ ・ sin θ cos θ Koeffizient
6 Stück x Anzahl der Variablen
Entspricht dem Zeichenausdruck
Jedoch,
'''
class Formula2 :
def __init__(self,L,num,coeff):
self.L = L
self.n = num #Anzahl der Variablentypen
self.coeff = coeff
#print self.coeff.dtype
def __add__(self,other):
self.L = np.append(self.L,other.L,axis =0)
#print self.coeff.dtype
#print other.coeff.dtype
expr,vari = self.L.shape
self.coeff = np.append(self.coeff,other.coeff,axis =0)
self.together()
#self.erase()
def __sub__(self,other):
expr,vari = other.L.shape
hoge = copy.deepcopy(other.L)
for i in range(expr):
hoge[i][5] = hoge[i][5] * -1
self.L = np.append(self.L,hoge,axis =0)
for i in range(expr):
other.coeff[i] = "-1*(" + other.coeff[i] +")"
self.coeff = np.append(self.coeff,other.coeff,axis =0)
self.together()
#self.erase()
def together(self): #Fassen Sie die Koeffizienten desselben Terms zusammen
expr,vari = self.L.shape
if expr >=2:
for i in range(expr-1):
for j in range(i+1,expr):
m=0
while self.L[i][0+m*6:5+m*6].tolist() == self.L[j][0+m*6:5+m*6].tolist():
m += 1
if m == self.n:
break
if m== self.n: #Alle stimmen mit der Addition von Koeffizienten überein
self.L[i][5] += self.L[j][5]
self.coeff[i] = self.coeff[i] +" + "+ self.coeff[j]
for k in range(vari): #j Zeile gelöscht
self.L[j][k] = 0
self.coeff[j] = "0"
def erase(self): #Löschen Sie den Term mit einem Koeffizienten von 0
expr,vari = self.L.shape
for i in list(reversed(range(expr))):
if self.coeff[i] == "0":
self.L = np.delete(self.L,i,axis = 0)
self.coeff =np.delete(self.coeff,i,axis = 0)
def partial(self,moji,kind): #Unterscheiden Sie Begriffe teilweise mit Art
expr,vari = self.L.shape
if kind == 0:
'''
sin
cos
+moji*6
'''
#print self.L
for i in range(expr):
if self.L[i][3+moji*6] !=0: #sin
hoge = copy.deepcopy(self.L[i])
hogestr = copy.deepcopy(self.coeff[i])
hoge[5] = hoge[5]*hoge[3+moji*6]
#Hoge auf Koeffizienten[3+moji*6]multiplizieren
hogestr = " ("+hogestr+"*"+str(hoge[3+moji*6])+") "
hoge[3+moji*6] = hoge[3+moji*6]-1
hoge[4+moji*6] = hoge[4+moji*6]+1
self.L = np.append(self.L,hoge[np.newaxis,:],axis =0)
self.coeff = np.append(self.coeff,[hogestr],axis =0)
#print self.L
for i in range(expr):
if self.L[i][4+moji*6] !=0: #cos
hoge = copy.deepcopy(self.L[i])
hogestr = copy.deepcopy(self.coeff[i])
hoge[5] = hoge[5]*hoge[4+moji*6]*-1
#Hoge auf Koeffizienten[3+moji*6]multiplizieren
hogestr = " ("+hogestr+"*(-1)*"+str(hoge[4+moji*6])+") "
hoge[4+moji*6]= hoge[4+moji*6]-1
hoge[3+moji*6] = hoge[3+moji*6]+1
self.L = np.append(self.L,hoge[np.newaxis,:],axis =0)
self.coeff = np.append(self.coeff,[hogestr],axis =0)
#print self.L
'''
Wie es ist
Verringern Sie die Reihenfolge um eins, wenn 0, nichts tun
Multiplizieren Sie die ursprüngliche Reihenfolge mit dem Koeffizienten
sin,Es ist in Ordnung, weil der cos-Verarbeitungsteil hinter und nicht im Verarbeitungsbereich liegt
'''
for i in range(expr):
if self.L[i][kind] !=0:
#print kind,self.L[i][5]
self.L[i][5] = self.L[i][5] * self.L[i][kind+moji*6]
self.coeff[i] = " ("+self.coeff[i]+"*"+str(self.L[i][kind+moji*6])+") "
#print self.L[i][5]
self.L[i][kind+moji*6] = self.L[i][kind+moji*6] -1
else:
self.L[i][5] = 0 #0 wenn nicht enthalten
self.coeff[i] = "0"
if kind == 1:
'''
Wie es ist
Verringern Sie die Reihenfolge um eins, wenn 0, nichts tun
Multiplizieren Sie die ursprüngliche Reihenfolge mit dem Koeffizienten
'''
for i in range(expr):
if self.L[i][kind+moji*6] !=0:
self.L[i][5] = self.L[i][5] * self.L[i][kind+moji*6]
self.coeff[i] = " ("+self.coeff[i]+"*"+str(self.L[i][kind+moji*6])+") "
self.L[i][kind+moji*6] = self.L[i][kind+moji*6] -1
else:
self.L[i][5] = 0 #0 wenn nicht enthalten
self.coeff[i] = "0"
self.together()
self.erase()
def diff(self): #Zeitdifferenzierung
'''
Nach teilweiser Differenzierung mit 4 Mustern
Erhöhen Sie die Reihenfolge von θ
'''
L1=copy.deepcopy(self.L)
L1_coeff = copy.deepcopy(self.coeff)
for i in range(self.n):
self.L =copy.deepcopy( L1) #L1:Vor der Differenzierung L2:Zum Speichern von Ergebnissen
self.coeff = copy.deepcopy(L1_coeff)
self.partial(i,0)
expr,vari = self.L.shape
#print expr
#print self.L
for j in range(expr):
self.L[j][1+i*6] = self.L[j][1+i*6] + 1
if i == 0:
L2 = copy.deepcopy( self.L)
L2_coeff = copy.deepcopy(self.coeff)
else:
L2 = np.append(L2,self.L,axis =0)
L2_coeff = np.append(L2_coeff,self.coeff,axis =0)
self.L =copy.deepcopy( L1)
self.coeff = copy.deepcopy(L1_coeff)
self.partial(i,1)
expr,vari = self.L.shape
for j in range(expr):
self.L[j][2+i*6] += 1
L2 = np.append(L2,self.L,axis =0)
L2_coeff = np.append(L2_coeff,self.coeff,axis =0)
self.L = copy.deepcopy( L2)
self.coeff = copy.deepcopy(L2_coeff)
self.together()
self.erase()
def disp2(self): #Polymeranzeige
print "-------------------Formula--------------"
expr,vari = self.L.shape
for j in range(expr):
hoge =""
for i in range(self.n):
hoge += "["+str(i+1)+" "
if i == 0:
hoge += self.coeff[j]+" "
if self.L[j][0+i*6]!= 0:
hoge += " θ^"
hoge += str(self.L[j][0+i*6])
if self.L[j][1+i*6]!= 0:
hoge += "θ ・^"
hoge += str(self.L[j][1+i*6])
if self.L[j][2+i*6]!= 0:
hoge += "θ ...^"
hoge += str(self.L[j][2+i*6])
if self.L[j][3+i*6]!= 0:
hoge += " sin^"
hoge += str(self.L[j][3+i*6])
if self.L[j][4+i*6]!= 0:
hoge += " cos^"
hoge += str(self.L[j][4+i*6])
hoge += " ] "
hoge += " + "
print hoge
#print self.coeff
def latex(self):
variable = " x"
k=1
print "-------------------Formula_latex---------"
expr,vari = self.L.shape
for j in range(expr):
hoge =""
for i in range(self.n):
hoge += " "
if i == 0:
hoge += self.coeff[j]
if self.L[j][0+i*6]!= 0:
hoge += variable+"_"+str(i+k)+"^{"
hoge += str(self.L[j][0+i*6])+"}"
if self.L[j][1+i*6]!= 0:
hoge += " \dot{"+variable+"_"+str(i+k)+"}^{"
hoge += str(self.L[j][1+i*6])+"}"
if self.L[j][2+i*6]!= 0:
hoge += " \ddot {"+variable+"_"+str(i+k)+"}^{"
hoge += str(self.L[j][2+i*6])+"}"
if self.L[j][3+i*6]!= 0:
hoge += " sin("+variable+"_"+str(i+k)+")^{"
hoge += str(self.L[j][3+i*6])+"}"
if self.L[j][4+i*6]!= 0:
hoge += " cos("+variable+"_"+str(i+k)+")^{"
hoge += str(self.L[j][4+i*6])+"}"
hoge += " "
hoge += " + "
print hoge
def separete(self): #Stellen Sie sich das als 2 Variablen vor und teilen Sie es in 3 Teile, 1 und 2b.
a1=np.asarray([[]])
a2=np.asarray([[]])
coeff1 = np.asarray([],dtype = object)
coeff2 = np.asarray([],dtype = object)
expr,vari = self.L.shape
for j in range(expr)[::-1]:
if self.L[j][2] == 1:
if len(a1[0]) == 0:
a1=(self.L[j]) [np.newaxis,:]
coeff1=np.asarray([self.coeff[j]],dtype = object)
else:
a1 = np.append(a1,(self.L[j]) [np.newaxis,:],axis =0)
coeff1=np.append(coeff1,[self.coeff[j]],axis = 0)
self.L = np.delete(self.L,j,0)
self.coeff=np.delete(self.coeff,j,0)
continue
if self.L[j][8] == 1:
if len(a2[0]) == 0:
a2=(self.L[j]) [np.newaxis,:]
coeff2=np.asarray([self.coeff[j]],dtype = object)
else:
a2 = np.append(a2,(self.L[j]) [np.newaxis,:],axis =0)
coeff2=np.append(coeff2,[self.coeff[j]],axis = 0)
self.L = np.delete(self.L,j,0)
self.coeff=np.delete(self.coeff,j,0)
continue
if len(a1[0]) == 0:
a1 =np.asarray([ [0,0,0,0,0,0 ,0,0,0,0,0,0]])
coeff1 = np.asarray(["0"],dtype = object)
if len(a2[0]) == 0:
a2 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0]])
coeff2 = np.asarray(["0"],dtype = object)
L1 = Formula2(a1,2,coeff1)
L2 = Formula2(a2,2,coeff2)
return L1,L2
def separete3(self): #Stellen Sie es sich als 3 Variablen vor und teilen Sie es in 4 Teile: 1 2 3 b
a1=np.asarray([[]])
a2=np.asarray([[]])
a3=np.asarray([[]])
coeff1 = np.asarray([],dtype = object)
coeff2 = np.asarray([],dtype = object)
coeff3 = np.asarray([],dtype = object)
expr,vari = self.L.shape
for j in range(expr)[::-1]:
if self.L[j][2] == 1:
if len(a1[0]) == 0:
a1=(self.L[j]) [np.newaxis,:]
coeff1=np.asarray([self.coeff[j]],dtype = object)
else:
a1 = np.append(a1,(self.L[j]) [np.newaxis,:],axis =0)
coeff1=np.append(coeff1,[self.coeff[j]],axis = 0)
self.L = np.delete(self.L,j,0)
self.coeff=np.delete(self.coeff,j,0)
continue
if self.L[j][8] == 1:
if len(a2[0]) == 0:
a2=(self.L[j]) [np.newaxis,:]
coeff2=np.asarray([self.coeff[j]],dtype = object)
else:
a2 = np.append(a2,(self.L[j]) [np.newaxis,:],axis =0)
coeff2=np.append(coeff2,[self.coeff[j]],axis = 0)
self.L = np.delete(self.L,j,0)
self.coeff=np.delete(self.coeff,j,0)
continue
if self.L[j][14] == 1:
if len(a3[0]) == 0:
a3=(self.L[j]) [np.newaxis,:]
coeff3=np.asarray([self.coeff[j]],dtype = object)
else:
a3 = np.append(a3,(self.L[j]) [np.newaxis,:],axis =0)
coeff3=np.append(coeff3,[self.coeff[j]],axis = 0)
self.L = np.delete(self.L,j,0)
self.coeff=np.delete(self.coeff,j,0)
continue
if len(a1[0]) == 0:
a1 =np.asarray([ [0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
coeff1 = np.asarray(["0"],dtype = object)
if len(a2[0]) == 0:
a2 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
coeff2 = np.asarray(["0"],dtype = object)
if len(a3[0]) == 0:
a3 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
coeff3 = np.asarray(["0"],dtype = object)
L1 = Formula2(a1,3,coeff1)
L2 = Formula2(a2,3,coeff2)
L3 = Formula2(a3,3,coeff3)
return L1,L2,L3
if __name__ == '__main__':
print "This is not for main TEST RESULT"
T=np.asarray([[1,2,0,4,5,1 ,1,2,0,4,5,1]])
coeff = np.asarray(["m "], dtype=object)
print coeff.dtype
sai =Formula2(T,2,coeff)
sai.disp2()
sai.diff()
sai.disp2()
#sai.latex()
L1,L2= sai.separete()
L1.disp2()
'''
L1+L2
L1.disp2()
'''
Ausführungscode
import numpy as np #Numpy Bibliothek
import copy
import math
#import matplotlib.pyplot as plt
import matplotlib
#matplotlib.use('Agg') # -----(1)
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
import Formula2
#L= np.asarray([[7,2,3,4,5,6, 7,2,3,4,5,6]])
#L= np.asarray([[2,2,0,5,4,5],[2,2,0,5,4,5]])
'''
L1,Betrachten wir den Fall, in dem es nur zwei Ausdrücke von L2 und zwei Variablen gibt.
Als Formel ≒ 0
Extrahiere nur den Gegenstand von θ ... und mache eine neue Instanziierung
Der entfernte wird als b bezeichnet[L1_b,L2_b]Machen Sie einen Vektor
Der Term von θ ...
[L1_1 L1_2
L2_1 L2_2 ]Betrachten Sie die Instanzmatrix A.
L1_2 =Punkt von θ2 in der Formel von L1
θ ...= A^-Berechnen Sie 1 b
Θ ・ θ erscheint auch in der expliziten Methode der Euler-Methode
Lass uns das machen
L1_b,L1_1,L1_Überlegen Sie, wie Sie 2 generieren
'''
m1=1.0
m2=1.0
r1=1.0
r2=1.0
g=9.81
T=np.asarray([[0,2,0,0,0,m1*r1**2/2 + m2*r2**2 /2,
0,0,0,0,0,1],
[0,0,0,0,0,m2*r2**2/2 ,
0,2,0,0,0,1],
[0,1,0,1,0,m2*r1*r2,
0,1,0,1,0,1],
[0,1,0,0,1,m2*r1*r2,
0,1,0,0,1,1]]
)
T_n=np.asarray(["m_1*r_1^2/2 + m2*r_2^2 /2",
"m_2*r_2^2/2" ,
"m_2*r_1*r_2",
"m_2*r_1*r_2"],dtype=object )
U= np.asarray([[0,0,0,0,0,g*(m1*r1+m2*(r1+r2)),
0,0,0,0,0,1],
[0,0,0,0,1,-1*g*(m1*r1+m2*r1),
0,0,0,0,0,1],
[0,0,0,0,0,-1*g*r2*m2,
0,0,0,0,1,1]
])
U_n= np.asarray(["g*(m_1*r_1+m_2*(r_1+r_2))",
"-1*g*(m_1*r_1+m_2*r_1)",
"-1*g*r_2*m_2"],dtype=object)
T1=Formula2.Formula2(T,2,T_n)
SaveT = copy.deepcopy(T1)
U1=Formula2.Formula2(U,2,U_n)
T1-U1
Lt_1=copy.deepcopy(T1)
Lt_2=copy.deepcopy(T1)
Lt_1.partial(0,0)
Lt_2.partial(0,1)
Lt_2.diff()
Lt_1-Lt_2
Ls_1=copy.deepcopy(T1)
Ls_2=copy.deepcopy(T1)
Ls_1.partial(1,0)
Ls_2.partial(1,1)
Ls_2.diff()
Ls_1-Ls_2
Lt_1.disp2()
Ls_1.disp2()
T1.disp2()
L1 = copy.deepcopy(Ls_1)
L2 = copy.deepcopy(Lt_1)
L1_1,L1_2 = L1.separete()
L2_1,L2_2 = L2.separete()
'''
T1:Lagrange
Lt_1 ,L1 :Bewegungsgleichung für θ1
Ls_2 ,L2 :Bewegungsgleichung für θ2
'''
Simulationscode
Analyseklasse
Formula.py
import numpy as np #Numpy Bibliothek
import copy
import math
'''
Formelverarbeitungsprogramm
Spalte: Polygon
Zeile: θ θ θ θ ・ sin θ cos θ Koeffizient
6 Stück x Anzahl der Variablen
'''
class Formula :
def __init__(self,L,num):
self.L = L
self.n = num #Anzahl der Variablentypen
def __add__(self,other):
self.L = np.append(self.L,other.L,axis =0)
self.together()
self.erase()
def __sub__(self,other):
expr,vari = other.L.shape
hoge = copy.deepcopy(other.L)
for i in range(expr):
hoge[i][5] = hoge[i][5] * -1
self.L = np.append(self.L,hoge,axis =0)
self.together()
self.erase()
def latex(self):
variable = " x"
k=1
print "-------------------Formula--------------"
expr,vari = self.L.shape
for j in range(expr):
hoge =""
for i in range(self.n):
hoge += " "
if self.L[j][5]!= 0:
if i == 0:
hoge += str(self.L[j][5+i*6])
if self.L[j][0+i*6]!= 0:
hoge += variable+"_"+str(i+k)+"^{"
hoge += str(self.L[j][0+i*6])+"}"
if self.L[j][1+i*6]!= 0:
hoge += " \dot{"+variable+"_"+str(i+k)+"}^{"
hoge += str(self.L[j][1+i*6])+"}"
if self.L[j][2+i*6]!= 0:
hoge += " \ddot {"+variable+"_"+str(i+k)+"}^{"
hoge += str(self.L[j][2+i*6])+"}"
if self.L[j][3+i*6]!= 0:
hoge += " sin("+variable+"_"+str(i+k)+")^{"
hoge += str(self.L[j][3+i*6])+"}"
if self.L[j][4+i*6]!= 0:
hoge += " cos("+variable+"_"+str(i+k)+")^{"
hoge += str(self.L[j][4+i*6])+"}"
hoge += " "
hoge += " + "
print hoge
def together(self): #Fassen Sie die Koeffizienten desselben Terms zusammen
expr,vari = self.L.shape
if expr >=2:
for i in range(expr-1):
for j in range(i+1,expr):
m=0
while self.L[i][0+m*6:5+m*6].tolist() == self.L[j][0+m*6:5+m*6].tolist():
m += 1
if m == self.n:
break
if m== self.n:
self.L[i][5] += self.L[j][5]
for k in range(vari):
self.L[j][k] = 0
def erase(self): #Löschen Sie den Term mit einem Koeffizienten von 0
expr,vari = self.L.shape
for i in list(reversed(range(expr))):
if self.L[i][5] == 0:
self.L = np.delete(self.L,i,axis = 0)
def partial(self,moji,kind): #Unterscheiden Sie Begriffe teilweise mit Art
expr,vari = self.L.shape
if kind == 0:
'''
sin
cos
+moji*6
'''
#print self.L
for i in range(expr):
if self.L[i][3+moji*6] !=0: #sin
hoge = copy.deepcopy(self.L[i])
hoge[5] = hoge[5]*hoge[3+moji*6]
hoge[3+moji*6] = hoge[3+moji*6]-1
hoge[4+moji*6] = hoge[4+moji*6]+1
self.L = np.append(self.L,hoge[np.newaxis,:],axis =0)
#print self.L
for i in range(expr):
if self.L[i][4+moji*6] !=0: #cos
hoge = copy.deepcopy(self.L[i])
hoge[5] = hoge[5]*hoge[4+moji*6]*-1
hoge[4+moji*6]= hoge[4+moji*6]-1
hoge[3+moji*6] = hoge[3+moji*6]+1
self.L = np.append(self.L,hoge[np.newaxis,:],axis =0)
#print self.L
'''
Wie es ist
Verringern Sie die Reihenfolge um eins, wenn 0, nichts tun
Multiplizieren Sie die ursprüngliche Reihenfolge mit dem Koeffizienten
'''
#print self.L
for i in range(expr):
if self.L[i][kind] !=0:
#print kind,self.L[i][5]
self.L[i][5] = self.L[i][5] * self.L[i][kind+moji*6]
#print self.L[i][5]
self.L[i][kind+moji*6] = self.L[i][kind+moji*6] -1
else:
self.L[i][5] = 0 #0 wenn nicht enthalten
if kind == 1:
'''
Wie es ist
Verringern Sie die Reihenfolge um eins, wenn 0, nichts tun
Multiplizieren Sie die ursprüngliche Reihenfolge mit dem Koeffizienten
'''
for i in range(expr):
if self.L[i][kind+moji*6] !=0:
self.L[i][5] = self.L[i][5] * self.L[i][kind+moji*6]
self.L[i][kind+moji*6] = self.L[i][kind+moji*6] -1
else:
self.L[i][5] = 0 #0 wenn nicht enthalten
self.together()
self.erase()
def diff(self): #Zeitdifferenzierung
'''
Nach teilweiser Differenzierung mit 4 Mustern
Erhöhen Sie die Reihenfolge von θ
'''
L1=copy.deepcopy(self.L)
for i in range(self.n):
self.L =copy.deepcopy( L1)
self.partial(i,0)
expr,vari = self.L.shape
#print expr
#print self.L
for j in range(expr):
self.L[j][1+i*6] = self.L[j][1+i*6] + 1
if i == 0:
L2 = copy.deepcopy( self.L)
else:
L2 = np.append(L2,self.L,axis =0)
self.L =copy.deepcopy( L1)
self.partial(i,1)
expr,vari = self.L.shape
for j in range(expr):
self.L[j][2+i*6] += 1
L2 = np.append(L2,self.L,axis =0)
self.L = copy.deepcopy( L2)
self.together()
self.erase()
def disp(self): #Polymeranzeige
print "-------------------Formula--------------"
expr,vari = self.L.shape
for j in range(expr):
hoge =""
for i in range(self.n):
hoge += str(self.L[j][5+i*6])
hoge += " θ^"
hoge += str(self.L[j][0+i*6])
hoge += "θ ・^"
hoge += str(self.L[j][1+i*6])
hoge += "θ ...^"
hoge += str(self.L[j][2+i*6])
hoge += " sin^"
hoge += str(self.L[j][3+i*6])
hoge += " cos^"
hoge += str(self.L[j][4+i*6])
hoge += " "
hoge += " + "
print hoge
def disp2(self): #Polymeranzeige
print "-------------------Formula--------------"
expr,vari = self.L.shape
for j in range(expr):
hoge =""
for i in range(self.n):
hoge += "["+str(i+1)+" "
if self.L[j][5]!= 0:
if i == 0:
hoge += str(self.L[j][5+i*6])
if self.L[j][0+i*6]!= 0:
hoge += " θ^"
hoge += str(self.L[j][0+i*6])
if self.L[j][1+i*6]!= 0:
hoge += "θ ・^"
hoge += str(self.L[j][1+i*6])
if self.L[j][2+i*6]!= 0:
hoge += "θ ...^"
hoge += str(self.L[j][2+i*6])
if self.L[j][3+i*6]!= 0:
hoge += " sin^"
hoge += str(self.L[j][3+i*6])
if self.L[j][4+i*6]!= 0:
hoge += " cos^"
hoge += str(self.L[j][4+i*6])
hoge += " ] "
hoge += " + "
print hoge
def subst(self,x): #Substitution x,In der Reihenfolge von x setzen
vari =np.array(x).shape
if vari[0] != self.n*2:
print "cannot subst"
else:
'''
Generieren Sie einen Zeilenvektor
Machen Sie L mit Punkt und fügen Sie alle Elemente des resultierenden Vektors hinzu
'''
expr,vari = self.L.shape
sum1 = 0
hoge=0
for j in range(expr):
hoge=self.L[j][5] #Koeffizient
for i in range(self.n):
hoge = hoge*x[0+i*2]**self.L[j][0+i*6] #θ
hoge = hoge*x[1+i*2]**self.L[j][1+i*6] #θ ・
hoge = hoge*math.sin(x[0+i*2])**self.L[j][3+i*6] #sin
hoge = hoge*math.cos(x[0+i*2])**self.L[j][4+i*6] #cos
sum1 += hoge
return sum1
def separete(self): #Stellen Sie sich das als 2 Variablen vor und teilen Sie es in 3 Teile, 1 und 2b.
a1=np.asarray([[]])
a2=np.asarray([[]])
expr,vari = self.L.shape
for j in range(expr)[::-1]:
if self.L[j][2] == 1:
if len(a1[0]) == 0:
a1=(self.L[j]) [np.newaxis,:]
else:
a1 = np.append(a1,(self.L[j]) [np.newaxis,:],axis =0)
self.L = np.delete(self.L,j,0)
continue
if self.L[j][8] == 1:
if len(a2[0]) == 0:
a2=(self.L[j]) [np.newaxis,:]
else:
a2 = np.append(a2,(self.L[j]) [np.newaxis,:],axis =0)
self.L = np.delete(self.L,j,0)
continue
if len(a1[0]) == 0:
a1 =np.asarray([ [0,0,0,0,0,0 ,0,0,0,0,0,0]])
if len(a2[0]) == 0:
a2 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0]])
L1 = Formula(a1,2)
L2 = Formula(a2,2)
return L1,L2
def separete3(self): #Stellen Sie es sich als 3 Variablen vor und teilen Sie es in 4 Teile: 1 2 3 b
a1=np.asarray([[]])
a2=np.asarray([[]])
a3=np.asarray([[]])
expr,vari = self.L.shape
for j in range(expr)[::-1]:
print "a1"
print a1.shape
print "a3"
print a3.shape
print "a2"
print a2.shape
if self.L[j][2] == 1:
if len(a1[0]) == 0:
a1=(self.L[j]) [np.newaxis,:]
else:
a1 = np.append(a1,(self.L[j]) [np.newaxis,:],axis =0)
self.L = np.delete(self.L,j,0)
continue
if self.L[j][8] == 1:
if len(a2[0]) == 0:
a2=(self.L[j]) [np.newaxis,:]
else:
a2 = np.append(a2,(self.L[j]) [np.newaxis,:],axis =0)
self.L = np.delete(self.L,j,0)
continue
if self.L[j][14] == 1:
if len(a3[0]) == 0:
a3=(self.L[j]) [np.newaxis,:]
else:
a3 = np.append(a3,(self.L[j]) [np.newaxis,:],axis =0)
self.L = np.delete(self.L,j,0)
continue
if len(a1[0]) == 0:
a1 =np.asarray([ [0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
if len(a2[0]) == 0:
a2 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
if len(a3[0]) == 0:
a3 = np.asarray([[0,0,0,0,0,0 ,0,0,0,0,0,0 ,0,0,0,0,0,0]])
L1 = Formula(a1,3)
L2 = Formula(a2,3)
L3 = Formula(a3,3)
return L1,L2,L3
if __name__ == '__main__':
print "This is not for main"
Ausführungsprogramm
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 21 14:11:03 2020
@author: kisim
"""
import numpy as np #Numpy Bibliothek
import copy
import math
#import matplotlib.pyplot as plt
import matplotlib
#matplotlib.use('Agg') # -----(1)
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
import Formula
#L= np.asarray([[7,2,3,4,5,6, 7,2,3,4,5,6]])
#L= np.asarray([[2,2,0,5,4,5],[2,2,0,5,4,5]])
'''
L1,Betrachten wir den Fall, in dem es nur zwei Ausdrücke von L2 und zwei Variablen gibt.
Als Formel ≒ 0
Extrahiere nur den Gegenstand von θ ... und mache eine neue Instanziierung
Der entfernte wird als b bezeichnet[L1_b,L2_b]Machen Sie einen Vektor
Der Term von θ ...
[L1_1 L1_2
L2_1 L2_2 ]Betrachten Sie die Instanzmatrix A.
L1_2 =Punkt von θ2 in der Formel von L1
θ ...= A^-Berechnen Sie 1 b
Θ ・ θ erscheint auch in der expliziten Methode der Euler-Methode
Lass uns das machen
L1_b,L1_1,L1_Überlegen Sie, wie Sie 2 generieren
'''
m1=1.0
m2=1.0
l1=1.0
l2=1.0
r1=1.0
r2=1.0
R1=1.0
g=9.81
T=np.asarray([[0,2,0,0,0,m1*r1**2/2 + m2*r2**2 /2,
0,0,0,0,0,1],
[0,0,0,0,0,m2*r2**2/2 ,
0,2,0,0,0,1],
[0,1,0,1,0,m2*r1*r2,
0,1,0,1,0,1],
[0,1,0,0,1,m2*r1*r2,
0,1,0,0,1,1]]
)
T_n=np.asarray(["m_1*r_1^2/2 + m2*r_2^2 /2",
"m_2*r_2^2/2" ,
"m_2*r_1*r_2",
"m_2*r_1*r_2"],dtype=object )
U= np.asarray([[0,0,0,0,0,g*(m1*r1+m2*(r1+r2)),
0,0,0,0,0,1],
[0,0,0,0,1,-1*g*(m1*r1+m2*r1),
0,0,0,0,0,1],
[0,0,0,0,0,-1*g*r2*m2,
0,0,0,0,1,1]
])
U_n= np.asarray(["g*(m_1*r_1+m_2*(r_1+r_2))",
"-1*g*(m_1*r_1+m_2*r_1)",
"-1*g*r_2*m_2"],dtype=object)
T1=Formula.Formula(T,2)
U1=Formula.Formula(U,2)
T1-U1
Lt_1=copy.deepcopy(T1)
Lt_2=copy.deepcopy(T1)
Lt_1.partial(0,0)
Lt_2.partial(0,1)
Lt_2.diff()
Lt_1-Lt_2
Ls_1=copy.deepcopy(T1)
Ls_2=copy.deepcopy(T1)
Ls_1.partial(1,0)
Ls_2.partial(1,1)
Ls_2.diff()
Ls_1-Ls_2
Lt_1.disp2()
Ls_1.disp2()
T1.disp2()
L1 = copy.deepcopy(Ls_1)
L2 = copy.deepcopy(Lt_1)
L1_1,L1_2 = L1.separete()
L2_1,L2_2 = L2.separete()
'''
Lt_1 :Bewegungsgleichung für θ1
Ls_2 :Bewegungsgleichung für θ2
'''
'''
Führen Sie eine numerische Lösung nach der Euler-Methode durch
Anfangswerteinstellung
Lösung
Anzeige
Müssen 3 Prozesse implementieren
'''
DT = 0.001 #Tageszeiten
x=[math.pi,2,math.pi,2] #Ursprünglicher Wert
saveX = x
T=10
N=int(T/DT) #Berechnungszeit/DT
for i in range(N):
b1 = L1.subst(x) * -1
b2 = L2.subst(x) * -1
a11 = L1_1.subst(x)
a12 = L1_2.subst(x)
a21 = L2_1.subst(x)
a22 = L2_2.subst(x)
hoge = a11*a22-a12*a21
if hoge == 0:
print "hoge = 0 break"
break
alpha1 = (a22*b1-a12*b2)/hoge #A^-1 b
alpha2 = (a11*b2-a21*b1)/hoge
x2 =np.array( x ) + np.array( [ x[1] ,alpha1 ,x[3],alpha2 ] ) * DT
x = x2
saveX = np.append(saveX,x, axis=0)
saveX = saveX.reshape([4,N+1],order ='F')
t=[float(i)*DT for i in range(N+1)]
t= np.asarray(t)
'''
#Punktverlauf
position_x1=r1*np.sin(saveX[0])
position_y1=-1*r1*np.cos(saveX[0])
position_x2=R1*np.sin(saveX[0])+r2*np.sin(saveX[2])
position_y2=-1*R1*np.cos(saveX[0]) - r2*np.cos(saveX[2])
plt.figure(figsize=(28, 20)) # (width, height)
plt.rcParams["font.size"] = 25
plt.plot(position_x1,position_y1,label ="1")
plt.plot(position_x2,position_y2,label="2")
plt.legend(fontsize = 25)
plt.savefig('figure0.png')
'''
'''
#Grafik
z=np.sin(saveX[0])*R1+np.sin(saveX[2])*r2
plt.figure(figsize=(28, 20)) # (width, height)
plt.rcParams["font.size"] = 25
plt.plot(t,z,label ="z")
plt.plot(t,saveX[0],label="x0")
plt.plot(t,saveX[1],label="dt x0")
plt.plot(t,saveX[2],label="x1")
plt.plot(t,saveX[3],label="dt x1")
plt.legend(fontsize = 25)
plt.savefig('figure0.png')
'''
#Animation
print "start animation"
ims = []
fig,ax = plt.subplots()
ax.set_xlim([-2, 2])
ax.set_ylim([-2, 2])
Step = 100 #Anzahl der zu überspringenden Frames
for i in range(int(N/Step)):
position_x1=r1*np.sin(saveX[0][i*Step])
position_y1=-1*r1*np.cos(saveX[0][i*Step])
position_x2=R1*np.sin(saveX[0][i*Step])+r2*np.sin(saveX[2][i*Step])
position_y2=-1*R1*np.cos(saveX[0][i*Step]) - r2*np.cos(saveX[2][i*Step])
im = plt.plot([0,position_x1,position_x2],[0,position_y1,position_y2],color='blue') #Zeichnen Sie Zufallszahlen
ims.append(im) #Fügen Sie dem Array ims ein Diagramm hinzu
#Zeigen Sie alle 100 ms 10 Diagramme an
ani = animation.ArtistAnimation(fig, ims, interval=DT*Step)
ani.save("output2.gif", writer='pillow')