Laissez Python calculer l'équation de mouvement d'Euler Lagrange

Objectif et aperçu

Comme la dérivation de l'équation cinétique par l'équation d'Euler-Lagrange utilise de l'énergie, elle est caractérisée par de petites erreurs de calcul. C'est une méthode pratique car elle est facile à appliquer même à l'équation de mouvement compliquée de 3 dimensions. Cependant, il est difficile de calculer la différenciation partielle après avoir dérivé le lagrangien et de résoudre les équations simultanées. Le contenu à publier cette fois est qu'après l'introduction du lagrangien, une différenciation partielle et une différenciation temporelle sont effectuées, et la simulation est effectuée à l'aide de la méthode d'Euler. Pas encore compatible avec les expressions de caractères. De plus, il ne prend en charge que les polynômes sin, cos, d / dt ^ 2 / dt ^ 2.

Pour ceux qui peuvent utiliser matlab et scilab, le retour du navigateur est recommandé car cela n'a rien à voir avec cela.

Pendule 3D en forme de Y https://www.youtube.com/watch?v=eeizdZscRjU https://www.youtube.com/watch?v=llBik9FsEDo Le résultat de la simulation Je suis assez satisfait de la courbe qui correspond à l'expérience. Si vous faites une expérience au lieu d'une simulation informatique, cela fonctionnera bien si vous utilisez la poudre utilisée pour le lavage.

Cependant, quand il s'agit de 3D, cela devient un état particulier en termes d'angle d'Euler, et il est difficile pour la valeur calculée d'augmenter de manière explosive ou de s'arrêter à nan et de ne pas fonctionner. Pour le moment, je la traite comme la matrice inverse générale de Moore, mais je ne sais pas quoi faire car je n'ai pas suffisamment étudié. Quelqu'un me dit Pour l'instant, nous utilisons la méthode Euler pour la simulation, mais il peut être préférable d'utiliser la méthode Lungelkutta. Cependant, j'estime que l'explosion numérique est différente du problème de la précision des calculs.

Double pendule 2D https://www.youtube.com/watch?v=25feOUNQB2Y Le résultat de la simulation (animation Jif, donc cliquez pour voir) Pour le moment, ça ressemble à ça.

Équation d'Euler-Lagrange

Lagrangien $ L $

L=T-U\\
T:Énergie physique\\
U:Énergie de position

D'autre part

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

Est l'équation du mouvement. x est la position généralisée, comme les coordonnées de position et l'angle de rotation.

Branchez le bélier

Je cours avec python2 dans l'environnement anaconda. 1 Comme mentionné ci-dessus, la difficulté de l'équation d'Euler-Lagrange est

• Comme il existe de nombreuses différenciations partielles, des erreurs de calcul sont susceptibles de se produire.
• Il est difficile de rassembler les mêmes termes après une différenciation partielle
• Lorsqu'il s'agit de coordonnées générales multivariées, les coordonnées générales sont entrelacées et il est difficile de les diviser sous la forme $ d ^ 2 x_1 / dt ^ 2 = $ et $ d ^ 2 x_2 / dt ^ 2 = $ (pour la simulation) Peut être résolu sans se diviser en

Lorsque vous résolvez réellement le problème, vous en avez probablement assez.

2 Ceux qui peuvent utiliser matlab et scilab ne seront pas en difficulté. De plus, lors de la résolution de l'équation différentielle après avoir dérivé l'équation de mouvement, il y a toujours une solution telle que ode45, donc je pense que cela devrait être fait. Cependant, je n'ai pas trouvé de moyen de le convertir en $ d ^ 2 x_1 / dt ^ 2 = $ et $ d ^ 2 x_2 / dt ^ 2 = $, j'ai donc décidé de le faire calculer de force.

3 code

class Formula :
    def __init__(self,L,num):
        self.L = L
        self.n = num                      #Nombre de types de variables
        
    
    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):               #Résumer les coefficients du même terme
        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):    #Effacez le terme avec un coefficient de 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):  #Différencier partiellement les termes avec le genre
        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
            '''
Tel quel
Diminuer la commande de un, si 0, ne rien faire
Multipliez l'ordre d'origine par le coefficient
            '''
            #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 si non inclus
            
        if kind == 1:   
            '''
Tel quel
Diminuer la commande de un, si 0, ne rien faire
Multipliez l'ordre d'origine par le coefficient
            '''
            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 si non inclus
                    
        self.together()
        self.erase()
        
    def diff(self):  #Différenciation temporelle
        '''
Après différenciation partielle avec 4 motifs
Augmenter l'ordre de θ
        '''
        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):  #Affichage polymère
        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,Mettez dans l'ordre de x
        vari =np.array(x).shape
        if vari[0] != self.n*2:
            print "cannot subst"
        else:
            '''
Générer un vecteur de ligne
Faites L avec un point et ajoutez tous les éléments du vecteur résultant
            '''
            
            expr,vari = self.L.shape
            sum1 = 0
            hoge=0
            for j in range(expr):
                hoge=self.L[j][5]           #coefficient
                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 =Informations sur les polygones
