from sympy import init_printing
init_printing(use_latex=True)
from sympy import Symbol, symbols
from sympy.abc import a, b, c
Documentación: http://docs.sympy.org/dev/modules/physics/mechanics/
Álgebra vectorial: http://docs.sympy.org/dev/modules/physics/mechanics/vectors.html#using-vectors-and-reference-frames-in-sympy-s-mechanics
Con SymPy podemos definir sistemas de referencia vectoriales:
from sympy.physics.mechanics import ReferenceFrame
A = ReferenceFrame('A', latexs=['\mathbf{i}', '\mathbf{j}', '\mathbf{k}'])
A.x + A.y + A.z
Ya podemos hacer operaciones básicas con los tres vectores coordenados:
type(A.x)
sympy.physics.mechanics.essential.Vector
2 * A.x
# También usando símbolos
a * (A.x + 2 * A.y)
4 + A.z # Type Error
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-8-ac698b78348d> in <module>() ----> 1 4 + A.z # Type Error /home/juanlu/.local/lib/python3.3/site-packages/sympy/physics/mechanics/essential.py in __add__(self, other) 1098 def __add__(self, other): 1099 """The add operator for Vector. """ -> 1100 other = _check_vector(other) 1101 return Vector(self.args + other.args) 1102 /home/juanlu/.local/lib/python3.3/site-packages/sympy/physics/mechanics/essential.py in _check_vector(other) 1899 other = sympify(other) 1900 if other != S(0): -> 1901 raise TypeError('A Vector must be supplied') 1902 else: 1903 other = Vector([]) TypeError: A Vector must be supplied
Los productos escalar y vectorial se pueden efectuar de varias maneras:
from sympy.physics.mechanics import dot, cross
dot(A.x, A.y)
dot(A.x, A.x)
A.x & A.y
A.x & A.x
cross(A.x, A.x)
cross(A.x, A.y)
A.x ^ A.x
A.x ^ A.y
Y podemos también obtener la magnitud de vectores o normalizarlos:
v = -A.x + a * A.y
v
v.magnitude()
v.normalize()
También podemos efectuar cálculo con vectores, utilizando sus métodos:
v.diff(a, A) # Se especifica el sistema de referencia
Nota: El método diff
de SymPy no funciona sobre estos vectores porque se debe especificar el sistema de referencia.
Podemos operar con vectores de sistemas de referencia distintos:
A1 = ReferenceFrame('A_1', latexs=['\mathbf{i_1}', '\mathbf{j_1}', '\mathbf{k_1}'])
A.x + A1.x
La información de los vectores se almacena en args
: ¡me encanta la implementación!
(A.x + A1.x).args
[(Matrix([ [1], [0], [0]]), A), (Matrix([ [1], [0], [0]]), A_1)]
Es el momento de orientar nuestro recién definido sistema de referencia con respecto a los ejes que ya estábamos manejando. Para eso usaremos el método orient
, y recuperaremos la matriz de cosenos directores usando el método dcm
:
phi = Symbol('phi')
A1.orient(A, 'Axis', [phi, A.z]) # Rotación phi alrededor del eje A.z
A1.dcm(A) # "Direct Cosine Matrix"
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$.Defino unos ejes fijos y unos ejes que rotan con velocidad $\Omega$. Para ello uso la función dynamicsymbols
, que crea una función indeterminada de t
. Es lo más adecuado a la hora de representar cualquier magnitud variable con el tiempo: coordenadas, fuerzas...
from sympy.physics.mechanics import dynamicsymbols
Omega = dynamicsymbols('Omega')
XYZ = ReferenceFrame("u", latexs=['\mathbf{i}', '\mathbf{j}', '\mathbf{k}'])
X1Y1Z = XYZ.orientnew("v", 'Axis', [Omega, XYZ.z], latexs=['\mathbf{i}_1', '\mathbf{j}_1', '\mathbf{k}_1'])
X1Y1Z.dcm(XYZ)
Vector de posición en ejes móviles:
R = Symbol('R')
r = R * X1Y1Z.x
r
Derivada del vector en ejes móviles (nula)
r.dt(X1Y1Z)
# Equivalente, mejor no usar
t = Symbol('t')
r.diff(t, X1Y1Z)
Derivada del vector en ejes fijos
v = r.dt(XYZ)
v.simplify()
v
Expresamos el vector $v$ en los ejes $XYZ$:
v.express(XYZ)
v.express(XYZ).args
[(Matrix([ [-R*sin(Omega(t))*Derivative(Omega(t), t)], [ R*cos(Omega(t))*Derivative(Omega(t), t)], [ 0]]), u)]
¡Perfecto!
Nota: se aplicaron algunos parches
Cinemática de la pala de un rotor:
Vamos a considerar ese orden para los movimientos aunque no tiene porqué ser el caso siempre.
¿Cuál será la velocidad de un punto $P$ de esa línea de cuartos del perfil?
Algunas variables:
r = Symbol('r')
e_b, e_a = symbols('e_b, e_a')
from sympy.physics.mechanics import ReferenceFrame, Point
# Definimos nuestra propia clase para que los versores sean IJK
class IJKReferenceFrame(ReferenceFrame):
def __init__(self, name):
super().__init__(name, latexs=['\mathbf{{{}_{{{}}}}}'.format(idx, name) for idx in ("i", "j", "k")])
Sistema de referencia árbol:
A = IJKReferenceFrame("A")
O = Point("O")
O.set_vel(A, 0)
Vamos introduciendo sistemas de referencia y sus orígenes:
A1 = IJKReferenceFrame("A1")
psi = dynamicsymbols('psi')
A1.orient(A, 'Axis', [psi, A.z])
A1.dcm(A) # T_{A1A}
A1.ang_vel_in(A)
A2 = IJKReferenceFrame("A2")
beta = dynamicsymbols('beta')
A2.orient(A1, 'Axis', [beta, -A1.y])
A2.dcm(A1) # T_{A2A1}
E_b = O.locatenew('E_b', e_b * A1.x)
E_b.pos_from(O)
A2.ang_vel_in(A1)
A2.dcm(A) # T_{A2A}
A2.ang_vel_in(A)
A3 = IJKReferenceFrame("A3")
zeta = dynamicsymbols('zeta')
A3.orient(A2, 'Axis', [zeta, A2.z])
A3.dcm(A2) # T_{A3A2}
E_a = E_b.locatenew('E_a', (e_a - e_b) * A2.x)
E_a.pos_from(E_b)
A3.ang_vel_in(A2)
Definimos ahora nuestro punto $P$:
P = E_a.locatenew('P', (r - e_a) * A2.x)
P.pos_from(O)
Vayamos viendo incrementalmente cómo conseguir componer los movimientos en la manera que nos interese:
# En SymPy, v2pt_theory hace referencia al campo de velocidades
# de un sólido rígido. Los argumentos son:
# * `O`, punto del que conocemos la velocidad respecto a A
# * `A`, sistema de referencia donde queremos expresar la velocidad
# * `A1`, sistema de referencia donde están fijos ambos puntos
E_b.v2pt_theory(O, A, A1)