

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

La manera usual de empezar a trabajar con SymPy es importar la función `init_session`:
```
from sympy import init_session
init_session(use_latex=True)```

 Esta función ya se encarga de importar todas las funciones básicas y preparar las salidas gráficas. Sin embargo, en este momento, esta función se encuentra en mantenimiento para su uso dentro de los notebooks por lo que activaremos la salida gráfica e importaremos las funciones de la manera usual. Puedes consultar el estado de la corrección en: https://github.com/sympy/sympy/pull/13300 y https://github.com/sympy/sympy/issues/13319 .


```python
from sympy import init_printing
init_printing()
```


```python
from sympy import symbols
```

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


```python
from sympy.physics.mechanics import ReferenceFrame
```


```python
A = ReferenceFrame("A")
A.x
```

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


```python
2 * A.x - 1 * 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$.


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

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


```python
# Definimos nuestra propia clase para que los versores sean IJK
# aeropython: preserve
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
```


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

### Álgebra vectorial

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


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


```python
r1
```


```python
v1
```


```python
from sympy.physics.mechanics import dot, cross
```


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


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

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


```python
(r1 ^ v1).magnitude()
```


```python
(r1 ^ v1).normalize()
```

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


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


```python
a = R * A0.i
a
```


```python
omega01 = Omega * A0.k
omega01
```


```python
da = omega01 ^ a
da
```


```python
da.magnitude()
```

<div class="alert alert-warning">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
</div>

### 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`.


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


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




$$\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$.

<div class="alert alert-success">¿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.</div>

#### 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`:


```python
A0.x.express(A1)
```


```python
A0.x.to_matrix(A1)
```




$$\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`:


```python
from sympy.physics.mechanics import dynamicsymbols
```


```python
alpha = dynamicsymbols('alpha')
alpha
```

Y pedir su derivada con el método `diff`:


```python
alpha.diff()
```

##### Ejercicio



(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$.**


```python
A = IJKReferenceFrame("A")
```


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




$$\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]$$




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




$$\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]$$




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




$$\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]$$




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




$$\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]$$




```python
B.dcm(A2)
```




$$\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]$$




```python
B.dcm(A1)
```




$$\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`:


```python
B.ang_vel_in(A2)
```


```python
B.ang_vel_in(A)
```


```python
B.ang_vel_in(A).express(A)
```

En ocasiones, la representación gráfica puede fallar, pero se puede volver a desactivar y activar llamando a la función`init_printing(pretty_print=True)` con diferentes valores (True/False) para `pretty_print`

### Derivada en ejes móviles

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


```python
v1 = A1.x
v1
```


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


```python
dv1.to_matrix(A1)
```




$$\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]$$




```python
(dv1 & A1.j).simplify()
```

### 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`.


```python
from sympy.physics.mechanics import Point
```


```python
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`:


```python
O.set_vel(A, 0)
```

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


```python
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`:


```python
E_b.pos_from(O)
```

<div class="alert alert-info">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</div>

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:


```python
E_b.v2pt_theory(O, A, A1)
```

##### Ejercicio



(Apuntes de Óscar López Rebollal)

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


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


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


```python
# 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])
```


```python
A2.dcm(A1)
```




$$\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]$$




```python
# 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)
```


```python
# 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)
```


```python
# 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)
```

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

<a href="https://twitter.com/share" class="twitter-share-button" data-url="https://github.com/AeroPython/Curso_AeroPython" data-text="Aprendiendo Python con" data-via="pybonacci" data-size="large" data-hashtags="AeroPython">Tweet</a>


---

#### <h4 align="right">¡Síguenos en Twitter!

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

##### <a rel="license" href="http://creativecommons.org/licenses/by/4.0/deed.es"></a><br /><span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Curso AeroPython</span> por <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo</span> se distribuye bajo una <a rel="license" href="http://creativecommons.org/licenses/by/4.0/deed.es">Licencia Creative Commons Atribución 4.0 Internacional</a>.

