Zweck und Übersicht

Da die Ableitung der kinetischen Gleichung durch die Euler-Lagrange-Gleichung Energie verbraucht, ist sie durch kleine Berechnungsfehler gekennzeichnet. Dies ist eine bequeme Methode, da sie auch auf die komplizierte Bewegungsgleichung mit drei Dimensionen leicht anzuwenden ist. Es ist jedoch schwierig, die partielle Differenzierung nach Ableitung des Lagrange zu berechnen und die simultanen Gleichungen zu lösen. Der Inhalt, der dieses Mal veröffentlicht werden soll, ist, dass nach der Einführung von Lagrange eine teilweise Differenzierung und eine zeitliche Differenzierung durchgeführt werden und die Simulation unter Verwendung der Euler-Methode durchgeführt wird. Noch nicht kompatibel mit Zeichenausdrücken. Außerdem werden nur sin, cos, d / dt ^ 2 / dt ^ 2-Polynome unterstützt.

Für diejenigen, die Matlab und Scilab verwenden können, wird die Zurücksetzung des Browsers empfohlen, da dies völlig unabhängig ist.

3D Y-förmiges Pendel https://www.youtube.com/watch?v=eeizdZscRjU https://www.youtube.com/watch?v=llBik9FsEDo Das Ergebnis der Simulation Ich bin sehr zufrieden mit der Kurve, die zum Experiment passt. Wenn Sie ein Experiment anstelle einer Computersimulation durchführen, funktioniert es gut, wenn Sie das zum Waschen verwendete Pulver verwenden.

Wenn es jedoch um 3D geht, wird es in Bezug auf den Euler-Winkel zu einem besonderen Zustand, und es ist schwierig, den berechneten Wert explosionsartig zu erhöhen oder bei nan anzuhalten und nicht zu arbeiten. Im Moment verarbeite ich es als Moores allgemeine inverse Matrix, aber ich weiß nicht, was ich tun soll, weil ich nicht genug gelernt habe. Jemand sagt es mir Im Moment verwenden wir die Euler-Methode für die Simulation, aber es ist möglicherweise besser, die Lungelkutta-Methode zu verwenden. Ich bin jedoch der Meinung, dass sich die numerische Explosion vom Problem der Berechnungsgenauigkeit unterscheidet.

2D Doppelpendel https://www.youtube.com/watch?v=25feOUNQB2Y Das Ergebnis der Simulation (Jif-Animation, also klicken, um zu sehen) Im Moment sieht es so aus.

Euler-Lagrange-Gleichung

Lagrange $ L $

L=T-U\\
T:Physikalische Energie\\
U:Positionieren Sie die Energie

Auf der anderen Seite

\frac {\partial L}{\partial x} - \frac{d}{dt} \frac{\partial L}{\partial \dot{x}} = 0

Ist die Bewegungsgleichung. x ist die verallgemeinerte Position, wie z. B. Positionskoordinaten und Drehwinkel.

RAM stecken

Ich laufe mit Python2 in der Anaconda-Umgebung. 1 Wie oben erwähnt, ist die Schwierigkeit der Euler-Lagrange-Gleichung

Da es viele teilweise Unterscheidungen gibt, treten wahrscheinlich Berechnungsfehler auf.
Es ist schwierig, nach teilweiser Differenzierung dieselben Begriffe zusammenzustellen
Wenn es um multivariate allgemeine Koordinaten geht, sind die allgemeinen Koordinaten miteinander verflochten und es ist schwierig, sie in die Form $ d ^ 2 x_1 / dt ^ 2 = $ und $ d ^ 2 x_2 / dt ^ 2 = $ (für die Simulation) zu unterteilen. Kann gelöst werden, ohne sich zu teilen

Wenn Sie das Problem tatsächlich lösen, haben Sie es wahrscheinlich satt.

2 Diejenigen, die Matlab und Scilab verwenden können, werden keine Probleme haben. Wenn Sie die Differentialgleichung nach Ableitung der kinetischen Gleichung lösen, gibt es immer eine Lösung wie ode45, daher denke ich, dass dies getan werden sollte. Ich konnte jedoch keinen Weg finden, es in $ d ^ 2 x_1 / dt ^ 2 = $ und $ d ^ 2 x_2 / dt ^ 2 = $ umzuwandeln, also entschied ich mich, es mit Gewalt berechnen zu lassen.

3 Code

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):
        self.L = np.append(self.L,-1*other.L,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:
                        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 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


self.L =Polygoninformationen
self.n =Anzahl der Variablen
self.L = [ [1, 2 ,3 , 4, 5, 6]
           [7, 8, 9, 10,11,12]]

Im Falle von

L = x^1*(\frac{d}{dt} x) ^2*(\frac{d^2}{dt^2} x) ^3  sin(x)^4 cos(x)^5 *6\\
+\\
x^7*(\frac{d}{dt} x) ^8*(\frac{d^2}{dt^2} x) ^9  sin(x)^{10} cos(x)^{11} *12\\

Repräsentiert.

Die Funktion ist

L1 + L2: L1 = L1 + L2 Es ist etwas verwirrend L1-L2 : L1 = L1-L2 *: Nicht implementiert

L1.partial (i, 0): Partielle Differenzierung mit $ x_i $ L1.partial (i, 1): Partielle Differenzierung mit $ \ frac {d} {dt} x_i $

L1.diff (): Zeitdifferenzierung

L1.disp (): Polypoly anzeigen L1.disp2 (): Leicht verständliche Anzeige

Ist implementiert.

Y-förmige Pendelbewegungsgleichung

Als nächstes werde ich erklären, wie man ein Y-förmiges Pendel simuliert. Das Y-förmige Pendel bewegt den V-förmigen Winkel nicht. Mit anderen Worten, es handelt sich um eine Stiftverbindung, die nur eine Drehung in eine Richtung zulässt (sie kann sich nur in der Z-X-Ebene bewegen). Ein Y-förmiges Pendelmodell kann erstellt werden, indem ein Gewicht mit einem Stiftgelenk und einem Kugelgelenk an der Decke befestigt wird (wobei drei Umdrehungen möglich sind). Stellen Sie auf diese Weise den Drehwinkel $ \ theta 123 $ ein.

\vec{R_{m1}} = 
\left(
    \begin{array}{ccc}
      R_1*cos( \theta 1)\\
0\\
-R_1*sin(\theta2)
    \end{array}
  \right)
,
\vec{R_{m2}} = \vec{R_{m1}} +
\left(
    \begin{array}{ccc}
     R_3*sin(\theta3) * cos(\theta2)\\
R_3*sin(\theta3) * sin(\theta2)\\
-R_3*cos(\theta3)
    \end{array}
  \right)

Es wird sein. Dies ist zeitdifferenziert, der Geschwindigkeitsvektor wird berechnet und die kinetische Energie wird aus der Größe des Geschwindigkeitsvektors berechnet.

T=\frac{1}{2} m_1 (R_1\dot{\theta})^2+\frac{1}{2} m_2 (R_1^2 \dot{\theta}^2 + R_3^2 sin(\theta 3)^2 \dot{\theta 2}^2 + R_3^2 \dot{\theta 3}^2
\\
-2*cos(\theta 1) sin(\theta 2) R_1 R_3 sin(\theta 3) \dot{\theta 1} \dot{\theta 2}\\
+2*(cos(\theta 1) cos(\theta 2) cos(\theta 3) + sin(\theta 1) sin(\theta 2)     )R_1 R_3 \dot{\theta 1} \dot{\theta 3}

Positionieren Sie die Energie

U=m_1 g R_1 (1- cos(\theta 1)) + m_2 g (R_1 (1- cos(\theta 1)) + R_3(1-cos(\theta 3))

Wird sein.

Wenn Sie dies verwenden, um Lagrange zu berechnen

  L=   m_1*r_1^2 /2 + m_2*r_1^2/2 \dot{ x_1}^{2.0}               +  
  m_2 /2 *r_3^2      \dot{ x_2}^{2.0}      sin( x_3)^{2.0}     +  
  m_2 /2 *r_3^2           \dot{ x_3}^{2.0}     +  
  m_2 /2 *(-2*r_1*r_3) \dot{ x_1}^{1.0} cos( x_1)^{1.0}      \dot{ x_2}^{1.0} sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  m_2 /2 *(2*r_1*r_3)  \dot{ x_1}^{1.0} cos( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
  m_2 /2 *(2*r_1*r_3)  \dot{ x_1}^{1.0} sin( x_1)^{1.0}           \dot{ x_3}^{1.0} sin( x_3)^{1.0}     +  
  -1*(g*(m_1*r1+m_2*(r_1+r_3)))               +  
  -1*(-1*g*(m_1*r1+m_2*r_1)) cos( x_1)^{1.0}               +  
  -1*(-1*g*m_2*r_3)           cos( x_3)^{1.0}

Es wird sein. θ und x werden vertauscht und angezeigt

x_1 = \theta 1 
x_2 = \theta 2
x_3 = \theta 3

Bei der Ableitung der Bewegungsgleichung

   (m_2 /2 *(2*r_1*r_3) *1.0)  + -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0}           \dot{ x_3}^{1.0} sin( x_3)^{1.0}     +  
   (m_2 /2 *(-2*r_1*r_3)*(-1)*1.0)  + -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *(-1)*1.0) ) \dot{ x_1}^{1.0} sin( x_1)^{1.0}      \dot{ x_2}^{1.0} sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
   (m_2 /2 *(2*r_1*r_3) *(-1)*1.0)  + -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *(-1)*1.0) ) \dot{ x_1}^{1.0} sin( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
   (-1*(-1*g*(m_1*r1+m_2*r_1))*(-1)*1.0)  sin( x_1)^{1.0}               +  
  0 + 0 + 0               +  
  0 + 0               +  
  0               +  
  -1*( ( (m_1*r_1^2 /2 + m_2*r_1^2/2*2.0) *1.0) ) \ddot { x_1}^{1.0}               +  
  -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) cos( x_1)^{1.0}      \dot{ x_2}^{2.0} cos( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *(-1)*1.0)  +  ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) cos( x_1)^{1.0}      \dot{ x_2}^{1.0} sin( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) cos( x_1)^{1.0}      \ddot { x_2}^{1.0} sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) sin( x_1)^{1.0}           \dot{ x_3}^{2.0} cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *(-1)*1.0) ) cos( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{2.0} sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) cos( x_1)^{1.0}      cos( x_2)^{1.0}      \ddot { x_3}^{1.0} cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) sin( x_1)^{1.0}           \ddot { x_3}^{1.0} sin( x_3)^{1.0}     =0\\

  (m_2 /2 *(-2*r_1*r_3)*1.0)  + -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0}      \dot{ x_2}^{1.0} cos( x_2)^{1.0}      sin( x_3)^{1.0}     +  
   (m_2 /2 *(2*r_1*r_3) *(-1)*1.0)  + -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0}      sin( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *(-1)*1.0) ) \dot{ x_1}^{2.0} sin( x_1)^{1.0}      sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(-2*r_1*r_3)*1.0) *1.0) ) \ddot { x_1}^{1.0} cos( x_1)^{1.0}      sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  0 + 0               +  
  -1*( ( (m_2 /2 *r_3^2*2.0) *1.0) )      \ddot { x_2}^{1.0}      sin( x_3)^{2.0}     +  
  -1*( ( (m_2 /2 *r_3^2*2.0) *2.0) )      \dot{ x_2}^{1.0}      \dot{ x_3}^{1.0} sin( x_3)^{1.0} cos( x_3)^{1.0}   =0\\


   (m_2 /2 *r_3^2*2.0)       \dot{ x_2}^{2.0}      sin( x_3)^{1.0} cos( x_3)^{1.0}     +  
   (m_2 /2 *(-2*r_1*r_3)*1.0)  + -1*( ( (m_2 /2 *(2*r_1*r_3) *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}      cos( x_3)^{1.0}     +  
   (m_2 /2 *(2*r_1*r_3) *1.0)  + -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) \dot{ x_1}^{1.0} sin( x_1)^{1.0}           \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
   (m_2 /2 *(2*r_1*r_3) *(-1)*1.0)  + -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *(-1)*1.0) ) \dot{ x_1}^{1.0} cos( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{1.0} sin( x_3)^{1.0}     +  
   (-1*(-1*g*m_2*r_3)*(-1)*1.0)            sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) \dot{ x_1}^{2.0} cos( x_1)^{1.0}           sin( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *(-1)*1.0) ) \dot{ x_1}^{2.0} sin( x_1)^{1.0}      cos( x_2)^{1.0}      cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) \ddot { x_1}^{1.0} cos( x_1)^{1.0}      cos( x_2)^{1.0}      cos( x_3)^{1.0}     +  
  -1*( ( (m_2 /2 *(2*r_1*r_3) *1.0) *1.0) ) \ddot { x_1}^{1.0} sin( x_1)^{1.0}           sin( x_3)^{1.0}     +  
  0 + 0 + 0               +  
  0 + 0               +  
  0               +  
  -1*( ( (m_2 /2 *r_3^2*2.0) *1.0) )           \ddot { x_3}^{1.0}   =0

