Clase 5b: Mecánica con SymPy¶

Si SymPy te ha parecido hasta ahora un CAS decente e incluso interesante (nada como tener los resultados en $\LaTeX$ incrustados en el notebook y la sintaxis de Python para hacer cálculo simbólico) entonces espera a ver el paquete mechanics. Con él, podremos manipular velocidades y aceleraciones de sólidos expresadas en distintos sistemas de referencia con una facilidad impresionante.

Tienes disponible la documentación de mechanics en http://docs.sympy.org/0.7.5/modules/physics/mechanics/index.html.

Sistemas de referencia¶

El objeto primordial que vamos a manejar van a ser los sistemas de referencia. Podremos definir relaciones geométricas entre ellos y de esta forma las transformaciones de vectores entre un sistema y otro serán triviales.

In [1]:

from sympy import init_session
init_session(use_latex=True)

IPython console for SymPy 1.0 (Python 3.6.1-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

Todo lo que necesitamos está en sympy.physics.mechanics, incluyendo la clase ReferenceFrame. Nada más crear un sistema de referencia podemos acceder a sus versores unitarios: x, y y z.

http://docs.sympy.org/0.7.5/modules/physics/vector/vectors.html

In [2]:

from sympy.physics.mechanics import ReferenceFrame

In [3]:

A = ReferenceFrame("A")
A.x

Out[3]:

$$\mathbf{\hat{a}_x}$$

Y para definir vectores solo tenemos que multiplicar cada componente por su versor:

In [4]:

2 * A.x - 1 * A.y

Out[4]:

$$2\mathbf{\hat{a}_x} - \mathbf{\hat{a}_y}$$

De ahora en adelante, para trabajar como si nos enfrentáramos a un problema de la escuela, vamos a hacer dos cosas:

Definir un sistema inercial $1$ del que partir, para así poder referir todos los demás sistemas a él.
Que los versores de ese sistema sean $i, j, k$.

In [5]:

A = ReferenceFrame("1", latexs=['\mathbf{i}', '\mathbf{j}', '\mathbf{k}'])
A.x + A.y + A.z

Out[5]:

$$\mathbf{i} + \mathbf{j} + \mathbf{k}$$

Y para no tener que hacerlo siempre, un pequeño truco de magia:

In [6]:

# Definimos nuestra propia clase para que los versores sean IJK
class IJKReferenceFrame(ReferenceFrame):
    def __init__(self, name):
        super().__init__(name, latexs=['\mathbf{%s}_{%s}' % (idx, name) for idx in ("i", "j", "k")])
        self.i = self.x
        self.j = self.y
        self.k = self.z

In [7]:

A = IJKReferenceFrame("1")
A.i + A.j + A.k

Out[7]:

$$\mathbf{i}_{1} + \mathbf{j}_{1} + \mathbf{k}_{1}$$

Álgebra vectorial¶

Nuestros vectores funcionan también con símbolos, y podemos realizar las operaciones de producto escalar y producto vectorial con ellos.

In [8]:

R, V = symbols('R, V', positive=True)
r1 = R * (A.x + A.y + A.z)
v1 = V * (A.x - 2 * A.z)

In [9]:

r1

Out[9]:

$$R\mathbf{i}_{1} + R\mathbf{j}_{1} + R\mathbf{k}_{1}$$

In [10]:

v1

Out[10]:

$$V\mathbf{i}_{1} - 2 V\mathbf{k}_{1}$$

In [11]:

from sympy.physics.mechanics import dot, cross

In [12]:

r1.dot(v1)
dot(r1, v1)
r1 & v1

Out[12]:

$$- R V$$

In [14]:

r1.cross(v1)
cross(r1, v1)
r1 ^ v1

Out[14]:

$$- 2 R V\mathbf{i}_{1} + 3 R V\mathbf{j}_{1} - R V\mathbf{k}_{1}$$

Podemos hallar también la norma de los vectores con su método magnitude e incluso normalizarlos con normalize:

In [15]:

(r1 ^ v1).magnitude()

Out[15]:

$$\sqrt{14} R V$$

In [16]:

(r1 ^ v1).normalize()

Out[16]:

$$- \frac{\sqrt{14}}{7}\mathbf{i}_{1} + \frac{3 \sqrt{14}}{14}\mathbf{j}_{1} - \frac{\sqrt{14}}{14}\mathbf{k}_{1}$$

Ejercicio¶

Usando directamente la fórmula para la derivada en ejes móviles:

$$\left(\frac{\operatorname{d}\!\mathbf{a}}{\operatorname{d}\!t}\right)_1 = \left(\frac{\operatorname{d}\!\mathbf{a}}{\operatorname{d}\!t}\right)_0 + \mathbf{\omega}_{01}\! \times \mathbf{a}$$

Calcula la derivada del vector de posición $R \mathbf{i}_0$, siendo $A_0$ un sistema de referencia que gira respecto al inercial con velocidad angular $\mathbf{\omega}_{01}=\Omega \mathbf{k}_0$. ¿Cuál es el módulo de la derivada?

In [17]:

R, Omega = symbols('R, Omega', positive=True)
A0 = IJKReferenceFrame('0')

In [18]:

a = R * A0.i
a

Out[18]:

$$R\mathbf{i}_{0}$$

In [19]:

omega01 = Omega * A0.k
omega01

Out[19]:

$$\Omega\mathbf{k}_{0}$$

In [20]:

da = omega01 ^ a
da

Out[20]:

$$\Omega R\mathbf{j}_{0}$$

In [21]:

da.magnitude()

Out[21]:

$$\Omega R$$

Si no especificaste `positive=True` vas a ver algo como $\sqrt{\Omega^2 R^2}$. Debería haber una forma de simplificar esta expresión _a posteriori_, pero de momento no funciona del todo bien. Preparando este notebook nos hemos dado cuenta y ya les hemos avisado :) https://github.com/sympy/sympy/issues/8326