self.n =Nombre de variables
self.L = [ [1, 2 ,3 , 4, 5, 6]
           [7, 8, 9, 10,11,12]]

dans le cas de

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ésente.

la fonction est

L1 + L2: L1 = L1 + L2 C'est un peu déroutant L1-L2 : L1 = L1-L2 *: Pas mis en œuvre

L1.partial (i, 0): Différenciation partielle avec $ x_i $ L1.partial (i, 1): Différenciation partielle avec $ \ frac {d} {dt} x_i $

L1.diff (): Différenciation temporelle

L1.disp (): Afficher le polypoly L1.disp2 (): affichage facile à comprendre

Est implémenté.

Équation du mouvement du pendule en forme de Y

Ensuite, j'expliquerai comment simuler un pendule en forme de Y. Le pendule en forme de Y ne déplace pas l'angle en forme de V. En d'autres termes, il s'agit d'une jonction de broche qui ne permet qu'une rotation à sens unique (elle ne peut se déplacer que dans le plan Z-X). Un modèle de pendule en forme de Y peut être créé en attachant un poids au plafond avec une articulation à goupille et une rotule (permettant trois rotations). De cette façon, définissez l'angle de rotation $ \ theta 123 $.

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

Ce sera. Ceci est différencié dans le temps, le vecteur vitesse est calculé et l'énergie cinétique est calculée à partir de l'amplitude du vecteur vitesse.

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}  

Énergie de position

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

Sera.

Si vous utilisez ceci pour calculer Lagrangean

  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}         

Ce sera. θ et x sont interchangés et affichés

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

Lors de la dérivation de l'équation du mouvement

   (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

Il existe trois équations de mouvement.

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

Puis

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}   

Lors de la dérivation de l'équation du mouvement

  -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

Il existe trois équations de mouvement. Simulons en utilisant cette équation.

Code de simulation (méthode Euler)

Position à partir de la vitesse $ r_ {n + 1} = r_ {n} + v {n} * DT $ Vitesse d'accélération $ v_ {n + 1} = v_ {n} + a {n} * DT $ Simulez en mettant à jour.

# -*- coding: utf-8 -*-
"""
Created on Sat Mar 21 14:11:03 2020

@author: kisim
"""
import numpy as np		#Bibliothèque Numpy
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
'''
Programme de traitement de formule
Colonne: Polygone
Ligne: θ θ ・ θ ・ ・ sin θ cos θ coefficient
6 pièces x nombre de variables


'''
class Formula :
    def __init__(self,L,num):
        self.L = L
        self.n = num                      #Nombre de types de variables
        
    
    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):               #Résumer les coefficients du même terme
        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):    #Effacez le terme avec un coefficient de 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):  #Différencier partiellement les termes avec le genre
        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
            '''
Tel quel
Diminuer la commande de un, si 0, ne rien faire
Multipliez l'ordre d'origine par le coefficient
            '''
            #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 si non inclus
            
        if kind == 1:   
            '''
Tel quel
Diminuer la commande de un, si 0, ne rien faire
Multipliez l'ordre d'origine par le coefficient
            '''
            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 si non inclus
                    
        self.together()
        self.erase()
        
    def diff(self):  #Différenciation temporelle
        '''
Après différenciation partielle avec 4 motifs
Augmenter l'ordre de θ
        '''
        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):  #Affichage polymère
        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):  #Affichage polymère
        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,Mettez dans l'ordre de x
        vari =np.array(x).shape
        if vari[0] != self.n*2:
            print "cannot subst"
        else:
            '''
Générer un vecteur de ligne
Faites L avec un point et ajoutez tous les éléments du vecteur résultant
            '''
            
            expr,vari = self.L.shape
            sum1 = 0
            hoge=0
            for j in range(expr):
                hoge=self.L[j][5]           #coefficient
                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):              #Considérez-le comme 2 variables et divisez-le en 3 parties, 1 et 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):              #Considérez-le comme 3 variables et divisez-le en 4 parties: 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,Considérons le cas où il n'y a que deux expressions de L2 et il y a deux variables.

Comme la formule ≒ 0
Extraire uniquement l'élément de θ ... et faire une nouvelle instanciation
Le supprimé est désigné comme b[L1_b,L2_b]Faire un vecteur
Le terme de θ ...
[L1_1   L1_2
 L2_1   L2_2 ]Considérons la matrice d'instance A

 L1_2 =Élément de θ2 dans la formule de L1
θ ...=　A^-Calculer 1 b
 
Θ ・ θ apparaît également dans la méthode explicite de la méthode d'Euler
 
Faisons cela
 
 L1_b,L1_1,L1_Pensez à comment générer 2
 
'''

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)

#Construction de l'équation de Lagrange
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()



'''
Effectuer une solution numérique par la méthode Euler
Réglage de la valeur initiale
Solution
afficher
Besoin de mettre en œuvre 3 processus
'''


DT = 0.0001              #Heures du jour
x=np.array([0.1,0,0.1,0,0.1,0])            #valeur initiale
saveX = x
T=10
N=int(T/DT)                 #Temps de calcul/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)




#Historique des points
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')




#Graphisme
#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')

Comme ça, la courbe de Lissajous est appliquée.