Introduction

Start with a base symbol what can either be:

• A yacas command, e.g. x, 2*a or something similar
• An R matrix or vector.

Here, we keep it simple.

x <- ysym("x")
2*x^2 - 5

## y: 2*x^2-5

c(-2, 5)*x

## {(-2)*x, 5*x}

c(2*x, -x^3)

## {2*x, -x^3}

as_r(c(-2, 5)*x) # or yac_expr(c(-2, 5)*x)

## expression(c(-2 * x, 5 * x))

Then consider an R matrix and vector:

A <- outer(0:3, 1:4, "-") + diag(2:5)
a <- 1:4
A

##      [,1] [,2] [,3] [,4]
## [1,]    1   -2   -3   -4
## [2,]    0    2   -2   -3
## [3,]    1    0    3   -2
## [4,]    2    1    0    4

a

## [1] 1 2 3 4

They are now considered yacas-enabled:

B <- ysym(A)
B

## {{ 1, -2, -3, -4},
##  { 0,  2, -2, -3},
##  { 1,  0,  3, -2},
##  { 2,  1,  0,  4}}

as_r(B)

##      [,1] [,2] [,3] [,4]
## [1,]    1   -2   -3   -4
## [2,]    0    2   -2   -3
## [3,]    1    0    3   -2
## [4,]    2    1    0    4

b <- ysym(a)
b

## {1, 2, 3, 4}

as_r(b)

## [1] 1 2 3 4

Notice how they are printed using yacas’s syntax.

We can apply yacas functions using y_fn():

y_fn(B, "Transpose")

## {{ 1,  0,  1,  2},
##  {-2,  2,  0,  1},
##  {-3, -2,  3,  0},
##  {-4, -3, -2,  4}}

y_fn(B, "Inverse")

## {{   37/202,     3/101,    41/202,    31/101},
##  {(-17)/101,    30/101,     3/101,     7/101},
##  {(-19)/202,  (-7)/101,    39/202,  (-5)/101},
##  { (-5)/101,  (-9)/101, (-11)/101,     8/101}}

y_fn(B, "Trace")

## y: 10

Standard R commands are available (see the section “Ryacas high-level reference” at the end of this document):

A %*% a

##      [,1]
## [1,]  -28
## [2,]  -14
## [3,]    2
## [4,]   20

B %*% b

## {-28, -14, 2, 20}

t(A)

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    1    2
## [2,]   -2    2    0    1
## [3,]   -3   -2    3    0
## [4,]   -4   -3   -2    4

t(B)

## {{ 1,  0,  1,  2},
##  {-2,  2,  0,  1},
##  {-3, -2,  3,  0},
##  {-4, -3, -2,  4}}

exp(B)

## {{ Exp(1), Exp(-2), Exp(-3), Exp(-4)},
##  {      1,  Exp(2), Exp(-2), Exp(-3)},
##  { Exp(1),       1,  Exp(3), Exp(-2)},
##  { Exp(2),  Exp(1),       1,  Exp(4)}}

as_r(exp(B))

##          [,1]      [,2]        [,3]        [,4]
## [1,] 2.718282 0.1353353  0.04978707  0.01831564
## [2,] 1.000000 7.3890561  0.13533528  0.04978707
## [3,] 2.718282 1.0000000 20.08553692  0.13533528
## [4,] 7.389056 2.7182818  1.00000000 54.59815003

A[, 2:3]

##      [,1] [,2]
## [1,]   -2   -3
## [2,]    2   -2
## [3,]    0    3
## [4,]    1    0

B[, 2:3]

## {{-2, -3},
##  { 2, -2},
##  { 0,  3},
##  { 1,  0}}

A[upper.tri(A)] <- 1
B[upper.tri(B)] <- 1
A

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    0    2    1    1
## [3,]    1    0    3    1
## [4,]    2    1    0    4

B

## {{1, 1, 1, 1},
##  {0, 2, 1, 1},
##  {1, 0, 3, 1},
##  {2, 1, 0, 4}}