#####    

---
_Las siguientes celdas contienen configuración del Notebook_

_Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como [seguro](http://ipython.org/ipython-doc/dev/notebook/security.html)_

    File > Trusted Notebook


```python
%%html
<a href="https://twitter.com/Pybonacci" class="twitter-follow-button" data-show-count="false">Follow @Pybonacci</a>

```


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




```python
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/aeropython.css'
HTML(open(css_file, "r").read())
```




/* This template is inspired in the one used by Lorena Barba
in the numerical-mooc repository: https://github.com/numerical-mooc/numerical-mooc
We thank her work and hope you also enjoy the look of the notobooks with this style */

<link href='http://fonts.googleapis.com/css?family=Source+Sans+Pro|Josefin+Sans:400,700,400italic|Ubuntu+Condensed' rel='stylesheet' type='text/css'>

El estilo se ha aplicado =)

<style>



#notebook_panel { /* main background */
    background: #f7f7f7;
}

div.cell { /* set cell width */
    width: 900px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    width: 950px;
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
    margin-top:0.7em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running { 
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    font-family: 'Source Sans Pro', sans-serif;
    background-color: rgb(256,256,256);
    font-size: 110%;
    border-radius: 0px; 
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Josefin Sans', serif;
    line-height: 145%;
    font-size: 125%;
    font-weight: 500;
    width:750px;
    margin-left:auto;
    margin-right:auto;
}


/* Formatting for header cells */
.text_cell_render h1, .text_cell_render h2, .text_cell_render h3,
.text_cell_render h4, .text_cell_render h5 {
    font-family: 'Ubuntu Condensed', sans-serif;
}
/*
.text_cell_render h1 {
    font-family: Flux, 'Ubuntu Condensed', serif;
    font-style:regular;
    font-weight: 400;    
    font-size: 30pt;
    text-align: center;
    line-height: 100%;
    color: #335082;
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}
*/
.text_cell_render h1 {
    font-weight: 600;
    font-size: 35pt;
    line-height: 100%;
    color: #000000;
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}

.text_cell_render h2 {
    margin-top:16px;
    font-size: 27pt;
    font-weight: 550;
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    font-style: regular;
    color: #2c6391;
}	

.text_cell_render h3 {
    font-size: 20pt;
    font-weight: 550
    text-align: left;
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    font-style: regular;
    color:  #387eb8;
}

.text_cell_render h4 {    /*Use this for captions*/
    font-size: 18pt;
    font-weight: 450
    text-align: left;
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    font-style: regular;
    color:  #5797cc;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-size: 18pt;
    font-weight: 550;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
    color:  #b21c0d;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'Ubuntu Condensed', sans-serif;
    font-weight: 300;
    font-size: 14pt;
    line-height: 100%;
    color: #252525;
    text-align: right;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
        font-family: 'Duru Sans', sans-serif;
        font-size: 100%;
}

</style>





ceba149fa33d3793c297b0574322588892c7d0cb
ipynb
Jupyter Notebook
notebooks_completos/041-SymPy-Mecanica.ipynb
caorodriguez/aprendealgo-
b8e9692214280cc0c1981bbd414b19fd0fa7dcb8
CC-BY-4.0
2015-03-05T18:57:16.000Z
2021-11-18T00:19:22.000Z
notebooks_completos/041-SymPy-Mecanica.ipynb
caorodriguez/aprendealgo-
b8e9692214280cc0c1981bbd414b19fd0fa7dcb8
CC-BY-4.0
2015-01-26T19:08:20.000Z
2020-05-20T17:09:58.000Z
notebooks_completos/041-SymPy-Mecanica.ipynb
caorodriguez/aprendealgo-
b8e9692214280cc0c1981bbd414b19fd0fa7dcb8
CC-BY-4.0
2015-02-18T15:12:45.000Z
2022-03-29T06:18:39.000Z
Qwen/Qwen-72B
1. YES
2. YES