Computación Simbólica con R

Computación Simbólica con R

Ing. Luis Manfredo Reyes

Traducción y adaptación del documento:

http://www.di.fc.ul.pt/~jpn/r/symbolic/

R es actualmente el lenguaje de análisis   más utilizado en el mundo académico y de investigación, primero por su potencia de procesamiento, y en segundo lugar por ser un software libre.

R  puede ser descargado del siguiente sitio: http://www.r-project.org/

Existen versiones para Windows, Mac, Linux y Android

UN ERROR MUY COMUN POR TODOS LADOS ES: HAY GENTE QUE CREE QUE R ES UN PAQUETE ESTADÍSTICO (COMO SAS O SPSS).

Por ser un lenguaje y no un paquete, R tiene capacidades para muchos procesos, pero se debe reconocer que quizá donde más se le ha utilizado es en el campo de la estadística.

Una de las ventajas de R, es la enorme cantidad de librerías (o paquetes), que han sido desarrolladas para ampliar sus capacidades.

R no fue diseñado como un software matemático, por lo que no posee muchas capacidades de cálculo simbólico (Por ejemplo (x-1)⁷ ). Sin embargo, para los que desean realizar éste tipo de cálculos sin abandonar el ambiente de R, hay esperanzas.

Existen varios paquetes desarrollados para “salir” momentáneamente de R a realizar cálculos en programas especializados y “retornar” con el resultado

Para éste fin, los paquetes más utilizados son: rSymPy y Yacas.

rSympy permite realizar cálculos simbólicos usando el software Python

Yacas permite realizar cálculos simbólicos usando el software Yacas

En éste documento se discute el uso de rSymPy, debido a que posee más funciones disponibles.

REQUISITOS PREVIOS:

Debe estar instalado el software Python, su agregado Symphy y el Java

Referencias:

• Manual de Referencia (muy corto)
•  http://cran.r-project.org/web/packages/rSymPy/rSymPy.pdf
• Tutorial de SymPy: http://docs.sympy.org/latest/tutorial/

Cargar la Librería:

library(rSymPy) 

Definir x como variable simbólica

x <- Var("x")

Suma de dos veces x

x+x

## [1] "2*x"

Simplificación sencilla

x*x/x

## [1] "x"

Potencias

y <- Var("x**3")

x/y

## [1] "x**(-2)"

Expresiones

z <- sympy("2.5*x**2")

z + y

## [1] "2.5*x**2 + x**3"

Cálculos numéricos

Raíz cuadrada de 8

sympy("sqrt(8).evalf()") 

## [1] "2.82842712474619"

Raíz cuadrada de 8 con 50 decimales

sympy("sqrt(8).evalf(50)")

## [1] "2.8284271247461900976033774484193961571393437507539"

Pi con 120 decimales

sympy("pi.evalf(120)")

## [1] "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230665"

Trigonométricas

sympy("one = cos(1)**2 + sin(1)**2")

## [1] "cos(1)**2 + sin(1)**2"

Errores de redondeo

sympy("(one - 1).evalf()")

## [1] "-.0e-124"

Redondeo

sympy("(one - 1).evalf(chop=True)") 

## [1] "0"

Crear una ecuación

sympy("Eq(x**2+2*x+1,(x+1)**2)")

## [1] "1 + 2*x + x**2 == (1 + x)**2"

Asignar expression a variable

sympy("a = x**2+2*x+1")

## [1] "1 + 2*x + x**2"

sympy("b = (x+1)**2")

## [1] "(1 + x)**2"

Simplificar

"0" == sympy("simplify(a-b)") 

## [1] TRUE

Otras simplificaciones

sympy("simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))")

## [1] "-1 + x"

sympy("(x + 2)*(x - 3)")

## [1] "-(2 + x)*(3 - x)"

Desarrollo de productos

sympy("expand((x + 2)*(x - 3))")

## [1] "-6 - x + x**2"

Factorización

sympy("factor(x**3 - x**2 + x - 1)")

## [1] "-(1 + x**2)*(1 - x)"

Simplificaciones

y <- Var("y")

z <- Var("z")

sympy("collect(x*y + x - 3 + 2*x**2 - z*x**2 + x**3, x)")

## [1] "-3 + x*(1 + y) + x**2*(2 - z) + x**3"

sympy("(x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)")

## [1] "-(-2*y*z - 2*x*y*z + y**2 + z**2 + x*y**2 + x*z**2)/(1 - x**2)"