Es gibt drei Bewegungsgleichungen.

m1=1.0
m2=10
r1=1.0
r3=1.0
g=9.81

Dann

5.5 \dot{ x_1}^{2.0}               +  
  5.0      \dot{ x_2}^{2.0}      sin( x_3)^{2.0}     +  
  5.0           \dot{ x_3}^{2.0}     +  
  -10.0 \dot{ x_1}^{1.0} cos( x_1)^{1.0}      \dot{ x_2}^{1.0} sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  10.0 \dot{ x_1}^{1.0} cos( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
  10.0 \dot{ x_1}^{1.0} sin( x_1)^{1.0}           \dot{ x_3}^{1.0} sin( x_3)^{1.0}     +  
  206.01000000000002               +  
  -107.91000000000001 cos( x_1)^{1.0}               +  
  -98.10000000000001           cos( x_3)^{1.0}

Bei der Ableitung der Bewegungsgleichung

  -107.91000000000001 sin( x_1)^{1.0}               +  
  -11.0 \ddot { x_1}^{1.0}               +  
  10.0 cos( x_1)^{1.0}      \dot{ x_2}^{2.0} cos( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  20.0 cos( x_1)^{1.0}      \dot{ x_2}^{1.0} sin( x_2)^{1.0}      \dot{ x_3}^{1.0} cos( x_3)^{1.0}     +  
  10.0 cos( x_1)^{1.0}      \ddot { x_2}^{1.0} sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  -10.0 sin( x_1)^{1.0}           \dot{ x_3}^{2.0} cos( x_3)^{1.0}     +  
  10.0 cos( x_1)^{1.0}      cos( x_2)^{1.0}      \dot{ x_3}^{2.0} sin( x_3)^{1.0}     +  
  -10.0 cos( x_1)^{1.0}      cos( x_2)^{1.0}      \ddot { x_3}^{1.0} cos( x_3)^{1.0}     +  
  -10.0 sin( x_1)^{1.0}           \ddot { x_3}^{1.0} sin( x_3)^{1.0}   =0\\

  -10.0 \dot{ x_1}^{2.0} sin( x_1)^{1.0}      sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  10.0 \ddot { x_1}^{1.0} cos( x_1)^{1.0}      sin( x_2)^{1.0}      sin( x_3)^{1.0}     +  
  -10.0      \ddot { x_2}^{1.0}      sin( x_3)^{2.0}     +  
  -20.0      \dot{ x_2}^{1.0}      \dot{ x_3}^{1.0} sin( x_3)^{1.0} cos( x_3)^{1.0}     =0\\


10.0      \dot{ x_2}^{2.0}      sin( x_3)^{1.0} cos( x_3)^{1.0}     +  
  -98.10000000000001           sin( x_3)^{1.0}     +  
  -10.0 \dot{ x_1}^{2.0} cos( x_1)^{1.0}           sin( x_3)^{1.0}     +  
  10.0 \dot{ x_1}^{2.0} sin( x_1)^{1.0}      cos( x_2)^{1.0}      cos( x_3)^{1.0}     +  
  -10.0 \ddot { x_1}^{1.0} cos( x_1)^{1.0}      cos( x_2)^{1.0}      cos( x_3)^{1.0}     +  
  -10.0 \ddot { x_1}^{1.0} sin( x_1)^{1.0}           sin( x_3)^{1.0}     +  
  -10.0           \ddot { x_3}^{1.0}    = 0

Es gibt drei Bewegungsgleichungen. Simulieren wir mit dieser Gleichung.

Simulationscode (Euler-Methode)

Position von Geschwindigkeit $ r_ {n + 1} = r_ {n} + v {n} * DT $ Geschwindigkeit aus Beschleunigung $ v_ {n + 1} = v_ {n} + a {n} * DT $ Simulieren Sie durch Aktualisierung.

# -*- 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
'''
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 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 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 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(a1,(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
        
            

#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=10
r1=1.0
r3=1.0
g=9.81

T1=np.array([[0,2,0,0,0,m1*r1**2 /2 + m2*r1**2/2,
             0,0,0,0,0,1,
             0,0,0,0,0,1]])
T2=np.array([[0,0,0,0,0,m2 /2 *r3**2 ,
             0,2,0,0,0,1,
             0,0,0,2,0,1]])
T3=np.array([[0,0,0,0,0,m2 /2 *r3**2 ,
             0,0,0,0,0,1,
             0,2,0,0,0,1]])
T4=np.array([[0,1,0,0,1,m2 /2 *(-2*r1*r3),
             0,1,0,1,0,1,
             0,0,0,1,0,1]])
T5=np.array([[0,1,0,0,1,m2 /2 *(2*r1*r3) ,
             0,0,0,0,1,1,
             0,1,0,0,1,1]])
T6=np.array([[0,1,0,1,0,m2 /2 *(2*r1*r3) ,
             0,0,0,0,0,1,
             0,1,0,1,0,1]])


print "T1"
T=Formula(T1,3)
T.disp2()

print "T2"
R=Formula(T2,3)
R.disp2()
T+R

print "T3"
R=Formula(T3,3)
R.disp2()
T+R

print "T4"
R=Formula(T4,3)
R.disp2()
T+R

print "T5"
R=Formula(T5,3)
R.disp2()
T+R

print "T6"
R=Formula(T6,3)
R.disp2()
T+R

print "T"
T.disp2()

U1=np.array([[0,0,0,0,0,g*(m1*r1+m2*(r1+r3)),
              0,0,0,0,0,1,
              0,0,0,0,0,1]])
U2=np.array([[0,0,0,0,1,-1*g*(m1*r1+m2*r1),
              0,0,0,0,0,1,
              0,0,0,0,0,1]])
U3=np.array([[0,0,0,0,0,-1*g*m2*r3,
              0,0,0,0,0,1,
              0,0,0,0,1,1]])
U = Formula(U1,3)
hoge = Formula(U2,3)
U+hoge
hoge = Formula(U3,3)
U+hoge
print "U"
U.disp2()
T1=copy.deepcopy(T)
T1+U
T-U
L=copy.deepcopy(T)

#Konstruktion der Lagrange-Gleichung
Lt_1=copy.deepcopy(T)
Lt_2=copy.deepcopy(T)
Lt_1.partial(0,0)
Lt_2.partial(0,1)
Lt_2.diff()
Lt_1-Lt_2

Ls_1=copy.deepcopy(T)
Ls_2=copy.deepcopy(T)
Ls_1.partial(1,0)
Ls_2.partial(1,1)
Ls_1.disp2()
print "Ls_2"
Ls_2.disp2()
Ls_2.diff()
Ls_2.disp2()
Ls_1-Ls_2

Lr_1=copy.deepcopy(T)
Lr_2=copy.deepcopy(T)
Lr_1.partial(2,0)
Lr_2.partial(2,1)
Lr_2.diff()
Lr_1-Lr_2

L1t=copy.deepcopy(Lt_1)
L1t_1,L1t_2,L1t_3 = L1t.separete3()
L2t=copy.deepcopy(Ls_1)
L2t_1,L2t_2,L2t_3 = L2t.separete3()
L3t=copy.deepcopy(Lr_1)
L3t_1,L3t_2,L3t_3 = L3t.separete3()



'''
Führen Sie eine numerische Lösung nach der Euler-Methode durch
Anfangswerteinstellung
Lösung
Anzeige
Müssen 3 Prozesse implementieren
'''


DT = 0.0001              #Tageszeiten
x=np.array([0.1,0,0.1,0,0.1,0])            #Ursprünglicher Wert
saveX = x
T=10
N=int(T/DT)                 #Berechnungszeit/DT

saveEnergy = np.asarray([T1.subst(x)])

for i in range(N):
    b1 = L1t.subst(x) * -1
    b2 = L2t.subst(x) * -1
    b3 = L3t.subst(x) * -1
    a11 = L1t_1.subst(x)
    a12 = L1t_2.subst(x)
    a13 = L1t_3.subst(x)
    a21 = L2t_1.subst(x)
    a22 = L2t_2.subst(x)
    a23 = L2t_3.subst(x)
    a31 = L3t_1.subst(x)
    a32 = L3t_2.subst(x)
    a33 = L3t_3.subst(x)
    B=np.array([b1,b2,b3])
    B=B.reshape([3,1])
    A=np.array([a11,a12,a13,a21,a22,a23,a31,a32,a33])
    A=A.reshape([3,3])
    
    hoge = np.linalg.det(A)
    if hoge == 0:
        #N=i
        #print "hoge = 0 break"
        Aneg = np.linalg.pinv(A)
    
    else:
        Aneg = np.linalg.inv(A)
    
    '''
    if hoge < 0.001:
        N=i
        print "hoge > 10000 break"
        break
    '''
    
    
    K = np.dot(Aneg,B)
    
    x2 =np.array( x ) + np.array( [ x[1] ,K[0][0] ,x[3], K[1][0], x[5] ,K[2][0]  ] ) * DT
    
    x = x2
    saveX = np.append(saveX,x, axis=0)
    
    Energy = T1.subst(x)
    saveEnergy = np.append(saveEnergy,[Energy], axis=0)


saveX = saveX.reshape([6,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]) +r3*np.sin(saveX[4])*np.cos(saveX[2])
position_y1=r3*np.sin(saveX[4])*np.sin(saveX[2])
plt.figure(figsize=(28, 20))  # (width, height)
plt.rcParams["font.size"] = 25
plt.plot(position_x1,position_y1,label ="1")

plt.legend(fontsize = 25)

plt.savefig('Position_test.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.xlim(0,T)
#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.plot(t,saveX[4],label="x2")
plt.plot(t,saveX[5],label="dt x2")
plt.legend(fontsize = 25)

plt.savefig('Pendulum_test.png')

plt.figure(figsize=(28, 20))  # (width, height)
plt.rcParams["font.size"] = 25
plt.xlim(0,T)
#plt.plot(t,z,label ="z")

plt.plot(t,saveEnergy,label="Energy")
plt.legend(fontsize = 25)

plt.savefig('Energy_test.png')

Auf diese Weise wird die Lissajous-Kurve angewendet.