Curso de Python para Ingenieros Mecánicos

Por: Eduardo Vieira

Sympy - Sistema de algebra simbólica

SymPy es una biblioteca de Python para matemática simbólica. Apunta a convertirse en un sistema de algebra computacional (CAS) con todas sus prestaciones manteniendo el código tan simple como sea posible para manterlo comprensible y fácilmente extensible. SymPy está escrito totalmente en Python y no requiere bibliotecas adicionales. Este proyecto comenzó en 2005, fue lanzado al público en 2007 y a él han contribuido durante estos años cientos de personas.

Otros CAS conocidos son Mathematica y Maple, sin embargo ambos son software privativo y de pago. Aquí puedes encontrar una comparativa de SymPy con Maple.

Puede hacer cosas como:

Manipular expresiones (simplificación, expansión...)
Calcular derivadas e integrales.
Límites y desarrollos en serie.
Resolución de ecuaciones.
Resolción de EDOs.
Matrices

Sin embargo, SymPy no acaba aquí ni mucho menos...



In [2]:

    
import sympy as sym
%matplotlib inline

Para que los resultados sean formateados en $\LaTeX$ podemos usar:



In [3]:

    
sym.init_printing(use_latex=True)

Variables simbólicas

En SymPy podemos crear símbolos para las variables con las que deseamos trabajar. Podemos crear un nuevo símbolo usando la clase Symbol:



In [4]:

    
x = sym.Symbol('x')



In [5]:

    
(sym.pi + x)**2









    Out[5]:





$$\left(x + \pi\right)^{2}$$



In [6]:

    
# forma alternativa de definir (varios) símbolos
a, b, c = sym.symbols("a, b, c")



In [7]:

    
type(a)









    Out[7]:





sympy.core.symbol.Symbol

Podemos agregar algunas propiedades a los símbolos cuando son creados:



In [8]:

    
x = sym.Symbol('x', real=True)



In [9]:

    
x.is_imaginary









    Out[9]:





False



In [10]:

    
x = sym.Symbol('x', positive=True)



In [11]:

    
x > 0









    Out[11]:





$$\mathrm{True}$$

Números complejos

La unidad imaginaria es denotada por I en Sympy.



In [12]:

    
1+1*sym.I









    Out[12]:





$$1 + i$$



In [13]:

    
sym.I**2









    Out[13]:





$$-1$$



In [14]:

    
(x * sym.I + 1)**2









    Out[14]:





$$\left(i x + 1\right)^{2}$$

Números racionales

Existen tres tipos distintos de números en SymPy: Real, Rational, Integer:



In [15]:

    
r1 = sym.Rational(4,5)
r2 = sym.Rational(5,4)



In [16]:

    
r1









    Out[16]:





$$\frac{4}{5}$$



In [17]:

    
r1+r2









    Out[17]:





$$\frac{41}{20}$$



In [18]:

    
r1/r2









    Out[18]:





$$\frac{16}{25}$$

Evaluación numérica

SymPy usa una librería para trabajar con números con precisión arbitraria, y tiene expresiones SymPy predefinidas para varias constantes matemáticas, tales como: pi, e y oo para el infinito.

Para evaluar numéricamente una expresión podemos usar la función evalf (o bien N). Ésta usa un argumento n que especifíca el número de cifras significativas.



In [19]:

    
sym.pi.evalf(n=50)









    Out[19]:





$$3.1415926535897932384626433832795028841971693993751$$



In [20]:

    
y = (x + sym.pi)**2



In [21]:

    
sym.N(y, 5) # equivalente a evalf









    Out[21]:





$$\left(x + 3.1416\right)^{2}$$

Cuando evaluamos numéricamente expresiones a menudo deseamos substituir un símbolo por un valor numérico. En SymPy hacemos esto usando la función subs:



In [22]:

    
y.subs(x, 1.5)









    Out[22]:





$$\left(1.5 + \pi\right)^{2}$$



In [23]:

    
sym.N(y.subs(x, 1.5))









    Out[23]:





$$21.5443823618587$$

La función subs también puede ser usada para substituir símbolos y expresiones:



In [24]:

    
y.subs(x, a+sym.pi)









    Out[24]:





$$\left(a + 2 \pi\right)^{2}$$

También podemos combinar la evaluación numérica de expresiones con arreglos NumPy:



In [25]:

    
import numpy as np



In [26]:

    
x_vec = np.arange(0, 10, 0.1)



In [27]:

    
y_vec = np.array([sym.N(((x + sym.pi)**2).subs(x, xx)) for xx in x_vec])



In [28]:

    
import matplotlib.pyplot as plt
plt.plot(x_vec, y_vec);

Sin embargo, este tipo de evaluación numérica puede ser muy lenta, y existe una forma mucho más eficiente de realizar la misma tarea: Usar la función lambdify para "mapear" una expresión de Sympy a una función que es mucho más eficiente para la evaluación numérica:



