Taller Python Científico - PyConES 2014

Parte II: Introducción a SymPy (en 15 min)

Juan Luis Cano, Kiko Correoso y Álex Sáez (@Pybonacci)

Fran Navarro, Zuria Bauer y Daniel Domene (@CAChemEorg)

Licencia Creative Commons
<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Taller de Python Científico</span> por <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Pybonacci y CAChemE</span> que se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.

Introducción al cálculo simbólico

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:

  • SymPy - Biblioteca o librería de Python que se puede utilizar en cualquier programa de Python, o en una sesión de IPython y que proporciona poderosas funciones de CAS.
  • Sage - Sage es un entorno CAS completo y muy potente que tiene como objetivo proporcionar un sistema de código abierto compitiendo con Mathematica y Maple. Sage no es una librería de Python regular, sino más bien un entorno de CAS que utiliza Python como lenguaje de programación.

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:

  • Crear símbolos y expresiones.
  • 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

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.

Empezando con SymPy

Creación de símbolos

Lo primero, como siempre, es importar aquello que vayamos a necesitar. En SymPy hay un comando muy sencillo que nos facilita esta tarea:

In [1]:
from sympy import init_session
In [2]:
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:

  • Gracias a use_latex=True obtenemos la salida en $\LaTeX$.
  • Ha creado una serie de variables para que podamos ponernos a trabajar en el momento.
Nota: En Python, no se declaran las variables, sin embargo, no puedes usar una hasta que no le hayas asignado un valor. Si ahora intentamos crear una variable `a` que sea `a = 2 * b`, veamos qué ocurre:
In [3]:
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:

¿Variables simbólicas?

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:

In [4]:
#creamos la variable/incóginta/símbolo 'x'
x = Symbol('x')
x**2
Out[4]:
$$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:

In [5]:
(x + pi)**2
Out[5]:
$$\left(x + \pi\right)^{2}$$

Igual que con Python, se pueden definir simbolos de forma alternativa y más compacta de la siguiente forma:

In [6]:
x, y, z, t = symbols('x y z t')

Para ver qué tipo de variable es $x$:

In [7]:
type(x)
Out[7]:
sympy.core.symbol.Symbol

Recuerda que seguimos trabajando con Python, por lo que si vuelvo a asignar un valor nuevo a $x$:

In [8]:
x = 3
x**2
Out[8]:
$$9$$

Las conclusiones son:

  • Si quiero usar una variable como símbolo debo crearla previamente.
  • Las operaciones con símbolos devuelven símbolos.
  • Si una variable que almacenaba un símbolo recibe otra asignación, cambia de tipo.

Nota: Por defecto, SymPy entiende que los símbolos son números complejos. Esto puede producir resultados inesperados ante determinadas operaciones como, por ejemplo, lo logaritmos. Podemos indicar que la variable es real, entera... en el momento de la creación:

In [9]:
#creamos símbolos reales
x, y, z, t = symbols('x y z t', real=True)
In [10]:
#preguntamos si son reales
x.is_integer is None
Out[10]:
True
__Las variables de tipo `Symbol` actúan como contenedores en los que no sabemos qué hay (un real, un complejo, una lista...), imagínatelas como incógnitas__.

Expresiones y manipulaciones algebraicas

Comencemos por crear una expresión como: $\cos(x)^2+\sin(x)^2$

In [11]:
expresion = cos(x)**2 + sin(x)**2
expresion
Out[11]:
$$\sin^{2}{\left (x \right )} + \cos^{2}{\left (x \right )}$$

¿Qué está ocurriendo? Python detecta que a es una variable de tipo Symbol y al operar con ella nos devuelve otra variable de Sympy.

Nota: No debes olvidar en ningún momento que estás utilizando Python, por lo tanto __el operador `=` está realizando una asignación__. La variable `expr` pasa a contener `cos(x)**2 + sin(x)**2`, pero esto no es una ecuación, más adelante veremos cómo definirlas.

simplify()

Podemos pedirle que simplifique la expresión anterior:

In [12]:
simplify(expresion)
Out[12]:
$$1$$

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()

In [13]:
#ver pi con 25 decimales
pi.evalf(25)
Out[13]:
$$3.141592653589793238462643$$
In [14]:
#el mismo resultado se obtiene ocn la función N()
N(pi,25)
Out[14]:
$$3.141592653589793238462643$$

subs()

El método subs también puede utilizarse para sustituir los símbolos en expresiones:

In [15]:
expresion = (sin(x) + 3 * x)
expresion.subs(x, pi)
Out[15]:
$$3 \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:

In [16]:
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]$

In [17]:
plot(x**2, (x, -5, 5))
Out[17]:
<sympy.plotting.plot.Plot at 0x80947f0>

Podemos además representar funciones en 3D de forma análoga:

In [18]:
from sympy.plotting import plot3d

plot3d(x*y, (x, -5, 5), (y, -5, 5))
Out[18]:
<sympy.plotting.plot.Plot at 0x83e8c50>

Cálculo

Además de las manipulaciones algebraicas, el otro uso principal de CAS es hacer el cálculo, como derivadas, integrales o límites...

Derivadas

Puedes derivar una expresion usando el método .diff() y la función dif()

In [19]:
#creamos una expresión
expresion = cos(x)

#obtenemos la derivada primera con funcion
diff(expresion, x)
Out[19]:
$$- \sin{\left (x \right )}$$
In [20]:
#utilizando método
expresion.diff(x)
Out[20]:
$$- \sin{\left (x \right )}$$

¿varias variables?

In [21]:
expr_xy = y**3 *sin(x)**2 + x**2 * cos(y)
expr_xy
Out[21]:
$$x^{2} \cos{\left (y \right )} + y^{3} \sin^{2}{\left (x \right )}$$
In [22]:
diff(expr_xy,x,2,y,2)
Out[22]:
$$2 \left(- 6 y \sin^{2}{\left (x \right )} + 6 y \cos^{2}{\left (x \right )} - \cos{\left (y \right )}\right)$$

Si queremos que la deje indicada, usamos Derivative()

In [23]:
Derivative(expr_xy, x,2,y)
Out[23]:
$$\frac{\partial^{3}}{\partial x^{2}\partial y} \left(x^{2} \cos{\left (y \right )} + y^{3} \sin^{2}{\left (x \right )}\right)$$

Integración

La integración de expresiones se realiza de una manera similar:

In [24]:
expresion = x**2
integrate(expresion, x)
Out[24]:
$$\frac{x^{3}}{3}$$

Al proporcionar límites para la variable de integración podemos evaluar integrales definidas:

In [25]:
expresion = x**3 - x**2 + x - 1

integrate(expresion, (x, -1, 1))
Out[25]:
$$- \frac{8}{3}$$

así como integrales impropias:

In [26]:
integrate(exp(-x**2), (x, -oo, oo))
Out[26]:
$$\sqrt{\pi}$$
En SimPy, `oo` se utiliza como notación para indicar el infinito.

Series

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()

In [27]:
#creamos la expresión
expr = exp(x)
expr
Out[27]:
$$e^{x}$$
In [28]:
#la desarrollamos en serie
series(expr)
Out[28]:
$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$$

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.

In [29]:
series(expr, n=10)
Out[29]:
$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \frac{x^{6}}{720} + \frac{x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{x^{9}}{362880} + \mathcal{O}\left(x^{10}\right)$$

Si nos molesta el $\mathcal{O}(x^{10})$ lo podemos quitar con removeO():

In [30]:
series(expr, n=10).removeO()
Out[30]:
$$\frac{x^{9}}{362880} + \frac{x^{8}}{40320} + \frac{x^{7}}{5040} + \frac{x^{6}}{720} + \frac{x^{5}}{120} + \frac{x^{4}}{24} + \frac{x^{3}}{6} + \frac{x^{2}}{2} + x + 1$$

Ejercicio.

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.

In [41]:
# 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.

In [49]:
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).

In [50]:
interact(taylor_graf, n=(0,10));

Si estás viendo este IPython Notebook "offline", este es el resultado interactivo:


Matrices

In [51]:
#creamos una matriz llena de símbolos
a, b, c, d = symbols('a b c d')
A = Matrix([[a, b],
            [c, d]])
A
Out[51]:
$$\left[\begin{matrix}a & b\\c & d\end{matrix}\right]$$
In [52]:
#sacamos autovalores
A.eigenvals()
Out[52]:
$$\begin{Bmatrix}\frac{a}{2} + \frac{d}{2} - \frac{1}{2} \sqrt{a^{2} - 2 a d + 4 b c + d^{2}} : 1, & \frac{a}{2} + \frac{d}{2} + \frac{1}{2} \sqrt{a^{2} - 2 a d + 4 b c + d^{2}} : 1\end{Bmatrix}$$
In [53]:
#inversa
A.inv()
Out[53]:
$$\left[\begin{matrix}\frac{1}{a} + \frac{b c}{a^{2} \left(d - \frac{b c}{a}\right)} & - \frac{b}{a \left(d - \frac{b c}{a}\right)}\\- \frac{c}{a \left(d - \frac{b c}{a}\right)} & \frac{1}{d - \frac{b c}{a}}\end{matrix}\right]$$

Esto es sólo la punta del iceberg

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:

This website does not host notebooks, it only renders notebooks available on other websites.

Delivered by Fastly, Rendered by OVHcloud

nbviewer GitHub repository.

nbviewer version: 90c61cc

nbconvert version: 5.6.1

Rendered (Sun, 18 Apr 2021 12:21:44 UTC)