Simulation de double pendule

Qu'est-ce qu'un double pendule?

https://www.youtube.com/watch?v=25feOUNQB2Y C'est une expérience célèbre en tant que phénomène de chaos. Le phénomène de chaos est un phénomène dans lequel l'état initial a une grande influence sur l'état ultérieur. L'histoire de l'effet papillon, dans lequel les ailes d'un papillon génèrent un typhon en raison du flux d'air, est également évoquée dans l'explication du même phénomène de chaos. (Peut-être pas correct)

Nous dériverons et simulerons l'équation du mouvement en utilisant la méthode Euler-Lagrange créée par ici.

Définition variable

define.jpg

Définissez le problème dans le système de points de qualité / deux dimensions.

La vitesse de m1 et la vitesse de 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}\\


Énergie d'exercice T

T=\frac{1}{2} m_1 |\vec{V_{m_1}}|^2+\frac{1}{2} m_2 |\vec{V_{m_2}}|^2

Position énergie U

U=m_1gr_1(1-cos(\theta_1)) + m_2g(r_1(1-cos(\theta_1)) + r_2(1-cos(\theta_2))

Lorsque vous calculez Lagrangean

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

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

(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

Le code que j'ai utilisé est difficile à voir, je vais donc le joindre à la fin.

Simulation (méthode Euler)

Lorsque simulé par la méthode Euler en utilisant la formule obtenue Réglage

m1=1.0\\
m2=1.0\\
r1=1.0\\
r2=1.0\\
g=9.81\\

valeur initiale:

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

output2.gif

valeur initiale:

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

output2.gif Celui-ci est plus chaotique (sentiment moyen 2 moyen)

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

output2.gif Ça ne peut pas être chaotique sans vitesse

Code de dérivation de l'équation de mouvement

A partir de là, ce sera le programme que vous avez créé. Cela devrait fonctionner avec le copier-coller. Je cours sur anaconda python2. gif Cliquez ici pour animer https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2 Je l'ai mentionné.

Classe d'analyse

Formula2


import numpy as np		#Bibliothèque Numpy
import copy
import math
'''
Programme de traitement de formule
Colonne: Polygone
Ligne: θ θ ・ θ ・ ・ sin θ cos θ coefficient
6 pièces x nombre de variables


Correspond à l'expression de caractère
cependant,

'''
class Formula2 :
    def __init__(self,L,num,coeff):
        self.L = L
        self.n = num                      #Nombre de types de variables
        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):               #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:  #Tous concordent Ajout de coefficients
                        self.L[i][5] += self.L[j][5]
                        self.coeff[i] = self.coeff[i] +" + "+ self.coeff[j]
                        for k in range(vari):  #j ligne supprimée
                            self.L[j][k] = 0
                        self.coeff[j] = "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.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):  #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])
                    hogestr = copy.deepcopy(self.coeff[i])
                    hoge[5] = hoge[5]*hoge[3+moji*6]
                    #Hoge au coefficient[3+moji*6]multiplier
                    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 au coefficient[3+moji*6]multiplier
                    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
            '''
Tel quel
Diminuer la commande de un, si 0, ne rien faire
Multipliez l'ordre d'origine par le coefficient
            sin,Ce n'est pas grave car la partie de traitement cos est en retard et n'est pas dans la plage de traitement
            '''
            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 si non inclus
                    self.coeff[i] = "0"
            
        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.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 si non inclus
                    self.coeff[i] = "0"
                    
        
        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)
        L1_coeff = copy.deepcopy(self.coeff)
        
        for i in range(self.n):
            self.L =copy.deepcopy( L1)   #L1:Avant différenciation L2:Pour enregistrer les résultats
            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):  #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 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):              #Considérez-le comme 2 variables et divisez-le en 3 parties, 1 et 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):              #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([[]])
        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()
        '''

Code d'exécution

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

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,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_Réfléchissez à comment générer 2
'''

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:Lagrangien
Lt_1 ,L1 :Équation de mouvement pour θ1
Ls_2 ,L2 :Equation de mouvement pour θ2
'''

Code de simulation

Classe d'analyse

Formula.py


import numpy as np		#Bibliothèque Numpy
import copy
import math
'''
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 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):               #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 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(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):              #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
        
    if __name__ == '__main__':
        print "This is not for main"

Programme d'exécution

# -*- 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

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,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_Réfléchissez à comment générer 2
'''

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  :Équation de mouvement pour θ1
Ls_2  :Equation de mouvement pour θ2
'''

'''
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.001              #Heures du jour
x=[math.pi,2,math.pi,2]            #valeur initiale
saveX = x
T=10
N=int(T/DT)                 #Temps de calcul/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)

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

'''
#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.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   #Nombre d'images à ignorer

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')             #Graphique des nombres aléatoires
        ims.append(im)                  #Ajouter un graphique au tableau ims

#Afficher 10 graphiques toutes les 100 ms
ani = animation.ArtistAnimation(fig, ims, interval=DT*Step)
ani.save("output2.gif", writer='pillow')

Recommended Posts

Simulation de double pendule
Simulation d'escalator
Equation de mouvement à double pendule en python