2*A - A

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    0    2    1    1
## [3,]    1    0    3    1
## [4,]    2    1    0    4

2*B - B

## {{1, 1, 1, 1},
##  {0, 2, 1, 1},
##  {1, 0, 3, 1},
##  {2, 1, 0, 4}}

A %*% solve(A)

##               [,1]          [,2]         [,3]         [,4]
## [1,]  1.000000e+00 -1.110223e-16 5.551115e-17 5.551115e-17
## [2,] -1.110223e-16  1.000000e+00 5.551115e-17 5.551115e-17
## [3,]  2.220446e-16 -1.110223e-16 1.000000e+00 0.000000e+00
## [4,]  0.000000e+00 -2.220446e-16 0.000000e+00 1.000000e+00

B %*% solve(B)

## {{1, 0, 0, 0},
##  {0, 1, 0, 0},
##  {0, 0, 1, 0},
##  {0, 0, 0, 1}}

solve(A %*% t(A))

##            [,1]       [,2]       [,3]       [,4]
## [1,]  3.7040816 -1.2551020 -0.8877551 -0.6224490
## [2,] -1.2551020  0.6938776  0.2346939  0.1530612
## [3,] -0.8877551  0.2346939  0.3367347  0.1326531
## [4,] -0.6224490  0.1530612  0.1326531  0.1734694

solve(B %*% t(B))

## {{   363/98, (-123)/98,  (-87)/98,  (-61)/98},
##  {(-123)/98,     34/49,     23/98,     15/98},
##  { (-87)/98,     23/98,     33/98,     13/98},
##  { (-61)/98,     15/98,     13/98,     17/98}}

solve(A, a)

## [1] -1.2857143 -0.2857143  0.8571429  1.7142857

solve(B, b)

## {(-9)/7, (-2)/7, 6/7, 12/7}

We can also assign a yacas variable, but remember that this may be difficult to distinguish:

yac_str("W") # Get variable W if exists, or else just a symbol

## [1] "W"

yac_str("Variables()") # ...or list variables

## [1] "{t,rformBitwiseOps,TeXForm'FuncPrec,I,j,Arow}"

B

## {{1, 1, 1, 1},
##  {0, 2, 1, 1},
##  {1, 0, 3, 1},
##  {2, 1, 0, 4}}

yac_assign(B, "W") # assign B in R to W in yacas
yac_str("W") # Get variable W if exists, or else just a symbol

## [1] "{{1,1,1,1},{0,2,1,1},{1,0,3,1},{2,1,0,4}}"

yac_str("Variables()") # ...or list variables

## [1] "{W,t,rformBitwiseOps,TeXForm'FuncPrec,I,j,Arow}"

yac_silent("Clear(W)")
yac_str("Variables()") # List variables

## [1] "{t,rformBitwiseOps,TeXForm'FuncPrec,I,j,Arow}"

yac_str("W") # Get variable W if exists, or else just a symbol

## [1] "W"

Simplify and output to TeX

There are additional functions available:

• simplify()
• tex()

To demonstrate these and some additional benefit, we exploit yacas’s symbolic availabilities.

D <- diag(4) %>% ysym()
D

## {{1, 0, 0, 0},
##  {0, 1, 0, 0},
##  {0, 0, 1, 0},
##  {0, 0, 0, 1}}

D <- D/2
D

## {{1/2,   0,   0,   0},
##  {  0, 1/2,   0,   0},
##  {  0,   0, 1/2,   0},
##  {  0,   0,   0, 1/2}}

D[2:3, 1] <- "d"
D[3, 4] <- "2*d + 2"
D

## {{  1/2,     0,     0,     0},
##  {    d,   1/2,     0,     0},
##  {    d,     0,   1/2, 2*d+2},
##  {    0,     0,     0,   1/2}}

D %>% solve()

## {{           2,            0,            0,            0},
##  {      (-4)*d,            2,            0,            0},
##  {      (-4)*d,            0,            2, (-4)*(2*d+2)},
##  {           0,            0,            0,            2}}

D %>% solve() %>% simplify()

