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)")) ## [1] 0

Or *solve an equation*:

as_r(yac_str("Solve(a+x*y==z,x)")) ## [1] "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 ## [1] "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 ## [1] "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") ## [1] "Undefined" yac_str("Limit(x,0,Left) 1/x") ## [1] "-Infinity" yac_str("Limit(x,0,Right) 1/x") ## [1] "Infinity" # integration yac_str("Integrate(x) Cos(x)") ## [1] "Sin(x)" yac_str("Integrate(x,a,b) Cos(x)") ## [1] "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)") ## [1] "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!

Hi,

Thank you!

Interesting!

Vlad

You are very welcome (as always)!

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.Ok, why is that “evil”? And could you provide better code? I would appreciate that… Thank you, Carl!

(btw: sent you a LinkedIn invite)

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.Thank you, Adam. I see it the same way!

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,Should suffice.

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.…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).

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.