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.

.. code:: ipython3

    >>> 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 :math:`e` o :math:`\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:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

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

Y ahora ya podemos manipularlos como queramos:

.. code:: ipython3

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

.. code:: ipython3

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

.. code:: ipython3

    >>>  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 :math:`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:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

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

.. code:: ipython3

    >>> 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:

.. code:: ipython3

    >>> 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 :math:`\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:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

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


También es posible calcular integrales definidas:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

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

.. code:: ipython3

    >>> 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:

.. code:: ipython3

    >>> # 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:

.. code:: ipython3

    >>> # 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:

.. code:: ipython3

    >>> 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:

.. code:: ipython3

    >>> # 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:

.. code:: ipython3

    >>> # 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.

.. code:: ipython3

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


Exportando a :math:`\LaTeX`
---------------------------

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

.. code:: ipython3

    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'