Movimiento relativo¶

¿A quién no le gusta multiplicar matrices de rotación? Para esa minoría que lo detesta, existe SymPy. Para ello debemos especificar la orientación de nuestros sistemas de referencia usando el método orient, y recuperaremos la matriz de cosenos directores usando el método dcm.

In [22]:

A1 = IJKReferenceFrame("1")
A0 = IJKReferenceFrame("0")

In [23]:

phi = symbols('phi')
A0.orient(A1, 'Axis', [phi, A1.z])  # Rotación phi alrededor del eje A1.z
A0.dcm(A1)  # "Direct Cosine Matrix"

Out[23]:

$$\left[\begin{matrix}\cos{\left (\phi \right )} & \sin{\left (\phi \right )} & 0\\- \sin{\left (\phi \right )} & \cos{\left (\phi \right )} & 0\\0 & 0 & 1\end{matrix}\right]$$

Usando el argumento Axis hemos especificado que rotamos el sistema un ángulo especificado alrededor de un eje. Otros métodos son:

Body: se especifican los tres ángulos de Euler.
Space: igual que Body, pero las rotaciones se aplican en orden inverso.
Quaternion: utilizando cuaternios, rotación alrededor de un vector unitario $\lambda$ una cantidad $\theta$.

¿Qué es lo bueno de usar uno de estos métodos? ¡Que **siempre** tenemos la transformación bien definida! Es imposible meter "a capón" una matriz de rotación que sea incorrecta o absurda.

Diferente sistema de referencia¶

Para expresar un vector en otro sistema de referencia, no hay más que usar los métodos express o to_matrix:

In [24]:

A0.x.express(A1)

Out[24]:

$$\operatorname{cos}\left(\phi\right)\mathbf{i}_{1} + \operatorname{sin}\left(\phi\right)\mathbf{j}_{1}$$

In [25]:

A0.x.to_matrix(A1)

Out[25]:

$$\left[\begin{matrix}\cos{\left (\phi \right )}\\\sin{\left (\phi \right )}\\0\end{matrix}\right]$$

Símbolos dinámicos¶

Si queremos especificar que un símbolo puede variar con el tiempo, hay que usar la función dynamicsymbols:

In [26]:

from sympy.physics.mechanics import dynamicsymbols

In [27]:

alpha = dynamicsymbols('alpha')
alpha

Out[27]:

$$\alpha{\left (t \right )}$$

Y pedir su derivada con el método diff:

In [28]:

alpha.diff()

Out[28]:

$$\frac{d}{d t} \alpha{\left (t \right )}$$

Ejercicio¶

Sistemas de ejes asociados a la pala

(Sacado de Cuerva et al. "Teoría de los Helicópteros")

Obtener la matriz de rotación de la pala $B$ respecto a los ejes $A1$.

In [29]:

A = IJKReferenceFrame("A")

In [30]:

A1 = IJKReferenceFrame("A1")
psi = dynamicsymbols('psi')
A1.orient(A, 'Axis', [psi, A.z])
A1.dcm(A)  # T_{A1A}

Out[30]:

$$\left[\begin{matrix}\cos{\left (\psi{\left (t \right )} \right )} & \sin{\left (\psi{\left (t \right )} \right )} & 0\\- \sin{\left (\psi{\left (t \right )} \right )} & \cos{\left (\psi{\left (t \right )} \right )} & 0\\0 & 0 & 1\end{matrix}\right]$$

In [31]:

A2 = IJKReferenceFrame("A2")
beta = dynamicsymbols('beta')
A2.orient(A1, 'Axis', [beta, -A1.y])
A2.dcm(A1)  # T_{A2A1}

Out[31]:

$$\left[\begin{matrix}\cos{\left (\beta{\left (t \right )} \right )} & 0 & \sin{\left (\beta{\left (t \right )} \right )}\\0 & 1 & 0\\- \sin{\left (\beta{\left (t \right )} \right )} & 0 & \cos{\left (\beta{\left (t \right )} \right )}\end{matrix}\right]$$

In [32]:

A3 = IJKReferenceFrame("A3")
zeta = dynamicsymbols('zeta')
A3.orient(A2, 'Axis', [zeta, A2.z])
A3.dcm(A1)  # T_{A3A1}

Out[32]:

$$\left[\begin{matrix}\cos{\left (\beta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & \sin{\left (\zeta\left(t\right) \right )} & \sin{\left (\beta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )}\\- \sin{\left (\zeta\left(t\right) \right )} \cos{\left (\beta{\left (t \right )} \right )} & \cos{\left (\zeta\left(t\right) \right )} & - \sin{\left (\beta{\left (t \right )} \right )} \sin{\left (\zeta\left(t\right) \right )}\\- \sin{\left (\beta{\left (t \right )} \right )} & 0 & \cos{\left (\beta{\left (t \right )} \right )}\end{matrix}\right]$$

In [33]:

B = IJKReferenceFrame("B")
theta = dynamicsymbols('theta')
B.orient(A3, 'Axis', [theta, A3.x])
B.dcm(A3)  # T_{BA3}

Out[33]:

$$\left[\begin{matrix}1 & 0 & 0\\0 & \cos{\left (\theta{\left (t \right )} \right )} & \sin{\left (\theta{\left (t \right )} \right )}\\0 & - \sin{\left (\theta{\left (t \right )} \right )} & \cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

In [34]:

B.dcm(A2)

Out[34]:

$$\left[\begin{matrix}\cos{\left (\zeta\left(t\right) \right )} & \sin{\left (\zeta\left(t\right) \right )} & 0\\- \sin{\left (\zeta\left(t\right) \right )} \cos{\left (\theta{\left (t \right )} \right )} & \cos{\left (\theta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & \sin{\left (\theta{\left (t \right )} \right )}\\\sin{\left (\theta{\left (t \right )} \right )} \sin{\left (\zeta\left(t\right) \right )} & - \sin{\left (\theta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & \cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

In [35]:

B.dcm(A1)

Out[35]:

$$\left[\begin{matrix}\cos{\left (\beta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & \sin{\left (\zeta\left(t\right) \right )} & \sin{\left (\beta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )}\\- \sin{\left (\beta{\left (t \right )} \right )} \sin{\left (\theta{\left (t \right )} \right )} - \sin{\left (\zeta\left(t\right) \right )} \cos{\left (\beta{\left (t \right )} \right )} \cos{\left (\theta{\left (t \right )} \right )} & \cos{\left (\theta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & - \sin{\left (\beta{\left (t \right )} \right )} \sin{\left (\zeta\left(t\right) \right )} \cos{\left (\theta{\left (t \right )} \right )} + \sin{\left (\theta{\left (t \right )} \right )} \cos{\left (\beta{\left (t \right )} \right )}\\- \sin{\left (\beta{\left (t \right )} \right )} \cos{\left (\theta{\left (t \right )} \right )} + \sin{\left (\theta{\left (t \right )} \right )} \sin{\left (\zeta\left(t\right) \right )} \cos{\left (\beta{\left (t \right )} \right )} & - \sin{\left (\theta{\left (t \right )} \right )} \cos{\left (\zeta\left(t\right) \right )} & \sin{\left (\beta{\left (t \right )} \right )} \sin{\left (\theta{\left (t \right )} \right )} \sin{\left (\zeta\left(t\right) \right )} + \cos{\left (\beta{\left (t \right )} \right )} \cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

Velocidad angular¶

También podemos hallar la velocidad angular de un sistema respecto a otro usando el método ang_vel_in:

In [36]:

B.ang_vel_in(A2)

Out[36]:

$$\dot{\theta}\mathbf{i}_{A3} + \dot{\zeta\left(t\right)}\mathbf{k}_{A2}$$

In [37]:

B.ang_vel_in(A)

Out[37]:

$$\dot{\theta}\mathbf{i}_{A3} + \dot{\zeta\left(t\right)}\mathbf{k}_{A2} - \dot{\beta}\mathbf{j}_{A1} + \dot{\psi}\mathbf{k}_{A}$$

In [38]:

B.ang_vel_in(A).express(A)

Out[38]:

$$(\left(- \operatorname{sin}\left(\psi\right) \operatorname{sin}\left(\zeta\left(t\right)\right) + \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\psi\right) \operatorname{cos}\left(\zeta\left(t\right)\right)\right) \dot{\theta} - \operatorname{sin}\left(\beta\right) \operatorname{cos}\left(\psi\right) \dot{\zeta\left(t\right)} + \operatorname{sin}\left(\psi\right) \dot{\beta})\mathbf{i}_{A} + (\left(\operatorname{sin}\left(\psi\right) \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\zeta\left(t\right)\right) + \operatorname{sin}\left(\zeta\left(t\right)\right) \operatorname{cos}\left(\psi\right)\right) \dot{\theta} - \operatorname{sin}\left(\beta\right) \operatorname{sin}\left(\psi\right) \dot{\zeta\left(t\right)} - \operatorname{cos}\left(\psi\right) \dot{\beta})\mathbf{j}_{A} + (\operatorname{sin}\left(\beta\right) \operatorname{cos}\left(\zeta\left(t\right)\right) \dot{\theta} + \operatorname{cos}\left(\beta\right) \dot{\zeta\left(t\right)} + \dot{\psi})\mathbf{k}_{A}$$

Derivada en ejes móviles¶

Hacer una derivada con la fórmula lo hace cualquiera, pero SymPy puede encargarse automáticamente.

In [39]:

v1 = A1.x
v1

Out[39]:

$$\mathbf{i}_{A1}$$

In [40]:

#v1.diff(dynamicsymbols._t, A2)
dv1 = v1.diff(symbols('t'), A)
dv1

Out[40]:

$$(\operatorname{sin}^{2}\left(\psi\right) \dot{\psi} + \operatorname{cos}^{2}\left(\psi\right) \dot{\psi})\mathbf{j}_{A1}$$

In [41]:

dv1.to_matrix(A1)

Out[41]:

$$\left[\begin{matrix}0\\\sin^{2}{\left (\psi{\left (t \right )} \right )} \frac{d}{d t} \psi{\left (t \right )} + \cos^{2}{\left (\psi{\left (t \right )} \right )} \frac{d}{d t} \psi{\left (t \right )}\\0\end{matrix}\right]$$

In [42]:

(dv1 & A1.j).simplify()

Out[42]:

$$\frac{d}{d t} \psi{\left (t \right )}$$

Puntos, velocidades y la rueda que no desliza¶

El último paso que nos queda para completar la cinemática es la posibilidad de definir puntos en sólidos y aplicar su campo de velocidades. SymPy también permite esto, y para ello no tenemos más que importar la clase Point.

In [43]:

from sympy.physics.mechanics import Point

In [44]:

O = Point("O")

Para trabajar como lo haríamos en la escuela, vamos a especificar que $O$ es el origen de $A$, y para eso vamos a imponer que su velocidad es cero con el método set_vel:

In [45]:

O.set_vel(A, 0)

Para definir nuevos puntos, podemos utilizar el método locate_new:

In [46]:

e_b = symbols('e_b')
E_b = O.locatenew('E_b', e_b * A1.x)

Y para obtener vectores de un punto a otro, el método pos_from:

In [47]:

E_b.pos_from(O)

Out[47]:

$$e_{b}\mathbf{i}_{A1}$$

La notación de este paquete está influenciada por el libro Kane, T. R. & Levinson, D. A. "Dynamics, Theory and Applications". Es ligeramente distinto a como estudiamos nosotros en la escuela, pero ¡están abiertos a que les hagamos cualquier tipo de sugerencia! https://github.com/sympy/sympy/issues/2584#issuecomment-31552654

Por último, el campo de velocidades de un sólido rígido se formula usando el método v2pt_theory.

$$v^P_A = v^O_A + \omega_{A_1 A} \times \mathbf{OP}$$

Este método pertenece al punto del cual queremos conocer la velocidad y recibe tres parámetros:

O, punto de velocidad conocida respecto a A
A, sistema de referencia donde queremos calcular la velocidad
A1, sistema de referencia donde están fijos ambos puntos (sistema de arrastre)

Por tanto, para hallar la velocidad del punto que acabamos de crear:

In [48]:

E_b.v2pt_theory(O, A, A1)

Out[48]:

$$e_{b} \dot{\psi}\mathbf{j}_{A1}$$

Ejercicio¶

Disco que rueda

(Apuntes del profesor Dr. Óscar López Rebollal)

¡Halla la velocidad y la aceleración de $P$!

In [49]:

# Creamos nuestros sistemas de referencia
A1 = IJKReferenceFrame('1')
A0 = IJKReferenceFrame('0')
A2 = IJKReferenceFrame('2')

In [50]:

# Creamos los símbolos dinámicos necesarios
xi, theta = dynamicsymbols('xi, theta')
xi, theta

Out[50]:

$$\left ( \xi{\left (t \right )}, \quad \theta{\left (t \right )}\right )$$

In [51]:

# Orientamos los sistemas de referencia
A0.orient(A1, 'Axis', [0, A1.k])  # A0 no gira respecto a A1
A2.orient(A0, 'Axis', [theta, A0.k])

In [52]:

A2.dcm(A1)

Out[52]:

$$\left[\begin{matrix}\cos{\left (\theta{\left (t \right )} \right )} & \sin{\left (\theta{\left (t \right )} \right )} & 0\\- \sin{\left (\theta{\left (t \right )} \right )} & \cos{\left (\theta{\left (t \right )} \right )} & 0\\0 & 0 & 1\end{matrix}\right]$$

In [53]:

# Creamos el punto C, centro del disco, y especificamos su velocidad
# respecto a A1
C = Point('C')
C.set_vel(A1, xi.diff() * A1.x)

In [54]:

# Localizamos el punto P, punto fijo del disco, respecto a C, en
# el sistema A2 (que gira solidariamente con el disco)
R = symbols('R')
P = C.locatenew('P', -R * A2.j)
P.pos_from(C)

Out[54]:

$$- R\mathbf{j}_{2}$$

In [55]:

# Hallamos la velocidad de P en A1, expresada en A0
# ¡Con esta llamada ya estamos diciendo que C y P son fijos en A2!
P.v2pt_theory(C, A1, A2).express(A0)

Out[55]:

$$(R \operatorname{cos}\left(\theta\right) \dot{\theta} + \dot{\xi})\mathbf{i}_{0} + R \operatorname{sin}\left(\theta\right) \dot{\theta}\mathbf{j}_{0}$$

Misión cumplida :)

Hemos hecho un repaso bastante profundo de las posibilidades del paquete mechanics de SymPy. Nos hemos dejado algunas cosas en el tintero pero no demasiadas: esta funcionalidad aún se está expandiendo y necesita pulir algunos detalles.

Referencias

Capítulo de aeromecánica del libro de Cuerva y otros http://nbviewer.ipython.org/gist/Juanlu001/7711865
Estabilidad longitudinal de un Boeing 747 http://nbviewer.ipython.org/github/AlexS12/Mecanica_Vuelo/blob/master/MVII_MatrizSistema.ipynb

¿Serás tú el siguiente que publique un notebook usando SymPy? ;)

Si te ha gustado esta clase:

¡Síguenos en Twitter!

Follow @Pybonacci <a href="https://twitter.com/AlexS12" class="twitter-follow-button" data-show-count="false" align="right";>Follow @AlexS12 <a href="https://twitter.com/newlawrence" class="twitter-follow-button" data-show-count="false" align="right";>Follow @newlawrence

Clase 5b: Mecánica con SymPy¶

Sistemas de referencia¶

Álgebra vectorial¶

Ejercicio¶

Movimiento relativo¶

Diferente sistema de referencia¶

Símbolos dinámicos¶

Ejercicio¶

Velocidad angular¶

Derivada en ejes móviles¶

Puntos, velocidades y la rueda que no desliza¶

Ejercicio¶

¡Síguenos en Twitter!

Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.¶