## {{         2,          0,          0,          0},
##  {    (-4)*d,          2,          0,          0},
##  {    (-4)*d,          0,          2, (-8)*(d+1)},
##  {         0,          0,          0,          2}}

D %>% solve() %>% simplify() %>% tex()

## [1] "\\left( \\begin{array}{cccc} 2 & 0 & 0 & 0 \\\\ -4 d & 2 & 0 & 0 \\\\ -4 d & 0 & 2 & -8 \\left( d + 1\\right)  \\\\ 0 & 0 & 0 & 2 \\end{array} \\right)"

\[ \left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ -4 d & 2 & 0 & 0 \\ -4 d & 0 & 2 & -8 \left( d + 1\right) \\ 0 & 0 & 0 & 2 \end{array} \right) \]

yacas has a Simplify() function. This is made available via a simplify() function that also includes a time-out that prevents yacas in making the R session hang, but it requires that the unix package is available. The default timeout value used when unix is available is 2 seconds.

Derivatives

We illustrate using the example in https://mikl.dk/post/2019-pizza-frozen-yogurt-goebner/:

L <- ysym("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
L

## y: (x^2*y)/4-a*(3*x+(3*y)/2-45)

We can consider one variable only:

deriv(L, "x")

## y: (x*y)/2-3*a

Hessian(L, "x")

## {{y/2}}

Or multiple variables:

deriv(L, c("x", "y", "a"))

## {(x*y)/2-3*a, x^2/4-(3*a)/2, 45-(3*x+(3*y)/2)}

H <- Hessian(L, c("x", "y", "a"))
H

## {{   y/2,    x/2,     -3},
##  {   x/2,      0, (-3)/2},
##  {    -3, (-3)/2,      0}}

as_r(H)

## expression(rbind(c(y/2, x/2, -3), c(x/2, 0, -3/2), c(-3, -3/2, 
##     0)))

eval(as_r(H), list(x = 2, y = 2, a = 2))

##      [,1] [,2] [,3]
## [1,]    1  1.0 -3.0
## [2,]    1  0.0 -1.5
## [3,]   -3 -1.5  0.0

The Jacobian is taken on a vector function denoted by many functions:

L2 <- ysym(c("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)", 
                   "x^3 + 4*a^2")) # just some function
L2

## {(x^2*y)/4-a*(3*x+(3*y)/2-45), x^3+4*a^2}

Jacobian(L2, "x")

## {{(x*y)/2-3*a},
##  {      3*x^2}}

Jacobian(L2, c("x", "y", "a"))

## {{     (x*y)/2-3*a,    x^2/4-(3*a)/2, 45-(3*x+(3*y)/2)},
##  {           3*x^2,                0,              8*a}}

Solving equations

Say we want to find roots of a polynomial. We use the generic solve(a, b, ...) function.

Note the conventions are as follows:

• Linear system of equations (works as R’s solve() as demonstrated above):
  • When a is a matrix and b not provided, this finds the inverse of a.
  • When a is a matrix and a vector b is provided, the linear system of equations is solved.
• Else:
  • solve(a, b): find roots of a for variable b, i.e. yacas Solve(a == 0, b)
  • solve(a, b, v): find solutions to a == b for variable v, i.e. yacas Solve(a == b, v)

xs <- ysym("x")
poly <- xs^2 - xs - 6
poly

## y: x^2-x-6

zeroes <- solve(poly, "x") # Solve(x^2 - x - 6 == 0, x)
zeroes

## {x==(-2), x==3}

tex(zeroes)

## [1] "\\left( x = -2, x = 3\\right)"

zeroes %>% y_rmvars()

## {-2, 3}

We can also find values of x where the polynomial equals another constant. If we were working with strings via the low-level interface it would be easy via paste(), but as we are working with ysym()’s we use the solve() function directly:

solve(poly, 3, "x") # Solve(x^2 - x - 6 == 3, x)

## {x==(Sqrt(37)+1)/2, x==(1-Sqrt(37))/2}

solve(poly, 3, "x") %>% tex()

## [1] "\\left( x = \\frac{\\sqrt{37} + 1}{2} , x = \\frac{1 - \\sqrt{37}}{2} \\right)"

\[ \left( x = \frac{\sqrt{37} + 1}{2} , x = \frac{1 - \sqrt{37}}{2} \right) \]

x <- ysym("x")
y <- ysym("y")
lhs <- c(3*x*y - y, x)
rhs <- c(-5*x, y+4)

sol <- solve(lhs, rhs, c("x", "y"))
sol

## {{      x==2,    y==(-2)},
##  {    x==2/3, y==(-10)/3}}

sol_vals <- lapply(seq_len(nrow(sol)), function(sol_no) {
  y_rmvars(sol[sol_no, ])
})
sol_vals

## [[1]]
## {2, -2} 
## 
## [[2]]
## {2/3, (-10)/3}

sol_envir <- lapply(sol_vals, function(l) {
  list(x = as_r(l[1]), y = as_r(l[2]))
})
sol_envir

## [[1]]
## [[1]]$x
## [1] 2
## 
## [[1]]$y
## [1] -2
## 
## 
## [[2]]
## [[2]]$x
## [1] 0.6666667
## 
## [[2]]$y
## [1] -3.333333

do.call(rbind, lapply(seq_along(sol_envir), function(sol_no) {
  sol_val <- sol_envir[[sol_no]]
  data.frame(sol_no = sol_no,
             eq_no = seq_along(sol_val),
             lhs = eval(as_r(lhs), sol_val),
             rhs = eval(as_r(rhs), sol_val))
}))

##   sol_no eq_no         lhs         rhs
## 1      1     1 -10.0000000 -10.0000000
## 2      1     2   2.0000000   2.0000000
## 3      2     1  -3.3333333  -3.3333333
## 4      2     2   0.6666667   0.6666667

Ryacas high-level reference

Principle:

• ysym(x) converts x to a yac_symbol that automatically runs yacas when needed. x can both be a text string with yacas commands or an R vector/matrix.
• as_r(x): Is used to convert the yac_symbol back to an R representation.
• y_fn(x, fn, ...): Apply a yacas function fn to the yac_symbol x, i.e. fn(x, ...); note that this is evaluated immediately when x is a yac_symbol as opposed to when x is a string

Reference:

The following functions work with yac_symbols.

• ysym(): Create yac_symbol
  • ysym_make(): Make a vector of yac_symbols
  • ysym_ls(): List declared yac_symbols
• yac_*() functions (see the “Getting started” vignette)
  • yac_str(): Return yacas string
  • yac_expr(): Return R expression
  • yac_silent(): Do something silently
  • yac_assign(): Assign a variable
• Other
  • simplify(x, timeout = 2): Try yacas’s Simplify() function. When the unix package is available, the timeout (in seconds), stops trying after that amount of time to avoid making the R process hang.
  • tex(): Convert
  • y_fn(x, fn, ...): Apply a yacas function fn to the yac_symbol x, i.e. fn(x, ...)
  • y_rmvars(x): Remove variable names in x
  • Derivatives:
    • deriv(expr, vars): takes derivative of yac_symbol expr with respect to vars
    • Jacobian(expr, vars): finds Jacobian of yac_symbol expr (usually a vector of expressions) with respect to vars
    • Hessian(expr, vars): finds Hessian matrix of yac_symbol expr with respect to vars
  • Other
    • lim()
• Standard R that has been implemented for yac_symbols:
  • print()
  • c()
  • dim()
  • cbind()
  • rbind()
  • [ getter
  • [<- setter
  • [[ getter
  • %*% matrix/vector multiplication
  • diag() getter
  • diag<-() setter
  • upper.tri() getter
  • lower.tri() getter
  • t()
  • solve() (see above and in help page)
  • integrate()
  • sum()
  • prod()
  • Ops: +, -, *, /, ^
  • Math functions: sin(), cos(), tan(), asin(), acos(), atan(), asinh(), acosh(), atanh(), exp(), log(), sqrt()