15 Álgebra Simbólica

Julia permite realizar operaciones y cálculos simbólicos llamando, a través de PyCall, a la librería SymPy de matemática simbólica de Python. El paquete que realiza el puente entre Julia y Python es sympy. Para utilizar esta funcionalidad el primer paso es instalar dicho paquete:

#using Pkg
#Pkg.add("SymPy");

Para usar el paquete sympy tanto Python como SymPy deben estar instalados en el computador. Existen dos formas de instalar SymPy. La primera es instalar el paquete PyCall en Julia mediante Pkg.add("PyCall"), el mismo que instara los paquetes requeridos de Python vía Conda. La segunda alternativa, es instalar Anaconda en el computador que es un paquete de Python y varias librerias pesadas para su aplicación en las ciencias.

using SymPy;

15.1 Declarar variables simbólicas

Las variables simbólicas, a diferencia de las variables típicamente utilizadas en Julia, no son evaluadas inmediatamente. En lugar de ello, se utiliza el “símbolo” para referirse a ellas. Existen dos formas equivalentes de declarar variables simbólicas. La primera es usar las funciones Sym() o symbols(), mientras que la segunda usa la macro @vars o @syms.

x,y,z = Sym("x,y,z")

(x, y, z)

@vars a b c

(a, b, c)

Al declarar las variables es posible realizar ciertos supuestos respecto de su comportamiento. El formato de declaración se realiza mediante el uso de la opción supuesto = true de la función symbols (en este caso no usar Sym). La declaración también funciona con la macro @vars. cuando se declara un suspuesto se impone una restricción respecto de los valores que puede tomar una variable. Una lista, no exaustiva per útil, de dichos supuestos es la siguiente: real, odd, even, zero, nonzero, positive, negative, finite, prime, etc. Ver la lista completa de supuestos en el siguiente link. Algunos ejemplos:

x,y,z = symbols("x,y,z", real=true)  # x,y,z son reales

(x, y, z)

@syms a b c positive=true # a,b,c son todos positivos

(a, b, c)

15.2 Expresiones y ecuaciones simbólicas

Asignar expresiones matemáticas

Una vez declaradas las variables es posible combinarlas para crear ecuaciones. Al igual que antes se usa = para asignar (note que =en este caso es asignación y no igualdad en referencia a dos lados de una ecuación). Es posible utilizar la operaciones algebráicas +,´,*,\,^ asi como funciones, como sin(), cos(), log(), exp() entre otras, para crear expresiones.

f = 2*x + 3*y^2 + log(z) # Asignamos la expresión a la variable f

\[\begin{equation*}2 x + 3 y^{2} + \log{\left(z \right)}\end{equation*}\]

Además es posible crear expresiones a partir de otras expresiones simbólicas.

g = x*f

\[\begin{equation*}x \left(2 x + 3 y^{2} + \log{\left(z \right)}\right)\end{equation*}\]

s = exp(x+y)

\[\begin{equation*}e^{x + y}\end{equation*}\]

Sustitución

Sympy permite substituir variables simbólica en expresiones mediante la función subs(expresion, (subs)) o simplemente evaluando la expresión expresion(subs). La substitución se puede realizar entre símbolos, expresiones, o números.

fx = 2*x^2

\[\begin{equation*}2 x^{2}\end{equation*}\]

Por ejemplo, substituímos x por y en la expresión fx y luego asignamos la nueva expresión a fy.

fy = subs(fx, (x, y))

\[\begin{equation*}2 y^{2}\end{equation*}\]

Podemos crear a expresión g = 3*y+z y substituir toda esta expresión por x en fx:

g = 3*y+z
ff = subs(fx, (x, g))

\[\begin{equation*}2 \left(3 y + z\right)^{2}\end{equation*}\]

Alternativamente, para producir el mismo resultados podemos usar simplemente la sintaxis de una función:

ff = fx(g)

\[\begin{equation*}2 \left(3 y + z\right)^{2}\end{equation*}\]

Para sustitur más de una variable en una expresión se usa subs(expresion, (subs1), (subs2),...):

fxy = 2*x + 3*y^2
fz = subs(fxy, (x,z),(y,z^2))  # sustitución: z por x y z^2 por y  

\[\begin{equation*}3 z^{4} + 2 z\end{equation*}\]

Alternativamente se puede evaluar la expresión directamente. Es importante para esto tener en mente el orden de las variables. Para determinar dicho orden usar free_symbols(expresion).

fz = fxy(z,z^2)

\[\begin{equation*}3 z^{4} + 2 z\end{equation*}\]

Otra alternativa es usar => para asignar, por ejemplo:

fz = fxy(x=>z,y=>z^2)

\[\begin{equation*}3 z^{4} + 2 z\end{equation*}\]

En cualquier caso, es recomendable definir explícitamente la variable a reemplazar para evitar ambiguedades.

Evaluación numérica

Reemplazar variables por valores numéricos permite evaluar numéricamente la expresión. Por ejemplo:

fxy_num = subs(fxy, (x,2.0),(y,2.0))

\[\begin{equation*}16.0\end{equation*}\]

Este número es un objeto del typo sym (simbólico) y no un tipo numérico con el cual trabajamos usualmente en Julia, como Int64 o Float64. Esto implica que si pasamos este número por una función que no acepta o no trabaja con objetos tipo sym obtendremos un error.

typeof(fxy_num)

Sym

Para transformar la evaluación numérica de sympy usamos la función N(expresión).

N(fxy_num)

16.0

typeof(N(fxy_num))

Float64

Ecuaciones

Existen dos formas de declarar ecuaciones. Primero, la función Eq() permite declarar ecuaciones bajo el formato Eq(lado_izquierdo,lado_derecho). Por ejemplo:

eq = Eq(x^3+y,2)

\[\begin{equation*}x^{3} + y = 2\end{equation*}\]

Segundo, es posible utilizar operadores tipo unicode para diferenciar la declaración de una ecuación con respecto a la asignación de una variable. Los operadores son: \ll[tab] (menor), \leqq[tab] (menor o igual), \Equal[tab] (igual), \geqq[tab] (mayor o igual), \gg[tab] (mayor), y \neg[tab] (negación). Por ejemplo:

eq = x^3+y ⩵ 2 

\[\begin{equation*}x^{3} + y = 2\end{equation*}\]

15.3 Manipulando Polinomios

Las funciones expand(), factor(), collect(), y simplify()

La función expand(expresión) expande un polinomio como la suma de monomios (realizando en el proceso todas las simplificanciones necesarias). Por ejemplo:

fxy = (x+y)^2

\[\begin{equation*}\left(x + y\right)^{2}\end{equation*}\]

fxy_exp = expand(fxy)

\[\begin{equation*}x^{2} + 2 x y + y^{2}\end{equation*}\]

La función factor(expresión) toma un polinomio y lo factoriza en sus factores irreductibles. La función factor() es el opuesto de la función expand(). Por ejemplo tomemos el resultado anterior y apliquemos la función factor():

factor(fxy_exp)

\[\begin{equation*}\left(x + y\right)^{2}\end{equation*}\]

La función collect(expresión) agrupa las potencias comúnes de un término presentes en una expresión. Por ejemplo, agrupemos los términos comúnes de \(x\) en la expresión \(xy + xy^2 + x^2y + x\):

fxy = x*y + x*y^2 + x^2*y + x
collect(fxy, x)

\[\begin{equation*}x^{2} y + x \left(y^{2} + y + 1\right)\end{equation*}\]

Tomemos ahora los términos comunes de \(y\) de la misma expresión:

collect(fxy, y)

\[\begin{equation*}x y^{2} + x + y \left(x^{2} + x\right)\end{equation*}\]

En algunos contextos puede ser útil extraer el coeficiente asociado a un término común. Por ejemplo, en la anterior expresión podemos extraer el coeficiente que acompaña a \(y^2\) usando la notación expresión.coeff(término_común):

fxy_n = collect(fxy, y)
fxy_n.coeff(y^2)

\[\begin{equation*}x\end{equation*}\]

La función simplify(expresión) simplifica una expresión reduciendo el número de términos al mínimo posible.

fxy = 2*x^2 - x^2 - x*y
simplify(fxy)

\[\begin{equation*}x \left(x - y\right)\end{equation*}\]

La función apart(expresión) realiza la descompocición en fracciones parciales de una función racional. Por ejemplo:

fx_1 = x^2+x
fx_2 = x^2+2x+1
ratio = fx_1/fx_2

\[\begin{equation*}\frac{x^{2} + x}{x^{2} + 2 x + 1}\end{equation*}\]

apart(ratio)

\[\begin{equation*}1 - \frac{1}{x + 1}\end{equation*}\]

La función together(expresión) realiza la operación inversa a apart(expresión). Por ejemplo:

fx = 1/x^2 + 1/x^3

\[\begin{equation*}\frac{1}{x^{2}} + \frac{1}{x^{3}}\end{equation*}\]

together(fx)

\[\begin{equation*}\frac{x + 1}{x^{3}}\end{equation*}\]

La función cancel(expresión) toma una función racional y la expresa en su forma canónica estándar \(p/q\), con \(p\) y \(q\) polinomios sin factores comunes. Por ejemplo:

num = expand((x-1)*(x-2)*(x-3))
denom = expand((x-1)*(x-4))
fx = num/denom

\[\begin{equation*}\frac{x^{3} - 6 x^{2} + 11 x - 6}{x^{2} - 5 x + 4}\end{equation*}\]

cancel(fx)

\[\begin{equation*}\frac{x^{2} - 5 x + 6}{x - 4}\end{equation*}\]

Tambien es posible realizar operaciones y simplificaciones con potencias. Por ejemplo:

fx = (x^a)*(x^b)

\[\begin{equation*}x^{a} x^{b}\end{equation*}\]

sympy.powsimp(fx)

\[\begin{equation*}x^{a + b}\end{equation*}\]

Note que usamos las notaciones sympy.powsimp(expresión). Esto se debe a que la función powersim no ha sido exportada al entorno global en el paquete sympy, por lo que tenemos que indicarle a Julia que dicha función viene de ese paquete.

15.4 Cálculo

El paquete sympy es capaz de calcular derivadas, integrales y límites. La función diff(expresión,variable_a_derivar) toma la derivada de la función declarada en expresióncon respecto a la variable variable_a_derivar. Si la función tiene un solo argumento, no es necesario declarar la variable a derivar. Por ejemplo:

fx = cos(x)*x^2
diff(fx)

\[\begin{equation*}- x^{2} \sin{\left(x \right)} + 2 x \cos{\left(x \right)}\end{equation*}\]

Si la función tiene más de una variable, definimos explícitamente la variable con respecto a la cual queremos tomar la derivada parcial. Por ejemplo:

fxy = x*y + x*y^2 + x^2*y + x

# Parcial con respecto a x
diff(fxy,x)

\[\begin{equation*}2 x y + y^{2} + y + 1\end{equation*}\]

# Parcial con respecto a y
diff(fxy,y)

\[\begin{equation*}x^{2} + 2 x y + x\end{equation*}\]

Alternativamente, si definimos un vector con todas las variables a derivar, pordemos calcular el vector gradiente de una función:

grad = diff(fxy,[x, y])

\[\begin{equation*}\left[\begin{matrix}2 x y + y^{2} + y + 1\\x^{2} + 2 x y + x\end{matrix}\right]\end{equation*}\]

Nota: El tipo de objeto obtenido al calcular el vector gradiente es SymMatrix que es específico a sympy. Como la funcionalidad no ha sido trasladada al 100% entre Julia y sympy de Python, la forma de acceder en este caso a los elementos del vector es indexando sus elementos partiendo en 0 (en Python el primer elemento de un vector es el que tiene posición 0):

@show typeof(grad)
@show grad[0]
@show grad[1];

typeof(grad) = SymMatrix
grad[0] = 2*x*y + y^2 + y + 1
grad[1] = x^2 + 2*x*y + x

Para calcular derivada de orden superior usamos diff(expresión,variable_a_derivar,orden_derivada).

diff(fxy,x,2) # segunda derivada con respecto a x

\[\begin{equation*}2 y\end{equation*}\]

diff(fxy,x,3) # tercera derivada con respecto a x

\[\begin{equation*}0\end{equation*}\]

Para calcular derivadas cruzadas usamos diff(expresión,var_primera_der,var_segunda_der,...).

diff(fxy,x,y)

\[\begin{equation*}2 x + 2 y + 1\end{equation*}\]

El paquete sympy provee funcionalidad para calcular integrales tanto definidas como idefinidas. Para calcular una integral indefinida usamos integrate(expresión,variable_a_integrar). Por ejemplo:

integrate(cos(x), x)

\[\begin{equation*}\sin{\left(x \right)}\end{equation*}\]

Si la función tiene más de una variable, definimos la dimensión en la cual queremos integrar:

fxy = x*y + x*y^2 + x^2*y + x
integrate(fxy,x)

\[\begin{equation*}\frac{x^{3} y}{3} + x^{2} \left(\frac{y^{2}}{2} + \frac{y}{2} + \frac{1}{2}\right)\end{equation*}\]

Para integrales dobles usamos integrate(expresión,(var_a_integrar_interior,),(var_a_integrar_exterior,)). Por ejemplo por ejemplo para integrar:

\[\int \int (xy + xy^2 + x^2*y + x) dx dy\]

integrate(fxy,(x,),(y,))  # Note la notación tipo tuplas.

\[\begin{equation*}\frac{x^{2} y^{3}}{6} + \frac{x^{2} y}{2} + y^{2} \left(\frac{x^{3}}{6} + \frac{x^{2}}{4}\right)\end{equation*}\]

Para calcular una integral definida usamos integrate(expresión, (variable_a_integrar, limite_inf, limite_sup)). Note que la variable a integrar y sus límites son provisto en forma de tupla. Por ejemplo:

fx = exp(-x)
integrate(fx, (x, 0, 1))

\[\begin{equation*}1 - e^{-1}\end{equation*}\]

Los límites \(-\infty\) e \(\infty\) sin permitidos y son declarados con doble o: oo.

integrate(fx, (x, -oo, oo))

\[\begin{equation*}\infty\end{equation*}\]

Para calcular una integral definida doble usamos integrate(expresión, (variable_interior, limite_inf, limite_sup), (variable_exterior, limite_inf, limite_sup)). Por ejemplo:

integrate(exp(-x^2 - y^2), (x, -oo, oo), (y, -oo, oo))

\[\begin{equation*}\pi\end{equation*}\]

Finalmente, para calcular límites usamos limit(expresión, variable, punto_convergencia). Por ejemplo:

limit(sin(x)/x, x, 0)

\[\begin{equation*}1\end{equation*}\]

15.5 Resolviendo de Ecuaciones y Sistemas de Ecuaciones

Para resolver ecuaciones, lineales o no lineales, se usa la función solve(ecuación, variable) sobre una ecuación generada con Eq(lado_izq,lado_der). Por ejemplo, para resolver por \(x\) la ecuación \(x^2+y^2=2\) usamos:

fxy = Eq(x^2+y^2,2)
solucion = solve(fxy, x)

\[ \left[ \begin{array}{r}- \sqrt{2 - y^{2}}\\\sqrt{2 - y^{2}}\end{array} \right] \]

En este caso tenemos dos soluciones, por lo que solve nos entrega un arreglo con dos elementos los mismos que pueden ser recuperados usando su respectivo índice:

@show solucion[1]
@show solucion[2];

solucion[1] = -sqrt(2 - y^2)
solucion[2] = sqrt(2 - y^2)

En le caso de ecuaciones lineales, tenemos la misma sintaxis:

fxy = Eq(4*x+2*y,2)
solve(fxy, x)   # note que el resultado es una arreglo con 1 elemento

\[ \left[ \begin{array}{r}\frac{1}{2} - \frac{y}{2}\end{array} \right] \]

Si la ecuación no es declarada y solamente se pasa una expresión a solve, éste asume que la ecuación es expresión==0 y resuelve un problema de raíces.

fx = x^2-4
solve(fx, x)

\[ \left[ \begin{array}{r}-2\\2\end{array} \right] \]

Para el caso de sistemas de ecuaciones, lineales o no lineales, usamos la misma sintaxis pero declarando el sistema en un arreglo, esto es solve([ecu1, ..., ecuN], [var1,...,varN]):

f1 = Eq(a*x + b*y-3,0)
f2 = Eq(c*x + b*y - 1,0)
solucion = solve([f1, f2], [x,y])

Dict{Any,Any} with 2 entries:
  y => (a - 3*c)/(b*(a - c))
  x => 2/(a - c)

Cómo la solución es única para \(x\) e \(y\), sympy nos entrega un diccionario. Podemos acceder a cada una de las soluciones utilizando el identificador, en este caso el nombre de la variable.

solucion[x]

\[\begin{equation*}\frac{2}{a - c}\end{equation*}\]

solucion[y]

\[\begin{equation*}\frac{a - 3 c}{b \left(a - c\right)}\end{equation*}\]

Para sistema de ecuaciones no lineales usamos exactamente la misma sintaxis. La única diferencia es que, como en general los sistemas no lineales podrían tener más de una solución, el resultado es almacenado en tuplas dentro un arreglo. Por ejemplo:

f1nl = Eq(x^2 + x,a)
f2nl = Eq(x - y,b)
solucion = solve([f1nl, f2nl], [x, y])

2-element Array{Tuple{Sym,Sym},1}:
 (-sqrt(4*a + 1)/2 - 1/2, -b - sqrt(4*a + 1)/2 - 1/2)
 (sqrt(4*a + 1)/2 - 1/2, -b + sqrt(4*a + 1)/2 - 1/2)

En este caso obtenemos las soluciones individuales indexando primero el elemento de arreglo y luego el elemento de la tupla. Por ejemplo, las dos soluciones de \(x\) están en el primer elemento del arreglo y luego en los índices 1 y 2 de la tupla:

solucion[1][1]

\[\begin{equation*}- \frac{\sqrt{4 a + 1}}{2} - \frac{1}{2}\end{equation*}\]

solucion[1][2]

\[\begin{equation*}- b - \frac{\sqrt{4 a + 1}}{2} - \frac{1}{2}\end{equation*}\]

15.6 Matrices

Para declarar matrices simbólicas usamos arreglos en dos dimensiones, como lo hacíamos con matrices numéricas, pero cuyos elementos son simbólicos. Esto es, una matriz se declara matriz = [] y las filas y columnas se determinan con ; y (espacio), respectivamente.

@vars a b c d x y z w;

A = [a b; c d]

\[\left[ \begin{array}{rr}a&b\\c&d\end{array}\right]\]

B = [x y; z w]

\[\left[ \begin{array}{rr}x&y\\z&w\end{array}\right]\]

Las operaciones matriciales básicas (suma, resta, multiplicación y división por la izquierda) operan de la misma forma que en matrices numéricas:

suma = A + B

\[\left[ \begin{array}{rr}a + x&b + y\\c + z&d + w\end{array}\right]\]

resta = A - B

\[\left[ \begin{array}{rr}a - x&b - y\\c - z&d - w\end{array}\right]\]

multip = A*B

\[\left[ \begin{array}{rr}a x + b z&a y + b w\\c x + d z&c y + d w\end{array}\right]\]

div_izq = A\B  # equivalentemente (A^-1)*B

\[\left[ \begin{array}{rr}\frac{- b z + d x}{a d - b c}&\frac{- b w + d y}{a d - b c}\\\frac{a z - c x}{a d - b c}&\frac{a w - c y}{a d - b c}\end{array}\right]\]

La inversa de una matriz se obtienen ya sea usando la función inv(matriz) o usando matriz^-1:

inversa = A^-1

\[\left[ \begin{array}{rr}\frac{d}{a d - b c}&- \frac{b}{a d - b c}\\- \frac{c}{a d - b c}&\frac{a}{a d - b c}\end{array}\right]\]

Para calcular la transpuesta de una matriz usamos la sintaxis matriz.T.

transp = A.T

\[\left[ \begin{array}{rr}a&c\\b&d\end{array}\right]\]

ATB = inv(A.T*A)

\[\left[ \begin{array}{rr}\frac{b^{2} + d^{2}}{a^{2} d^{2} - 2 a b c d + b^{2} c^{2}}&- \frac{a b + c d}{a^{2} d^{2} - 2 a b c d + b^{2} c^{2}}\\- \frac{a b + c d}{a^{2} d^{2} - 2 a b c d + b^{2} c^{2}}&\frac{a^{2} + c^{2}}{a^{2} d^{2} - 2 a b c d + b^{2} c^{2}}\end{array}\right]\]

Para calcular el determinante de una matriz usamos la sintaxis matriz.det().

A.det()

\[\begin{equation*}a d - b c\end{equation*}\]

Para generar las matrices especiales identidad, ceros y unos usamos las funciones eye(dimension), zeros(nfilas,ncols) y ones(nfilas,ncols), respectivamente. Estas tres funciones no han sido exportada al entorno global en el paquete sympy, por lo que tenemos que indicarle a Julia que dicha función viene de ese paquete usando sympy.función().

sympy.eye(2)

\[\left[ \begin{array}{rr}1&0\\0&1\end{array}\right]\]

sympy.zeros(2,2)

\[\left[ \begin{array}{rr}0&0\\0&0\end{array}\right]\]

sympy.ones(2,2)

\[\left[ \begin{array}{rr}1&1\\1&1\end{array}\right]\]

Finalmente, podemos calcular la matriz Jacobiana de un sistema de ecuaciones usando los aprendido hasta aquí. Definamos el sistema como una arreglo de ecuaciones:

f1 = 2*x^2 + 3*y^2 + 5*z^2
f2 = 2*x + 3*x*y + x*z
f3 = 3*x*y + 2*x + z^3
sistema = [f1, f2, f3]

\[ \left[ \begin{array}{r}2 x^{2} + 3 y^{2} + 5 z^{2}\\3 x y + x z + 2 x\\3 x y + 2 x + z^{3}\end{array} \right] \]

Calculamos la matriz Jacobiana del sistema diferenciando el sistema respecto de sus variables:

J = sistema.diff([x,y,z])

\(\displaystyle \left[\begin{matrix}\left[\begin{matrix}4 x\\3 y + z + 2\\3 y + 2\end{matrix}\right]\\\left[\begin{matrix}6 y\\3 x\\3 x\end{matrix}\right]\\\left[\begin{matrix}10 z\\x\\3 z^{2}\end{matrix}\right]\end{matrix}\right]\)

Note que sympy retorna un arreglo de arreglos con los 3 gradientes. Esta no es la forma a la que estamos acostumbrados de representar la matriz Jacobiana. Podemos escribir una función que reciba el arreglo anterior y lo reorganice de tal forma que tengamos una matriz con el número de filas igual al número de ecuaciones y el número de columnas igual al número de variables. Esto es:

function jacobian(mat)
    nv, nn,ne, nn = mat.shape
    J = Array{Sym}(undef, ne, nv)
    
    for i in 1:ne, j in 1:nv
        J[i,j] = mat[j][1][i][1]
    end
    
    return J
end

jacobian (generic function with 1 method)

Ahora aplicamos la función:

jacobian(J)

\[\left[ \begin{array}{rrr}4 x&6 y&10 z\\3 y + z + 2&3 x&x\\3 y + 2&3 x&3 z^{2}\end{array}\right]\]