[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

define.jpg

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\\

output2.gif

Ursprünglicher Wert:

\theta_1 = 180deg\\
\theta_2 = 180deg\\
\dot{\theta_1 }= 0.1rad/s\\
\dot{\theta_2} = 0.1rad/s

output2.gif 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

output2.gif 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')

Recommended Posts

Doppelpendelsimulation
Rolltreppensimulation
Doppelte Pendelbewegungsgleichung in Python