sympy("cancel((x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1))")

## [1] "(2*y*z - y**2 - z**2)/(1 - x)"

sympy("expand_func(gamma(x + 3))")

## [1] "x*(1 + x)*(2 + x)*gamma(x)"

Crear una variable “persistente”

sympy("y = x*x")

## [1] "x**2"

Matrices

sympy("A = Matrix([[1,x], [y,1]])")

## [1] "[   1, x]\n[x**2, 1]"

sympy("A**2")

## [1] "[1 + x**3,      2*x]\n[  2*x**2, 1 + x**3]"

sympy("B = A.subs(x,1.1)") 

## [1] "[   1, 1.1]\n[1.21,   1]"

sympy("B**2")

## [1] "[2.331,   2.2]\n[ 2.42, 2.331]"

sympy("expr = sin(2*x) + cos(2*x)")

## [1] "cos(2*x) + sin(2*x)"

sympy("expr.subs(sin(2*x), 2*sin(x)*cos(x))")

## [1] "2*cos(x)*sin(x) + cos(2*x)"

sympy("expr.subs(x,pi/2)")

## [1] "-1"

a1 <- Var("a1")

a2 <- Var("a2")

a3 <- Var("a3")

a4 <- Var("a4")

Inversa y Determinante de matrices

A <- Matrix(List(a1, a2), List(a3, a4))

Inv <- function(x) Sym("(", x, ").inv()")

Det <- function(x) Sym("(", x, ").det()")

A

## [1] "[a1, a2]\n[a3, a4]"

cat(A,"\n")

## ( Matrix( ( [ a1,a2 ] ),( [ a3,a4 ] ) ) )

Inv(A)

## [1] "[1/a1 + a2*a3/(a1**2*(a4 - a2*a3/a1)), -a2/(a1*(a4 - a2*a3/a1))]\n[            -a3/(a1*(a4 - a2*a3/a1)),        1/(a4 - a2*a3/a1)]"

Det(A)

## [1] "a1*a4 - a2*a3"

Crear función exponencial

Exp <- function(x) Sym("exp(", x, ")")

Exp(-x) * Exp(x)

## [1] "1"

y <- Var("y")

sympy("sqrt(8)")           

## [1] "2*2**(1/2)"

Resolver ecuación

sympy("solve(x**2 - 2, x)") 

## [1] "[2**(1/2), -2**(1/2)]"

Límites

sympy("limit(1/x, x, oo)")  

## [1] "0"

sympy("limit(1/x, x, 0)") 

## [1] "oo"

Integración

sympy("integrate(exp(-x))")             

## [1] "-exp(-x)"

sympy("integrate(exp(-x*y),x)")         

## [1] "-exp(-x*y)/y"

Integral definida

sympy("integrate(exp(-x), (x, 0, oo))")

## [1] "1"

Integración definida con función nativa de R

integrate( function(x) exp(-x), 0, Inf) 

## 1 with absolute error < 5.7e-05

sympy("integrate(x**2 - y, (x, -5, 5), (y, -pi, pi))")

## [1] "500*pi/3"

Derivada

sympy("diff(sin(2*x), x, 1)")

## [1] "2*cos(2*x)"

Función nativa de R para derivada

D( expression(sin(2*x)), "x" ) # also possible in base R

## cos(2 * x) * 2

Segunda derivada

sympy("diff(sin(2*x), x, 2)")

## [1] "-4*sin(2*x)"

Tercera derivada

sympy("diff(sin(2*x), x, 3)")

## [1] "-8*cos(2*x)"

Derivadas parciales

sympy("diff(exp(x*y*z), x, y, y)")

## [1] "z*exp(x*y*z) + x*y*z**2*exp(x*y*z)"

sympy("diff(exp(x*y*z), x, z, 3)")

## [1] "y*exp(x*y*z) + x*z*y**2*exp(x*y*z)"

Serie de Taylor

sympy("(1/cos(x)).series(x, 0, 10)")

## [1] "1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + O(x**10)"

sympy("exp(x).series(x, 0, 5)")  

## [1] "1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)"

sympy("exp(x).series(x, 0, 5).removeO()")

## [1] "1 + x + x**2/2 + x**3/6 + x**4/24"

Calcular eigenvalores de una matriz

sympy("Matrix([[1, 2], [2, 2]]).eigenvals()") 

## [1] "{3/2 - 17**(1/2)/2: 1, 3/2 + 17**(1/2)/2: 1}"