# Doing Maths Symbolically: R as a Computer Algebra System (CAS) When I first saw the Computer Algebra System Mathematica in the nineties I was instantly fascinated by it: you could not just calculate things with it but solve equations, simplify, differentiate and integrate expressions and even solve simple differential equations… not just numerically but symbolically! It helped me a lot during my studies at the university back then. Normally you cannot do this kind of stuff with R but fear not, there is, of course, a package for that!

There are many so-called Computer Algebra Systems (CAS) out there, commercial but also open-source. One very mature one is called YACAS (for Yet Another Computer Algebra System). You find the documentation here: Yacas Documentation (many of the following examples are taken from there).

You can use the full power of it in R by installing the Ryacas package from CRAN. You can use Yacas commands directly by invoking the yac_str function, the as_r function converts the output to R. Let us first simplify a mathematical expression:

library(Ryacas)
##
## Attaching package: 'Ryacas'
## The following object is masked from 'package:stats':
##
##     deriv
## The following objects are masked from 'package:base':
##
##     %*%, determinant, diag, diag<-, I, lower.tri, upper.tri

# simplify expressions
as_r(yac_str("Simplify(a*b*a^2/b-a^3)"))
##  0


Or solve an equation:

as_r(yac_str("Solve(a+x*y==z,x)"))
##  "x==-(a-z)/y"


And you can do all kinds of tedious stuff that is quite error-prone when done differently, e.g. expanding expressions like by using the binomial theorem:

as_r(yac_str("Expand((x-2)^20)"))
## expression(x^20 - 40 * x^19 + 760 * x^18 - 9120 * x^17 + 77520 *
##     x^16 - 496128 * x^15 + 2480640 * x^14 - 9922560 * x^13 +
##     32248320 * x^12 - 85995520 * x^11 + 189190144 * x^10 - 343982080 *
##     x^9 + 515973120 * x^8 - 635043840 * x^7 + 635043840 * x^6 -
##     508035072 * x^5 + 317521920 * x^4 - 149422080 * x^3 + 49807360 *
##     x^2 - 10485760 * x + 1048576)


To demonstrate how easily the results can be integrated into R let us do some curve sketching on a function. First, we define two helper function for converting an expression into a function (which can then be used to plot it) and for determining the derivative of order n of some function (we redefine the D function for that):

as_function <- function(expr) {
as.function(alist(x =, eval(parse(text = expr))))
}

# redefine D function
D <- function(eq, order = 1) {
yac_str(paste("D(x,", order, ")", eq))
}


Now, we define the function (in this case a simple polynomial ), determine the first and second derivatives symbolically and plot everything:

xmin <- -5
xmax <- 5

eq <- "2*x^3 - 3*x^2 + 4*x - 5"
eq_f <- as_function(eq)
curve(eq_f, xmin, xmax, ylab = "y(x)")
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)

D_eq <- D(eq)
D_eq
##  "6*x^2-6*x+4"

D_eq_f <- as_function(D_eq)
curve(D_eq_f, xmin, xmax, add = TRUE, col = "blue")

D2_eq <- D(eq, 2)
D2_eq
##  "12*x-6"

D2_eq_f <- as_function(D2_eq)
curve(D2_eq_f, xmin, xmax, add = TRUE, col = "green") Impressive, isn’t it! Yacas can also determine limits and integrals:

# determine limits
yac_str("Limit(x,0) 1/x")
##  "Undefined"

yac_str("Limit(x,0,Left) 1/x")
##  "-Infinity"

yac_str("Limit(x,0,Right) 1/x")
##  "Infinity"

# integration
yac_str("Integrate(x) Cos(x)")
##  "Sin(x)"

yac_str("Integrate(x,a,b) Cos(x)")
##  "Sin(b)-Sin(a)"


As an example, we can prove in no-time that the famous approximation is actually too big (more details can be found here: Proof that 22/7 exceeds π): yac_str("Integrate(x,0,1) x^4*(1-x)^4/(1+x^2)")
##  "22/7-Pi"


And, as the grand finale of this post, Yacas is even able to solve ordinary differential equations symbolically! Let us first take the simplest of them all:

as_r(yac_str("OdeSolve(y' == y)"))
## expression(C115 * exp(x))


It correctly returns the exponential function (The C-term is just an arbitrary constant). And finally a more complex, higher-order one:

as_r(yac_str("OdeSolve(y'' - 4*y == 0)"))
## expression(C154 * exp(2 * x) + C158 * exp(-2 * x))


I still find CAS amazing and extremely useful… and an especially powerful one can be used from within R!

## 11 thoughts on “Doing Maths Symbolically: R as a Computer Algebra System (CAS)”

1. vlad says:

Hi,
Thank you!
Interesting!

1. Learning Machines says:

You are very welcome (as always)!

2. Carl G Witthoft says:

OMG you used eval(parse()) — the Voldemort of R programming 🙂 . In all seriousness, you’d be better off passing a closure (function) as the argument rather than sending in a string to be interpreted.

1. Learning Machines says:

Ok, why is that “evil”? And could you provide better code? I would appreciate that… Thank you, Carl!

(btw: sent you a LinkedIn invite)

1. adam Joseph ginensky says:

I’ve found eval(parse()) to be useful in many circumstances and I’ve always wondered why it is considered so horrible. There is a stackoverflow exchange which quotes some evidently well-known computer expert as saying, “if eval is the answer, you’re asking the wrong question”, however as far as I can tell, he is discussing things like eval(3+2). There are times when I think eval is the perfect tool. That’s my two cents.

1. Learning Machines says:

Thank you, Adam. I see it the same way!

3. Carl G Witthoft says:

fortune(106)
If the answer is parse() you should usually rethink the question.
— Thomas Lumley
R-help (February 2005)

fortune(181)
Personally I have never regretted trying not to underestimate my own future stupidity.
— Greg Snow (explaining why eval(parse(…)) is often suboptimal, answering a question triggered by
the infamous fortune(106))
R-help (January 2007)

Anyway, curve() is a bit of a pain, and I am trying here to allow you to stick with the objects of type “expression” that things like the “D” function work with. So,

eqx <- expression(2*x^3 - 3*x^2 + 4*x - 5)
do.call(curve,list(eqx,from=1,to=5))
do.call(curve,list(D(eqx,"x"),from=1,to=5))


Should suffice.

1. Learning Machines says:

Thank you, Carl.

I still don’t understand why eval(parse(…)) is suboptimal in this case…

Having said that I see that your code is elegantly approaching the matter. To include it in my example above I guess I would have to restructure as_function and my redefined D function.

2. Learning Machines says:

…ok, I found this on SE: What specifically are the dangers of eval(parse(…))?.

It seems to boil down to two main points: security (which is of no concern here because it is not a live application on the internet where people can inject malicious code via obscure formulas) and debuggability (although in this context the construct seems clear and easy enough).

4. adam Joseph ginensky says:

I’ve found eval(parse()) to be useful in many circumstances and I’ve always wondered why it is considered so horrible. There is a stackoverflow exchange which quotes some evidently well-known computer expert as saying, “if eval is the answer, you’re asking the wrong question”, however as far as I can tell, he is discussing things like eval(3+2). There are times when I think eval is the perfect tool. That’s my two cents.

This site uses Akismet to reduce spam. Learn how your comment data is processed.