Taller de Python Científico por Pybonacci y CAChemE que se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.
Con este documento vas a aprender a calcular derivadas, integrales, límites... y resolver ecuaciones con Python de forma libre (¡y gratuita!) tal y como trabajan Mathematica, Maple o Mupad.
Hay dos sistemas de álgebra computacional (CAS, computer algebraic system) a destacar para Python:
Sage es, en algunos aspectos, más poderoso que SymPy, pero ambos ofrecen una funcionalidad CAS muy completa. La ventaja de SymPy es que es una librería más de Python y se puede usar directamente con código o integrar en con el bloc de notas IPython. En esta lección, por lo tanto, vamos a ver cómo utilizar SymPy con Notebooks de IPython.
Hoy veremos cómo:
SymPy se puede usar directamente online sin registro alguno en SymPy Live, Sympy Gamma o mediante Notebooks con registro gratuito en Cloud Sage. Existen otras alternativas gratutias (pero comerciales) online como Wakari, además el futuro de IPython pasa por Jupyter (notebooks colaborativos para Julia, Python y R). Si lo que quieres instalártelo en tu ordenador recomendamos Anaconda de Continuum.
Lo primero, como siempre, es importar aquello que vayamos a necesitar. En SymPy hay un comando muy sencillo que nos facilita esta tarea:
from sympy import init_session
init_session(use_latex=True)
IPython console for SymPy 0.7.5 (Python 3.4.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) Documentation can be found at http://www.sympy.org
WARNING: Hook shutdown_hook is deprecated. Use the atexit module instead.
Podemos ver el comando init_session
ha llevado a cabo algunas acciones por nostros:
use_latex=True
obtenemos la salida en $\LaTeX$.a = 2 * b
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-3-d57516069ecc> in <module>() ----> 1 a = 2 * b NameError: name 'b' is not defined
Como en b
no había sido creada, Python no sabe qué es b
. Esto mismo nos ocurre con los símbolos de SymPy. Antes de usar una variable, debo decir que es un símbolo y asignárselo:
A diferencia de lo que estamos acostumbrados, las variables en SymPy son incógnitas. Podemos crear un nuevo símbolo con la instancia de clase Symbol
:
#creamos la variable/incóginta/símbolo 'x'
x = Symbol('x')
x**2
¿Qué ha pasado aquí arriba? En primer lugar creamos el símbolo con la variable $x$ y después vimos su valor. Si te has fijado, a diferencia de código normal en Python (o MATLAB), $x^2$ no devuelve ningún valor (número, texto, verdadero o falso) simplemente devuelve esa incógnita al cuadrado en el que a priori no sabemos que contiene.
Así, escribiendo una expresión con esta variable no habrá evaluación numérica:
(x + pi)**2
Igual que con Python, se pueden definir simbolos de forma alternativa y más compacta de la siguiente forma:
x, y, z, t = symbols('x y z t')
Para ver qué tipo de variable es $x$:
type(x)
sympy.core.symbol.Symbol
Recuerda que seguimos trabajando con Python, por lo que si vuelvo a asignar un valor nuevo a $x$:
x = 3
x**2
Las conclusiones son:
#creamos símbolos reales
x, y, z, t = symbols('x y z t', real=True)
#preguntamos si son reales
x.is_integer is None
True
Comencemos por crear una expresión como: $\cos(x)^2+\sin(x)^2$
expresion = cos(x)**2 + sin(x)**2
expresion
¿Qué está ocurriendo? Python detecta que a es una variable de tipo Symbol
y al operar con ella nos devuelve otra variable de Sympy.
simplify()
¶Podemos pedirle que simplifique la expresión anterior:
simplify(expresion)
En este caso parece estar claro, lo que quiere decir más simple, pero como en cualquier CAS el comando simplify
puede no devolvernos la expresión que nosotros queremos. Cuando esto ocurra necesitaremos usar otras instrucciones.
evalf()
¶Podemos obtener el valor numérico haciendo uso del método .evalf()
#ver pi con 25 decimales
pi.evalf(25)
#el mismo resultado se obtiene ocn la función N()
N(pi,25)
El método subs
también puede utilizarse para sustituir los símbolos en expresiones:
expresion = (sin(x) + 3 * x)
expresion.subs(x, pi)
Si queremos representar estas expresiones podemos hacerlo con la función plot
. Para ello, cargamos las librerías necesarias mediante el siguiente código:
from sympy.plotting import plot
# Representa las figuras en línea con el documento (no en una ventana emergente)
%matplotlib inline
Representamos $f(x) = x^2$ $\forall \: x \in [-5,5]$
plot(x**2, (x, -5, 5))
<sympy.plotting.plot.Plot at 0x80947f0>
Podemos además representar funciones en 3D de forma análoga:
from sympy.plotting import plot3d
plot3d(x*y, (x, -5, 5), (y, -5, 5))
<sympy.plotting.plot.Plot at 0x83e8c50>
Además de las manipulaciones algebraicas, el otro uso principal de CAS es hacer el cálculo, como derivadas, integrales o límites...
Puedes derivar una expresion usando el método .diff()
y la función dif()
#creamos una expresión
expresion = cos(x)
#obtenemos la derivada primera con funcion
diff(expresion, x)
#utilizando método
expresion.diff(x)
¿varias variables?
expr_xy = y**3 *sin(x)**2 + x**2 * cos(y)
expr_xy
diff(expr_xy,x,2,y,2)
Si queremos que la deje indicada, usamos Derivative()
Derivative(expr_xy, x,2,y)
La integración de expresiones se realiza de una manera similar:
expresion = x**2
integrate(expresion, x)
Al proporcionar límites para la variable de integración podemos evaluar integrales definidas:
expresion = x**3 - x**2 + x - 1
integrate(expresion, (x, -1, 1))
así como integrales impropias:
integrate(exp(-x**2), (x, -oo, oo))
La expansión de la serie es también una de las características más útiles de un CAS. En SymPy podemos realizar un desarrollo en serie de una expresión mediante el método .series()
o la función series()
#creamos la expresión
expr = exp(x)
expr
#la desarrollamos en serie
series(expr)
Se puede especificar el número de términos pasándole un argumento n=...
. El número que le pasemos será el primer término que desprecie.
series(expr, n=10)
Si nos molesta el $\mathcal{O}(x^{10})$ lo podemos quitar con removeO()
:
series(expr, n=10).removeO()
Como estás viendo, Python dispone de muchas librerías científicas. Vamos a comrpobar lo bien que se llevan los widgets de IPython Notebook 2.0 y SymPy.
# Cargamos widgets interactivos de IPython Notebook
from IPython.html.widgets import interact
from IPython.display import display
Obtén el desarrollo del $e^x$ y representa el resultado con la función real.
def taylor_graf(n):
expresion = #completar#
#calculamos la expansion elimnando el termino del error
expresion_aprox = series(expresion, x, n=n+1).removeO()
p1 = plot(expresion, (x, -3, 3),
show=False, line_color='b')
p2 = plot(expresion_aprox, (x, -3, 3),
show=False, line_color='r')
# Haz la segunda función parte de la primera
p1.extend(p2)
# Nota: Este código se debe a que representar dos funciones
# juntas con diferentes colores aún no está implementado
# http://stackoverflow.com/questions/21429866/change-color-implicit-plot-sympy
p1.show()
Función exponencial $e^x$ (en azul) y la suma de los primeros n+1 términos de su serie de Taylor en torno a cero (en rojo).
interact(taylor_graf, n=(0,10));
Si estás viendo este IPython Notebook "offline", este es el resultado interactivo:
#creamos una matriz llena de símbolos
a, b, c, d = symbols('a b c d')
A = Matrix([[a, b],
[c, d]])
A
#sacamos autovalores
A.eigenvals()
#inversa
A.inv()
SymPy es una herramienta con posibilidades infinitas, aquí hemos visto sus funcionalidades más básicas. Hay muchas, muchas cosas que se pueden hacer. Si buscas material en español: