Cálculo simbólico con Sympy

Sympy permite hacer operaciones analíticas o con símbolos en lugar de con valores numéricos Al igual que en Python existen varios tipos de datos numéricos como enteros (int), decimales (float) o booleanos (bool:True, False, etc.), Sympy posee tres tipos de datos propios: Real, Rational e Integer, es decir, números reales, racionales y enteros. Estoquiere decir que Rational(1,2) representa 1/2, Rational(5,2) a 5/2, etc. en lugar de 0.5 o 2.5.

>>> import sympy as sp
>>> a = sp.Rational(1,2)

>>> a
1/2

>>> a*2
1

>>> sp.Rational(2)**50/sp.Rational(10)**50
1/88817841970012523233890533447265625

También existen algunas constantes especiales, como el número \(e\) o \(\pi\), si embargo éstos se tratan con símbolos y no tienen un valor númerico determinado. Eso quiere decir que no se puede obtener un valor numérico de una operación usando el pi de Sympy, como (1+pi), como lo haríamos con el de Numpy, que es numérico:

>>> sp.pi**2
pi**2

>>> sp.pi.evalf()
3.14159265358979

>>> (sp.pi + sp.exp(1)).evalf()
5.85987448204884

como se ve, sin embargo, se puede usar el método evalf() para evaluar una expresión para tener un valor en punto flotante (float).

Para hacer operaciones simbólicas hay que definir explícitamente los símbolos que vamos a usar, que serán en general las variables y otros elementos de nuestras ecuaciones:

>>> x = sp.Symbol('x')
>>> y = sp.Symbol('y')

Y ahora ya podemos manipularlos como queramos:

>>> x+y+x-y
2*x

>>> (x+y)**2
(x + y)**2

>>> ((x+y)**2).expand()
2*x*y + x**2 + y**2

Es posible hacer una substitución de variables usando subs(viejo, nuevo):

>>> ((x+y)**2).subs(x, 1)
(1 + y)**2

>>> ((x+y)**2).subs(x, y)
4*y**2

Operaciones algebraicas

Podemos usar apart(expr, x) para hacer una descomposición parcial de fracciones:

>>>  1/( (x+2)*(x+1) )
       1
───────────────
(2 + x)*(1 + x)

>>>  sp.apart(1/( (x+2)*(x+1) ), x)
  1          1
─────── - ────────
1 + x      2 + x

>>>  (x+1)/(x-1)
-(1 + x)
─────────
 1 - x

>>>  sp.apart((x+1)/(x-1), x)
      2
1 - ──────
    1 - x

Cálculo de límites

Sympy puede calcular límites usando la función limit() con la siguiente sintaxis: limit(función, variable, punto), lo que calcularía el limite de \(f(x)\) cuando variable -> punto:

.. code:: ipython3
>>> x = sp.Symbol("x")
>>> sp.limit(sin(x)/x, x, 0)
1

es posible incluso usar límites infinitos:

>>> sp.limit(x, x, oo)
oo

>>> sp.limit(1/x, x, oo)
0

>>> sp.limit(x**x, x, 0)
1

Cálculo de derivadas

La función de Sympy para calcular la derivada de cualquier función es diff(func, var). Veamos algunos ejemplos:

>>> x = sp.Symbol('x')
>>> diff(sp.sin(x), x)
cos(x)
>>> diff(sp.sin(2*x), x)
2*cos(2*x)

>>> diff(sp.tan(x), x)
1 + tan(x)**2

Se puede comprobar que es correcto calculando el límite:

>>> dx = sp.Symbol('dx')
>>> sp.limit( (tan(x+dx)-tan(x) )/dx, dx, 0)
1 + tan(x)**2

También se pueden calcular derivadas de orden superior indicando el orden de la derivada como un tercer parámetro opcional de la función diff():

>>> sp.diff(sin(2*x), x, 1)         # Derivada de orden 1
2*cos(2*x)

>>> sp.diff(sin(2*x), x, 2)         # Derivada de orden 2
-4*sin(2*x)

>>> sp.diff(sin(2*x), x, 3)         # Derivada de orden 3
-8*cos(2*x)

Expansión de series

Para la expansión de series se aplica el método series(var, punto, orden) a la serie que se desea expandir:

>>> cos(x).series(x, 0, 10)
1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)
>>> (1/cos(x)).series(x, 0, 10)
1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + O(x**10)


e = 1/(x + y)
s = e.series(x, 0, 5)

pprint(s)

La función pprint de Sympy imprime el resultado de manera más legible:

     4    3    2
1   x    x    x    x
─ + ── - ── + ── - ── + O(x**5)
y    5    4    3    2
    y    y    y    y

Si usamos una qtconsole (iniciada con ipython qtconsole) e ejecutamos init_printing(), las fórmulas se mostrarán con \(\LaTeX\), si está instalado.

Integración simbólica

La integración definida e indefinida de funciones es una de las funcionalidades más interesantes de :mod:sympy`. Veamos algunos ejemplos:

>>> sp.integrate(6*x**5, x)
x**6
>>> sp.integrate(sp.sin(x), x)
-cos(x)
>>> sp.integrate(sp.log(x), x)
-x + x*log(x)
>>> sp.integrate(2*x + sp.sinh(x), x)
cosh(x) + x**2

También con funciones especiales:

>>> sp.integrate(exp(-x**2)*erf(x), x)
pi**(1/2)*erf(x)**2/4

También es posible calcular integrales definidas:

>>> sp.integrate(x**3, (x, -1, 1))
0
>>> sp.integrate(sin(x), (x, 0, pi/2))
1
>>> sp.integrate(cos(x), (x, -pi/2, pi/2))
2

Y también integrales impropias:

>>> sp.integrate(exp(-x), (x, 0, oo))
1
>>> sp.integrate(log(x), (x, 0, 1))
-1

Algunas integrales definidas complejas es necesario definirlas como objeto Integral() y luego evaluarlas con el método evalf():

>>> integ = sp.Integral(sin(x)**2/x**2, (x, 0, oo))
>>> integ.evalf()
>>> 1.6

Ecuaciones algebraicas y álgebra lineal

También se pueden resolver sistemas de ecuaciones de manera simbólica:

>>> # Una ecuación, resolver x
>>> sp.solve(x**4 - 1, x)
[I, 1, -1, -I]

# Sistema de dos ecuaciones. Resuelve x e y
>>>  sp.solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
{y: 1, x: -3}

Sympy tiene su propio tipo de dato Matriz, independiente del de Numpy/Scipy. Con él se pueden definir matrices numéricas o simbólicas y operar con ellas:

>>> # Matriz identidad 2x2
>>> sp.Matrix([[1,0], [0,1]])
[1, 0]
[0, 1]

>>> x = sp.Symbol('x')
>>> y = sp.Symbol('y')
>>> A = sp.Matrix([[1,x], [y,1]])
>>> A
[1, x]
[y, 1]

>>> A**2
[1 + x*y,     2*x]
[    2*y, 1 + x*y]

Hay que fijarse en que muchas de las funciones anteriores ya existen en Numpy con el mismo nombre (ones(), eye(), etc.), por lo que si queremos usar ambas debemos importar los paquetes con otro nombre, Por ejemplo:

>>> import numpy as np

>>> # Funcion eye de Sympy (matriz)
>>> sp.eye(3)
[1, 0, 0]
[0, 1, 0]
[0, 0, 1]

>>> # Funcion eye de Numpy (array)
>>> np.eye(3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

Es posible operar entre ellas, salvo que las matrices de Numpy no pueden operar con símbolos, algo que se puede hacer con Sympy. La selección de elementos de matrices de Sympy se hace de manera idéntica a los arrays o matrices de Numpy:

>>> # Multiplico la matriz identidad por x
>>> x = sp.Symbol('x')
>>> M = sp.eye(3) * x
>>> M
[x, 0, 0]
[0, x, 0]
[0, 0, x]

>>> # Substituyo x por 4 en la matriz
>>> M.subs(x, 4)
[4, 0, 0]
[0, 4, 0]
[0, 0, 4]

>>> # Substituyo la variable x por y en la matriz
>>> y = sp.Symbol('y')
>>> M.subs(x, y)
[y, 0, 0]
[0, y, 0]
[0, 0, y]

>>>  def f(x): return 1.5*x**2
....:

>>>  sp.eye(3).applyfunc(f)
[1.5,   0,   0]
[  0, 1.5,   0]
[  0,   0, 1.5]


>>> M = sp.Matrix(( [1, 2, 3], [3, 6, 2], [2, 0, 1] ))

>>> # Determinante de la matriz
>>> M.det()
-28
# Matriz inversa
>>>  M.inv()
[-3/14, 1/14,  1/2]
[-1/28, 5/28, -1/4]
[  3/7, -1/7,    0]

>>> # Substituyo algunos elementos de la matriz
>>> M[1,2] = x
>>> M[2,1] = 2*x
>>> M[1,1] = sqrt(x)
>>> M
[1,        2, 3]
[3, x**(1/2), x]
[2,      2*x, 1]

Podemos resolver un sistema de ecuaciones por el método LU:

>>> # Matriz 3x3
>>> A = sp.Matrix([ [2, 3, 5], [3, 6, 2], [8, 3, 6] ])
>>> # Matriz 1x3
>>> x = sp.Matrix(3,1,[3,7,5])
>>> b = A*x
>>> # Resuelvo el sistema por LU
>>> soln = A.LUsolve(b)
>>> soln
[3]
[7]
[5]

Gráficos con Sympy

Sympy trae su propio módulo para hacer gráficos, que en realidad utiliza matplotlib, pero con la ventaja de que no es necesario darle valores numéricos para representar las funciones, si bien se pueden indicar límites.

sp.plot(sp.sin(x)*sp.log(x))

sp.plot(sp.sin(2*sp.sin(2*sp.sin(x))))

Exportando a \(\LaTeX\)

Cuando hemos terminado los cálculos, podemos pasar a \(\LaTeX\) la ecuación que queremos:

In [5]: int1 = sp.Integral( sp.sin(x)**3 /(2*x), x)

In [6]: sp.latex(int1)
Out[6]: '\\int \\frac{\\sin^{3}{\\left (x \\right )}}{2 x}\\, dx'