Computación simbólica con Sympy
Contents
Computación simbólica con Sympy#
Important
Para tener presente
El paquete sympy nos brinda la posibilidad de usar expresiones aproximadas de manera simbólica privilegiando el símbolo frente a la evaluación.
Computación simbólica#
Hasta el momento hemos trabajado con varios elementos matemáticos que en ocasiones se toman de manera aproximada, por ejemplo, la raíz cuadrada de un número entero. Aunque hay algunos números cuya raíz es entera la mayoría tiene raices cuadradas irracionales y la máquina nos ofrece apenas una aproximación a esos valores.
import numpy as np
np.sqrt(8)
2.8284271247461903
np.sqrt(8)*np.sqrt(8)
8.000000000000002
import sympy
sympy.sqrt(8)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 import sympy
2 sympy.sqrt(8)
ModuleNotFoundError: No module named 'sympy'
sympy.sqrt(8)*sympy.sqrt(8)
Nota: para obtener una salida estética de sympy usamos la orden:
from sympy import init_printing
init_printing()
El poder de la computación simbólica
El poder real de un sistema de cálculo simbólico como SymPy es la capacidad de hacer todo tipo de cálculos simbólicamente. SymPy puede simplificar expresiones, calcular derivadas, integrales y límites, resolver ecuaciones, trabajar con matrices y mucho, mucho más, y hacerlo todo simbólicamente. Incluye módulos para trazar, imprimir (como salida impresa en 2D de fórmulas matemáticas, LATEX), generación de código, física, estadística, combinatoria, teoría de números, geometría, lógica y más.
Conceptos básicos#
Para iniciar el trabajo con la librería sympy debemos aclarar rutinas que son fundamentales para iniciar cualquier trabajo con este paquete. En primer lugar, por ejemplo, hay que aclarar que vamos a usar un símbolo así mismo, las ecuaciones tienen que poder expresarse y no podemos usar ni = ni ==. Veamos algunos de estos conceptos básicos.
Símbolos#
Para usar variables simbólicas debemos definirlas, usamos la función symbols:
import sympy as sp
x = sp.symbols("x")
x
Notemos que las operaciones las realiza de manera análoga a lo que solemos hacer con lápiz y papel:
(x+1)**2
sp.expand((x+1)**2)
(3*x-5)*(2*x+7)-x
x+2*x-8+5*x+7
comoescribo = sp.symbols("comoveo")
comoescribo+8
Sustitución#
Para remplazar valores en una expresión usamos subs
expr = (3*x-5)*(2*x+7)-x
expr
expr.subs(x,1)
expr.subs(x,0)
También se pueden sustituir expresiones
y=sp.symbols("y")
expr = x**y
expr
expr = expr.subs(y, x**y)
expr
from sympy import *
expr = sin(2*x)
expand_trig(expr)
expr.subs(sin(2*x), 2*sin(x)*cos(x)).subs(cos(2*x),2*cos(x)**2-1)
z = symbols("z")
expr = x**3 + 4*x*y - z
expr
expr.subs([(x, 2), (y, 4), (z, 0)])
Convertir cadenas a expresiones SymPy#
Usamos sympify para convertir cadenas en expresiones:
str_expr = "x**2 + 3*x - 1/2"
print(str_expr)
expr = sympify(str_expr)
expr
x**2 + 3*x - 1/2
import ipywidgets as widgets
from ipywidgets import interact
def modulo1(t,n):
z="Error"
n=sympify(n)
try:
t=sympify(t)
z=sp.expand(t**n)
except:
print('Error!!')
return z
widn=widgets.IntSlider(min=0,max=10,value=1)
interact(modulo1, t="sin(x)",n=widn)
<function __main__.modulo1(t, n)>
Obtener el valor en punto flotante#
Usamos evalf para obtener el valor en punto flotante.
expr = sqrt(8)
expr
expr.evalf()
subs y evalf son buenos si desea hacer una evaluación simple, pero si tiene la intención de evaluar una expresión en muchos puntos, hay formas más eficientes. Por ejemplo, si desea evaluar una expresión en mil puntos, usar SymPy sería mucho más lento de lo necesario, especialmente si solo le importa la precisión de la máquina. En su lugar, debe usar librerías como NumPy y SciPy .
La forma más fácil de convertir una expresión SymPy en una expresión que puede evaluarse numéricamente es usar la lambdifyfunción. lambdify actúa como una función lambda, excepto que convierte los nombres de SymPy a los nombres de la biblioteca numérica dada, generalmente NumPy. Por ejemplo:
import numpy as np
a = np.arange(10)
a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
expr = sin(x)
f = lambdify(x, expr, "numpy")
f(a)
array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ,
-0.95892427, -0.2794155 , 0.6569866 , 0.98935825, 0.41211849])
f = lambdify(x, expr, "math")
f(5)
-0.9589242746631385
Ecuaciones#
Como sabemos = y == ya fueron usado en el lenguaje para unas operaciones específicas.
x+1==4
False
sp.Eq(x + 1, 4)
Igualdad de expresiones y simplificación#
Ahora saltemos y hagamos algunas matemáticas interesantes. Una de las características más útiles de un sistema de manipulación simbólica es la capacidad de simplificar expresiones matemáticas.
SymPy tiene docenas de funciones para realizar varios tipos de simplificación, pero hay una función general llamada simplify() que intenta aplicar todas estas funciones de manera inteligente para llegar a la forma más simple de una expresión. Aquí hay unos ejemplos:
Es importante tener esta definición en cuenta cuando se presentan igualdades entre expresiones.
simplify(sin(x)**2 + cos(x)**2)
(x**3 + x**2 - x - 1)/(x**2 + 2*x + 1)
simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
simplify(gamma(x)/gamma(x - 2))
simplify (x**2+2*x+1 - (x+1)**2)
Un problema de simplify() es que puede ser innecesariamente lento, ya que intenta muchos tipos de simplificaciones antes de elegir la mejor. Si ya sabe exactamente qué tipo de simplificación está buscando, es mejor aplicar las funciones de simplificación específicas que aplican esas simplificaciones.
Simplificación de funciones polinómicas racionales#
expand() es una de las funciones de simplificación más comunes en SymPy. Aunque tiene muchos ámbitos, por ahora, consideraremos su función en la expansión de expresiones polinómicas. Por ejemplo:
expand((x+2)**3)
factor() toma un polinomio y lo factoriza en factores irreducibles sobre los números racionales. Por ejemplo:
factor(x**3-x**2+x-1)
factor_list(x**3-x**2+x-1)
(1, [(x - 1, 1), (x**2 + 1, 1)])
cancel() tomará cualquier función racional y la pondrá en la forma canónica estándar, 𝑝/𝑞, dónde 𝑝 y 𝑞 son polinomios expandidos sin factores comunes, y los coeficientes principales de 𝑝 y 𝑞 no tienen denominadores (es decir, son enteros).
cancel((x**2 + 2*x + 1)/(x**2 + x))
expr = (x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)
expr
cancel(expr)
apart() realiza una descomposición en fracciones parciales de una función racional.
expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)
expr
apart(expr)
Para simplificar expresiones usando identidades trigonométricas, es recomendable el uso de trigsimp().
trigsimp(sin(x)*tan(x)/sec(x))
cancel(sin(x)*tan(x)/sec(x))
Para expandir las funciones trigonométricas, es decir, aplicar la suma o las identidades de doble ángulo, usamos expand_trig().
expand_trig(sin(x + y))
expand_trig(tan(2*x))
powsimp() permite simplificaciones con leyes de exponentes.
x, y = symbols('x y', positive=True)
a, b = symbols('a b', real=True)
z, t, c = symbols('z t c')
powsimp(x**a*x**b)
powsimp(z**t*z**c)
Con rewrite() podemos reescribir algunas funciones especiales en términos de otra:
tan(x).rewrite(sin)
factorial(x).rewrite(gamma)
exp(I*x).rewrite(sin)
Fracciones Continuas
Una fracción continua es una expresión de la forma:
Realizar una función en sympy que admita una lista \([a_0,a_1,a_2,\cdots,a_n]\) y que retorne el valor obtenido por la fracción.
L=[1,2,3]
L.reverse()
L
[3, 2, 1]
def lista_a_fraccion(L):
L.reverse()
expr=Integer(0)
for i in L:
expr=1/(expr+i)
return expr
lista_a_fraccion([x,y,z])
lista_a_fraccion([1,2*x,3])
Cálculo#
SymPy también permite realizar tareas básicas del cálculo diferencial e integral.
from sympy import *
x, y, z = symbols('x y z')
Derivadas#
diff(cos(x),x) #derivadas en una sola variable
diff(exp(4*x**2+5*x),x) #derivadas en una sola variable
diff(y*exp(4*x**2+5*y),x) #derivadas parciales
diff(cos(x),x,x) #derivadas de orden superior
diff(cos(x),x,x,x)
diff(cos(x),x,3)
diff(y*exp(4*x**2+5*y),x,y) #derivadas de orden superior
Como en las ecuaciones, a veces requerimos derivadas expresadas, sin evaluarse:
deriv = Derivative( exp (x * y * z), x, y, y, z)
deriv
Para evaluarla:
deriv.doit()
diff( exp (x * y * z), x, y, y, z)
Integrales#
integrate(cos(x),x) # Integrales
integrate(exp (-x), (x, 0, oo))
Integral(exp (-x), (x, 0, oo))
Integral(exp (-x), (x, 0, oo)).doit()
integrate(x**2, (x, 0,1))
Integral(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))
integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))
Integral(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo)).doit()
Y ¿si la integral no es posible de calcular?
integrate(x**x, x)
Para expresar una integral no evaluada usamos Integral:
integ = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x -exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)), x)
integ
integ.doit()
Límites#
limit(sin(x)/x, x, 0) #Límites
expr = sin(x)/x
expr.subs(x, 0)
limit(expr, x, 0)
expr = Limit((cos(x) - 1)/x, x, 0,"-") #Límites sin evaluar
expr
expr.doit()
expr = exp(sin(x))
expr.series(x,0, 10) # Series
exp(x - 6).series(x, 6,10) # Series no centradas en 0
Solucionar ecuaciones#
Ahora hagamos una rápida exploración por las diferentes ecuaciones que se resuelven en Sympy.
Eq(x**2,1)
solveset(Eq(x**2,1),x)
solveset(Eq(x**2,1),x)[0]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-94-bf4836c147c4> in <module>
----> 1 solveset(Eq(x**2,1),x)[0]
TypeError: 'FiniteSet' object is not subscriptable
list(solveset(Eq(x**2,1),x))[0]
solveset(x ** 2 - x, x) # Raices
solveset(x - x, x)
solveset(x - x, x,domain=sp.Integers)
solveset(x - x, x,domain=Interval(0,10))
solveset(x - x, x,domain=Interval.Ropen(0,10))
solveset(x - x, x,domain=Interval.Lopen(0,10))
solveset(x - x, x,domain=Interval.open(0,10))
solveset(x - x, x,domain=FiniteSet(0,1,2,3))
solveset(sin(x) - 1, x, domain=sp.Reals)
solveset(sin(x) - 1, x, domain=Interval.open(5,34))
solveset(exp(x), x) # Cuando la solución no existe
solveset(cos(x) - x, x) # Cuando es incapaz de encontrar una solución
Sistemas de ecuaciones lineales#
Usamos linsolve para encontrar las soluciones:
linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))
linsolve([x + y + z - 1, x - y + z - 3,x+y+1 ], (x, y, z))
linsolve([x + y + z - 1, x + y + z - 3,x+y+1 ], (x, y, z))
Sistemas no lineales#
Usaremos nonlinsolve
a, b, c, d = symbols('a, b, c, d', real=True)
nonlinsolve([a**2 + a, a - b], [a, b])
nonlinsolve([x*y - 1, x - 2], x, y)