2. TP DAE : pendule sur une courbe#
Marc BUFFAT dpt mécanique, Université Lyon 1
vous devez écrire les fonctions dont le nom est fixé, mais dont vous devez spécifier les arguments en fonction de la question posée.
Chaque fonction validée rapporte des points
Attention: executer toutes les cellules depuis le début en utilisant le bouton run
# initialisation
import os,sys
import numpy as np
import matplotlib.pyplot as plt
from validation.validation import check_function,liste_functions,info_etudiant
from IPython.display import Markdown, display, HTML
def printmd(string):
display(Markdown(string))
# test si numero étudiant spécifier
try: NUMERO_ETUDIANT
except NameError: NUMERO_ETUDIANT = None
if type(NUMERO_ETUDIANT) is not int :
printmd("**ERREUR:** numéro d'étudiant non spécifié!!!")
NOM,PRENOM,NUMERO_ETUDIANT=info_etudiant()
# parametres spécifiques
_uid_ = NUMERO_ETUDIANT
_precis_ = 1.0e-5
printmd("**Login étudiant {} {} uid={}**".format(NOM,PRENOM,_uid_))
np.random.seed(_uid_)
_m_=np.round(1.0+2*np.random.rand(),2) #masse suivant la trajectoire
_l_=np.round(1.0+2*np.random.rand(),2) #longueur de la barre
_M_=np.round(3.0+4*np.random.rand(),2) #masse au bout du pendule
_a0_=np.round(2.0/(2+np.random.randint(4)),2)
printmd("**Parametres:** m={} M={} l={} a0={}".format(_m_,_M_,_l_,_a0_))
ERREUR: numéro d’étudiant non spécifié!!!
Login étudiant Marc BUFFAT uid=137764122
Parametres: m=2.7 M=6.55 l=2.81 a0=0.5
# bibliotheque scikit
from scikits.odes.dae import dae
from scikits.odes.ode import ode
2.1. Pendule sur une courbe#
2.1.1. modélisation#
Les coordonnées du système sont \(Q=[x(t),y(t),\theta(t)]\), (coordonnées de m et angle de rotation du pendule). Le lagrangien du système est \(L(Q,\dot{Q},t)\) , et les équations de Lagrange sous contrainte s’écrivent
\[\frac{d}{dt}(\frac{\partial L}{\partial\dot{Q}_{i}})-\frac{\partial L}{\partial Q_{i}}=H_{i}+\lambda G_{i}\]
où le vecteur \(\vec{H}\) représente la contribution du travail \(\vec{H}.\delta\vec{Q}\) de la force de frottement \(\vec{F}\), \(\lambda\) le multiplicateur de Lagrange (i.e. la force de liaison) associé à la contrainte \(f(x,y)=0\), et \(\vec{G}\) la contribution due travail de cette force de liaison \(\lambda\vec{G}.\delta\vec{Q}\). $\( G_i = \frac{\partial f}{\partial Q_i}\)$
La courbe \(f(x,y)=0\) est le chameau \(y=x^{2}+a_0\cos(\omega x)\)
A l’instant initial le pendule se trouve la courbe en \(x_0, y_0=f(x_0)\) et on le lâche sans vitesse initiale avec un angle \(\theta=0\). On cherche à déterminer la position d’équilibre du système en fonction de \(x_{0}\).
2.1.2. Définition des repères et des points#
un référenciel R0 et un point O de référence:
système de deux points:
Q(xp,yp) de masse M glissant sur la courbe auquel est accroché l’autre point P
P distant de l de Q
# sympy avec Lagrange
from sympy import init_printing, simplify, symbols, sin, cos, Matrix
from sympy.physics.mechanics import dynamicsymbols, Point, ReferenceFrame, \
Particle, LagrangesMethod
init_printing()
# parametres et ddl
t, M , m, l, g, a0, omega, K = symbols('t M m l g a0 omega K')
x, y, theta = dynamicsymbols('x y theta')
xp, yp, thetap = dynamicsymbols('x y theta',1)
# définition des points
O = Point('O')
R0 = ReferenceFrame('R_0')
P = Point('P')
P.set_pos(O,x*R0.x + y*R0.y)
Q = Point('Q')
Q.set_pos(P,l*sin(theta)*R0.x - l*cos(theta)*R0.y)
P.set_vel(R0,P.pos_from(O).diff(t,R0))
display("VP=",P.vel(R0),"VQ=",Q.vel(R0))
'VP='
'VQ='
# definition des points matériels
PaP = Particle('Pa_P',P,m)
PaQ = Particle('Pa_Q',Q,M)
PaP.potential_energy = m*g*P.pos_from(O).dot(R0.y)
display(PaP.potential_energy, PaP.kinetic_energy(R0))
PaQ.potential_energy = M*g*Q.pos_from(O).dot(R0.y)
display(PaQ.potential_energy, PaQ.kinetic_energy(R0))
2.1.3. Calcul du lagrangien#
# Lagrangien
L = PaP.kinetic_energy(R0) + PaQ.kinetic_energy(R0) - PaP.potential_energy - PaQ.potential_energy
display("L=",L)
'L='
# contrainte
f = y - (x**2 + a0*cos(omega*x))
display("contrainte=",f)
'contrainte='
# force frottement
FL = [(P, -K*xp*R0.x - K*yp*R0.y)]
display("Force de frottement=",FL)
'Force de frottement='
[(P, - K*Derivative(x(t), t)*R_0.x - K*Derivative(y(t), t)*R_0.y)]
2.2. Bilan des équation#
Énergie cinétique
\[T=\frac{1}{2}m(\dot{x}^{2}+\dot{y}^{2})+\frac{1}{2}M\left((\dot{x}+l\dot{\theta}\cos\theta)^{2}+(\dot{y}+l\dot{\theta}\sin\theta)^{2}\right)\]
Énergie potentielle
\[U=mgy+Mg(y-lcos\theta)\]
Lagrangien
\[L=T-U\]
Force de frottement généralisée
\[\begin{split}H=\left[\begin{array}{c}
-K\dot{x}\\
-K\dot{y}\\
0
\end{array}\right]\end{split}\]
contrainte $\( f(x,y) = y-(x^{2}+a_0\cos\omega x)\)$
gradient de la contrainte $\( G = \left[\begin{array}{c} \frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\\ 0 \end{array}\right]\)$
Ecrire les équations du mouvement, et les transformer en un système d’ordre 1 de la forme:
\[ M \dot{Y} = F(Y,t) \]
avec
\[Y=[x,y,\theta,\dot{x},\dot{y},\dot{\theta},\lambda]\]
on dérivera la contrainte autant de fois que nécessaire pour pouvoir résoudre.
LM = LagrangesMethod(L, [x,y,theta], hol_coneqs=[f], forcelist=FL, frame=R0)
eqs = LM.form_lagranges_equations()
display(eqs)
2.2.1. Mise sous forme matricielle#
\[ A X = B(Y) \mbox{ avec } Y=[x,y,\theta,\dot{x},\dot{y},\dot{\theta}] \mbox{ et } X=[\dot{Y},\lambda]\]
A=LM.mass_matrix_full
display(A)
B=LM.forcing_full
display(B)
2.2.2. linéarisation#
pour des cooordonnées \(q\) et des vitesses généralisées \(u\), dont \(q_i\) et \(u_i\) sont indépendantes et \(r\) un forcage externe
\begin{split}M \begin{bmatrix} \delta \dot{q} \ \delta \dot{u} \ \delta \lambda \end{bmatrix} = A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix}\end{split}
xe=0
ye=xe**2+a0*cos(omega*xe)
print(xe,ye)
P0 = { x:xe, y:ye, theta:0, xp:0, yp:0, thetap:0, xp.diff(t):0,yp.diff(t):0,thetap.diff(t):0}
l0 = LM.solve_multipliers(op_point=P0)
P0.update(l0)
LIN = LM.to_linearizer(q_ind=[x,theta],qd_ind=[xp,thetap], q_dep=[y], qd_dep=[yp])
ML,AL,BL = LIN.linearize(op_point = P0)
display(ML)
display(AL)
display(BL)
0 a0
lambda1 = symbols('lambda1')
X = Matrix([xp,yp,thetap,xp.diff(t),yp.diff(t),thetap.diff(t),lambda1])
Y = Matrix([x,theta,xp,thetap])
ML*X - AL*Y
AL,BL, inp_vec = LM.linearize([x,theta],[xp,thetap],[y],[yp],op_point=P0, A_and_B=True)
display(AL)
display(BL)
Yp = Matrix([xp,thetap,xp.diff(t),thetap.diff(t)])
Yp - AL*Y
# A.N
A1=np.array(AL.subs([(M,_M_),(m,_m_),(l,_l_),(g,10),(K,2),(omega,3.3),(a0,_a0_)])).astype(np.float64)
display(A1)
np.linalg.eig(A1)
array([[ 0. , 0. , 1. , 0. ],
[ 0. , 0. , 0. , 1. ],
[118.02314815, 24.25925926, -0.74074074, 0. ],
[-42.00112034, -12.19190721, 0.2636088 , 0. ]])
EigResult(eigenvalues=array([-1.08479916e+01+0.j , 1.01090216e+01+0.j ,
-8.85400092e-04+1.95704657j, -8.85400092e-04-1.95704657j]), eigenvectors=array([[ 8.67633393e-02+0.j , 9.30902155e-02+0.j ,
-1.01408097e-03+0.08883113j, -1.01408097e-03-0.08883113j],
[-2.99703036e-02+0.j , -3.20133654e-02+0.j ,
-2.01894326e-04-0.44625769j, -2.01894326e-04+0.44625769j],
[-9.41207973e-01+0.j , 9.41051003e-01+0.j ,
-1.73845769e-01-0.00206325j, -1.73845769e-01+0.00206325j],
[ 3.25117601e-01+0.j , -3.23623804e-01+0.j ,
8.73347255e-01+0.j , 8.73347255e-01-0.j ]]))
2.2.3. changement de variables Y#
from sympy import MatrixSymbol,ImmutableMatrix
Y = MatrixSymbol('Y',7,1)
# met derniere colonne de A (xcateur lagrange) dans B
B1 = B - Y[6]*A[:,6]
A1 = A[:,:]
A1[:,6] = [0,0,0,0,0,0,0]
# penalisation de la contrainte
beta = symbols('beta')
B1[6] = B1[6] - beta*f
# substitution
A1=A1.subs([(xp,Y[3]),(yp,Y[4]),(thetap,Y[5]),(x,Y[0]),(y,Y[1]),(theta,Y[2])])
A1 = ImmutableMatrix(A1)
display(A1)
B1=B1.subs([(xp,Y[3]),(yp,Y[4]),(thetap,Y[5]),(x,Y[0]),(y,Y[1]),(theta,Y[2])])
B1 = ImmutableMatrix(B1)
display(B1)
# generation automatique des fonctions
from sympy.utilities.autowrap import autowrap
!rm -rf wrapp*
matA = autowrap(A1,args=(Y,M,m,l,g,K,a0,omega),backend='cython',tempdir='./autowrap')
smbB = autowrap(B1,args=(Y,M,m,l,g,K,a0,omega,beta),backend='cython',tempdir='./autowrap')
2.3. Programmation#
2.3.1. parametres#
les valeurs de m,l,M et a0 sont fixés au début.
on choisira la valeur des parametres de penalisation \(\beta_1\) \(\beta_2\)
on prendra \(K=2.0\) , \(\omega=3.3\), \(g=10.0\)
le seul paramêtre variable est donc la position initiale \(x_0\)
Définir la valeurs de ces parametres et ecrire une fonction résidu qui calcul le résidu du système:
\[ res = M.\dot{Y} - F(Y,t) \]
# parametres
## BEGIN SOLUTION
m=_m_ #masse suivant la trajectoire
l=_l_ #longueur de la barre
M=_M_ #masse au bout du pendule
a0 = _a0_
K=2.0 #coefficient de frottement
g=10.0 #constante de pesenteur
omega=3.3
# penalisation
beta1=0.
beta2=1.e2
nit = 0
## END SOLUTION
def residu(t,Y,dY,res):
'''calcul res=M dY -F(Y,t)'''
global nit
nit += 1
YY = np.zeros((7,1))
YY[:,0] = Y[:]
A = matA(YY,M,m,l,g,K,a0,omega)
B = smbB(YY,M,m,l,g,K,a0,omega,beta2)
res[:] = A.dot(dY) - B[:,0]
return
2.3.2. Vérification#
vérifier que pour les 2 positions d’équilibres le résidu est bien nul. On tracera aussi la courbe f(x) et la position des 2 points d’équilibre stables.
## BEGIN SOLUTION
from scipy.optimize import fsolve
f = lambda x : x**2 + a0*np.cos(omega*x)
df = lambda x: 2*x - a0*omega*np.sin(omega*x)
xe = fsolve(df,1.0)[0]
print("racine {} err={}".format(xe,df(xe)))
X = np.linspace(-2,2,100)
# verification
dY0 = np.zeros(7)
Y0 = np.array([xe,f(xe),0.,0.,0.,0.,(m+M)*g])
res = np.zeros(7)
residu(0,Y0,dY0,res)
print("pour xe={} res={}".format(xe,res))
Y0 = [-xe,f(-xe),0.,0.,0.,0.,(m+M)*g]
res = np.zeros(7)
residu(0,Y0,dY0,res)
print("pour xe={} res={}".format(-xe,res))
# tracer
plt.plot(X,f(X))
plt.plot([xe,-xe],[f(xe),f(-xe)],'o')
## END SOLUTION
racine 0.666798402972126 err=5.773159728050814e-15
pour xe=0.666798402972126 res=[0.00000000e+00 0.00000000e+00 0.00000000e+00 5.34017275e-13
0.00000000e+00 0.00000000e+00 0.00000000e+00]
pour xe=-0.666798402972126 res=[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -5.34017275e-13
0.00000000e+00 0.00000000e+00 0.00000000e+00]
[<matplotlib.lines.Line2D at 0x7fa0fa3a9930>]
2.4. Résolution#
ecrire une fonction solution qui calcule la solution en fonction d’une position initiale x0
def solution(x0):
'''calcul solution avec CI x0. Renvoie T et Y'''
## BEGIN SOLUTION
global nit
# conditions initiales (recalculer par le solveur)
t0=0.0
Y0 = np.array([x0, f(x0), 0.0, 0.0, 0.0, 0.0, 0.0])
dY0= np.array([0., 0., 0., 0., -g , 0., 0.])
res= np.zeros(7)
residu(t0,Y0,dY0,res)
print("CI Y0=",Y0,"\n dY0=",dY0,"\n res=",res)
#choix du solveur
solver = dae('ida', residu, compute_initcond='yp0',
first_step_size=1e-10,
atol=1e-6,rtol=1e-6,
algebraic_vars_idx=[6],
old_api=False, max_steps=5000)
# integration
tfinal=30*np.pi/np.sqrt(g/l) # periode
N=500 # nombre de points
t = np.linspace(0,tfinal,N)
nit = 0
solution = solver.solve(t,Y0,dY0)
print(solution.message)
print("Nbre d'appel a rhs: ",nit)
T = solution.values.t
Y = solution.values.y
return T,Y
## END SOLUTION
2.4.1. Vérification#
en choissisant une valeur de x0 proche de la position d’équilibre, vérifier que le mouvement du pendule est celui attendu.
## BEGIN SOLUTION
x0 = xe + 0.1
T,Y = solution(x0)
#trace de la trajectoire
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
plt.plot(Y[:,0],Y[:,1],label='m')
plt.plot(Y[:,0]+l*np.sin(Y[:,2]),Y[:,1]-l*np.cos(Y[:,2]),lw=2,label='M')
plt.axis('equal')
plt.legend(loc=0)
plt.subplot(1,2,2)
plt.plot(T,Y[:,0],label='x')
plt.plot(T,Y[:,1],label='y')
plt.legend(loc=0)
## END SOLUTION
CI Y0= [0.7667984 0.17848772 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 8156
<matplotlib.legend.Legend at 0x7fa0f8aabbe0>
2.5. Analyse pour x0 grand#
pour une valeur de x0 assez grande vérifier que le pendule peut passer d’une position déquilibre à l’autre.
## BEGIN SOLUTION
x0 = 2.0
T,Y = solution(x0)
#trace de la trajectoire
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
plt.plot(Y[:,0],Y[:,1],label='m')
plt.plot(Y[:,0]+l*np.sin(Y[:,2]),Y[:,1]-l*np.cos(Y[:,2]),lw=2,label='M')
plt.axis('equal')
plt.legend(loc=0)
plt.subplot(1,2,2)
plt.plot(T,Y[:,0],label='x')
plt.plot(T,Y[:,1],label='y')
plt.legend(loc=0)
## END SOLUTION
CI Y0= [2. 4.4751163 0. 0. 0. 0. 0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 48738
<matplotlib.legend.Legend at 0x7fa0f8549270>
2.6. Etude en fonction de x0#
en faisant varier x0 de 0 2.0, déterminer la position finale en fonction de x0
## BEGIN SOLUTION
N = 21
XF=np.zeros(N)
X0=np.linspace(0.01,2.,N)
for i in range(N):
T,Y = solution(X0[i])
if np.abs(Y[-1,0]-xe) < 0.1 :
XF[i]=xe
if np.abs(Y[-1,0]+xe) < 0.1 :
XF[i]=-xe
## END SOLUTION
CI Y0= [0.01 0.49982777 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 19210
CI Y0= [0.1095 0.47970045 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 18640
CI Y0= [0.209 0.42939947 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 18255
CI Y0= [0.3085 0.35768553 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 16019
CI Y0= [0.408 0.27772293 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 14544
CI Y0= [0.5075 0.20567298 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 10538
CI Y0= [0.607 0.15896717 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 5742
CI Y0= [0.7065 0.1544447 0. 0. 0. 0. 0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 4511
CI Y0= [0.806 0.20655316 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 11086
CI Y0= [0.9055 0.32580487 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 16031
CI Y0= [1.005 0.51765367 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 18990
CI Y0= [1.1045 0.78191212 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 20432
CI Y0= [1.204 1.11277161 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 25320
CI Y0= [1.3035 1.49942293 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 32708
CI Y0= [1.403 1.92721127 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 26079
CI Y0= [1.5025 2.37920202 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 31995
CI Y0= [1.602 2.83799026 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 35436
CI Y0= [1.7015 3.2875604 0. 0. 0. 0. 0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 37012
CI Y0= [1.801 3.71499745 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 39771
CI Y0= [1.9005 4.11186698 0. 0. 0. 0.
0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 44507
CI Y0= [2. 4.4751163 0. 0. 0. 0. 0. ]
dY0= [ 0. 0. 0. 0. -10. 0. 0.]
res= [ 0. 0. 0. 0. 0. 0. -10.]
Successful function return.
Nbre d'appel a rhs: 48738
Tracer la position finale en fonction de x0
## BEGIN SOLUTION
plt.plot(X0,XF,'x')
plt.title('position finale')
## END SOLUTION
Text(0.5, 1.0, 'position finale')
2.7. Conclusion#
écrire vos commentaires et conclusion