In [29]:

    
f = sym.lambdify([x], (x + sym.pi)**2, 'numpy')  # el primer argumento es una lista de variables de las que la función f dependerá: en este caso sólo x -> f(x)



In [30]:

    
type(f)









    Out[30]:





function



In [31]:

    
y_vec = f(x_vec)  # ahora podemos pasar directamente un arreglo Numpy. Así f(x) es evaluado más eficientemente

La mayor eficiencia de usar funciones "lambdificadas" en lugar de usar evalación numérica directa puede ser significativa, a menudo de varios órdenes de magnitud. Aún en este sencillo ejemplo obtenemos un aumento de velocidad importante:



In [32]:

    
%%timeit

y_vec = np.array([sym.N(((x + sym.pi)**2).subs(x, xx)) for xx in x_vec])









    



32.5 ms ± 1.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)



In [33]:

    
%%timeit

y_vec = f(x_vec)









    



7.26 µs ± 30.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Manipulaciones algebráicas

Uno de los usos principales de un sistema de cálculo simbólico es que realiza manipulaciones algebráicas de expresiones. Por ejemplo, si queremos expandir un producto, factorizar una expresión, o simplificar un resultado. En esta sección presentamos las funciones para realizar estas operaciones básicas en SymPy.

Expand and factor

Primeros pasos en la manipulación algebráica



In [34]:

    
(x+1)*(x+2)*(x+3)









    Out[34]:





$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$



In [35]:

    
sym.expand((x+1)*(x+2)*(x+3))









    Out[35]:





$$x^{3} + 6 x^{2} + 11 x + 6$$

La función expand acepta varias argumentos clave con los que se puede indicar qué tipo de expansión deseamos realizar. Por ejemplo, para expandir expresiones trigonométricas, usamos el argumento clave trig=True:



In [36]:

    
sym.sin(a+b)









    Out[36]:





$$\sin{\left (a + b \right )}$$



In [37]:

    
sym.expand(sym.sin(a+b), trig=True)









    Out[37]:





$$\sin{\left (a \right )} \cos{\left (b \right )} + \sin{\left (b \right )} \cos{\left (a \right )}$$

Ver help(expand) para una descripción detallada de los distintos tipos de expansiones que la función expand puede realizar.

También podemos factorizar expresiones, usando la función factor de SymPy:



In [38]:

    
sym.factor(x**3 + 6 * x**2 + 11*x + 6)









    Out[38]:





$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$

Simplify

La función simplify intenta simplificar una expresión usando distinta técnicas. Existen también alternativas más específicas a la función simplify: trigsimp, powsimp, logcombine, etc.

El uso básico de estas funciones en el siguiente:



In [39]:

    
# simplify expande un producto
sym.simplify((x+1)*(x+2)*(x+3))









    Out[39]:





$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$



In [40]:

    
# simplify usa identidades trigonometricas
sym.simplify(sym.sin(a)**2 + sym.cos(a)**2)









    Out[40]:





$$1$$



In [41]:

    
sym.simplify(sym.cos(x)/sym.sin(x))









    Out[41]:





$$\frac{1}{\tan{\left (x \right )}}$$

apart and together

Podemos también manipular expresiones simbólicas que involucran fracciones usando las funciones apart y together. La primera de estas funciones separa una fracción en sus correspondientes fracciones parciales; la segunda hace todo lo contrario.



In [42]:

    
f1 = 1/((a+1)*(a+2))



In [43]:

    
f1









    Out[43]:





$$\frac{1}{\left(a + 1\right) \left(a + 2\right)}$$



In [44]:

    
sym.apart(f1)









    Out[44]:





$$- \frac{1}{a + 2} + \frac{1}{a + 1}$$



In [45]:

    
f2 = 1/(a+2) + 1/(a+3)



In [46]:

    
f2









    Out[46]:





$$\frac{1}{a + 3} + \frac{1}{a + 2}$$



In [47]:

    
sym.together(f2)









    Out[47]:





$$\frac{2 a + 5}{\left(a + 2\right) \left(a + 3\right)}$$

Podemos usar también Simplify:



In [48]:

    
sym.simplify(f2)









    Out[48]:





$$\frac{2 a + 5}{\left(a + 2\right) \left(a + 3\right)}$$

Cálculo

Además de realizar manipulaciones algebráicas, SimPy puede realizar operaciones de cálculo, tales como derivar y derivar expresiones.

Derivación

Derviar es usualmente algo simple. Usamos la función diff. El primer argumento es una expresión que será derivada, y el segundo argumento es el símbolo respecto al cual se realizará la derivada:



In [49]:

    
y









    Out[49]:





$$\left(x + \pi\right)^{2}$$



In [50]:

    
sym.diff(y**2, x)









    Out[50]:





$$4 \left(x + \pi\right)^{3}$$

Para calcular derivadas de orden superior podemos usar:



In [51]:

    
sym.diff(y**2, x, x)









    Out[51]:





$$12 \left(x + \pi\right)^{2}$$

o bien



In [52]:

    
sym.diff(y**2, x, 2) # hace lo mismo









    Out[52]:





$$12 \left(x + \pi\right)^{2}$$

Calculando la derivada de una función de varias variables:



In [53]:

    
x, y, z = sym.symbols("x,y,z")



In [54]:

    
f = sym.sin(x*y) + sym.cos(y*z)

$\frac{d^3f}{dxdy^2}$



In [55]:

    
sym.diff(f, x, 1, y, 2)









    Out[55]:





$$- x \left(x y \cos{\left (x y \right )} + 2 \sin{\left (x y \right )}\right)$$

Integracion

La integración se realiza de la misma forma:



In [56]:

    
f









    Out[56]:





$$\sin{\left (x y \right )} + \cos{\left (y z \right )}$$



In [57]:

    
sym.integrate(f, x)









    Out[57]:





$$x \cos{\left (y z \right )} + \begin{cases} 0 & \text{for}\: y = 0 \\- \frac{1}{y} \cos{\left (x y \right )} & \text{otherwise} \end{cases}$$

Si le damos los límites de integración calculamos la integral definida



In [58]:

    
sym.integrate(f, (x, -1, 1))









    Out[58]:





$$2 \cos{\left (y z \right )}$$

Tambien podemos calcular integrales impropias



In [59]:

    
from sympy import oo
# oo es el infinito



In [60]:

    
sym.integrate(sym.exp(-x**2), (x, -oo, oo))









    Out[60]:





$$\sqrt{\pi}$$

Sumatoria y productos

Podemos evaluar sumatorias usando: 'Sum'



In [61]:

    
n = sym.Symbol("n")



In [62]:

    
sym.Sum(1/n**2, (n, 1, 10))









    Out[62]:





$$\sum_{n=1}^{10} \frac{1}{n^{2}}$$



In [63]:

    
sym.Sum(1/n**2, (n,1, 10)).evalf()









    Out[63]:





$$1.54976773116654$$



In [64]:

    
sym.Sum(1/n**2, (n, 1, oo)).evalf()









    Out[64]:





$$1.64493406684823$$

Lo mismo con productos



In [65]:

    
sym.Product(n, (n, 1, 10)) # 10!









    Out[65]:





$$\prod_{n=1}^{10} n$$

Límites

Los límites se evalúan con la función limit



In [66]:

    
sym.limit(sym.sin(x)/x, x, 0)









    Out[66]:





$$1$$

Se puede usar el límite para verificar una derivada



In [67]:

    
f









    Out[67]:





$$\sin{\left (x y \right )} + \cos{\left (y z \right )}$$



In [68]:

    
sym.diff(f, x)









    Out[68]:





$$y \cos{\left (x y \right )}$$

La derivada por definicion es $$\frac{d}{{dx}}f\left( x \right) = \mathop {\lim }\limits_{h \to 0} \frac{{f\left( {x + h } \right) - f\left( x \right)}}{h }$$



In [69]:

    
h = sym.Symbol("h")



In [70]:

    
sym.limit((f.subs(x, x+h) - f)/h, h, 0)









    Out[70]:





$$y \cos{\left (x y \right )}$$

OK!

Podemos calcular los limites laterales usando el argumento dir:



In [71]:

    
sym.limit(1/x, x, 0, dir="+")









    Out[71]:





$$\infty$$



In [72]:

    
sym.limit(1/x, x, 0, dir="-")









    Out[72]:





$$-\infty$$

Series

Las series tambien son una de las funciones mas útiles de un CAS. En SymPy se puede calcular la serie de una funcion usando series



In [73]:

    
sym.series(sym.exp(x), x)









    Out[73]:





$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$$

Por defecto, se calcula la expansión alrededor de $x=0$ (serie de Taylor), pero podemos expandir la serie al rededor de otro punto



In [74]:

    
sym.series(sym.exp(x), x, 1)









    Out[74]:





$$e + e \left(x - 1\right) + \frac{e}{2} \left(x - 1\right)^{2} + \frac{e}{6} \left(x - 1\right)^{3} + \frac{e}{24} \left(x - 1\right)^{4} + \frac{e}{120} \left(x - 1\right)^{5} + \mathcal{O}\left(\left(x - 1\right)^{6}; x\rightarrow 1\right)$$

Tambien podemos especificar el orden



In [75]:

    
sym.series(sym.exp(x), x, 1, 10)









    Out[75]:





$$e + e \left(x - 1\right) + \frac{e}{2} \left(x - 1\right)^{2} + \frac{e}{6} \left(x - 1\right)^{3} + \frac{e}{24} \left(x - 1\right)^{4} + \frac{e}{120} \left(x - 1\right)^{5} + \frac{e}{720} \left(x - 1\right)^{6} + \frac{e}{5040} \left(x - 1\right)^{7} + \frac{e}{40320} \left(x - 1\right)^{8} + \frac{e}{362880} \left(x - 1\right)^{9} + \mathcal{O}\left(\left(x - 1\right)^{10}; x\rightarrow 1\right)$$

La expansion en serie incluye la aproximación de orden. Lo cual es útil para trabajar con series de diferente orden



In [76]:

    
s1 = sym.cos(x).series(x, 0, 5)
s1









    Out[76]:





$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} + \mathcal{O}\left(x^{5}\right)$$



In [77]:

    
s2 = sym.sin(x).series(x, 0, 2)
s2









    Out[77]:





$$x + \mathcal{O}\left(x^{2}\right)$$



In [78]:

    
sym.expand(s1 * s2)









    Out[78]:





$$x + \mathcal{O}\left(x^{2}\right)$$

Si queremos eliminar la aproximación de orden usamos el método removeO



In [79]:

    
sym.expand(s1.removeO() * s2.removeO())









    Out[79]:





$$\frac{x^{5}}{24} - \frac{x^{3}}{2} + x$$

Pero hay que tener en cuenta que esta no es la expansión correcta de $\cos(x)\sin(x)$ al $5$to orden:



In [80]:

    
(sym.cos(x)*sym.sin(x)).series(x, 0, 6)









    Out[80]:





$$x - \frac{2 x^{3}}{3} + \frac{2 x^{5}}{15} + \mathcal{O}\left(x^{6}\right)$$

Álgebra Lineal

Matrices

Las matrices se definen usando la clase Matrix:



In [81]:

    
m11, m12, m21, m22 = sym.symbols("m11, m12, m21, m22")
b1, b2 = sym.symbols("b1, b2")



In [82]:

    
A = sym.Matrix([[m11, m12],[m21, m22]])
A









    Out[82]:





$$\left[\begin{matrix}m_{11} & m_{12}\\m_{21} & m_{22}\end{matrix}\right]$$



In [83]:

    
b = sym.Matrix([[b1], [b2]])
b









    Out[83]:





$$\left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]$$

Con la clase Matrix podemos realizar las operaciones matriciales típicas



In [84]:

    
A**2









    Out[84]:





$$\left[\begin{matrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\end{matrix}\right]$$



In [85]:

    
A * b









    Out[85]:





$$\left[\begin{matrix}b_{1} m_{11} + b_{2} m_{12}\\b_{1} m_{21} + b_{2} m_{22}\end{matrix}\right]$$

Y tambien determinantes e inversas:



In [86]:

    
A.det()









    Out[86]:





$$m_{11} m_{22} - m_{12} m_{21}$$



In [87]:

    
A.inv()









    Out[87]:





$$\left[\begin{matrix}\frac{m_{22}}{m_{11} m_{22} - m_{12} m_{21}} & - \frac{m_{12}}{m_{11} m_{22} - m_{12} m_{21}}\\- \frac{m_{21}}{m_{11} m_{22} - m_{12} m_{21}} & \frac{m_{11}}{m_{11} m_{22} - m_{12} m_{21}}\end{matrix}\right]$$

resolviendo ecuaciones

Para resolver ecuaciones podemos utilizar la función solve



In [88]:

    
sym.solve(x**2 - 1, x)









    Out[88]:





$$\left [ -1, \quad 1\right ]$$



In [89]:

    
sym.solve(x**4 - x**2 - 1, x)









    Out[89]:





$$\left [ - i \sqrt{- \frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad i \sqrt{- \frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad - \sqrt{\frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad \sqrt{\frac{1}{2} + \frac{\sqrt{5}}{2}}\right ]$$

Sistemas de ecuaciones



In [90]:

    
sym.solve([x + y - 1, x - y - 1], [x,y])









    Out[90]:





$$\left \{ x : 1, \quad y : 0\right \}$$

La ecuacion puede estar en términos de otras variables simbólicas



In [91]:

    
sym.solve([x + y - a, x - y - c], [x,y])









    Out[91]:





$$\left \{ x : \frac{a}{2} + \frac{c}{2}, \quad y : \frac{a}{2} - \frac{c}{2}\right \}$$

Lecturas complementarias

http://sympy.org/en/index.html - The SymPy projects web page.
https://github.com/sympy/sympy - The source code of SymPy.
http://live.sympy.org - Online version of SymPy for testing and demonstrations.



In [92]:

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









    Out[92]:




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



El estilo se ha aplicado =)



In [ ]: