knitr::opts_chunk$set(echo = TRUE)
library(Ryacas)
Ryacas
makes the yacas
computer algebra system available from
within R
. The name yacas
is short for "Yet Another Computer Algebra
System". The yacas
program
is developed by Ayal Pinkhuis and others,
and is available at http://www.yacas.org/ for various
platforms. There is a comprehensive documentation (300+ pages) of
yacas
available at https://yacas.readthedocs.io and the
documentation contains many examples.
This version of Ryacas
is somewhat different to previous versions of Ryacas
because we have tried to make the interface a lot simpler.
The old version of Ryacas
is available as a legacy version called Ryacas0
at https://github.com/r-cas/ryacas0/ with documentation directly available at https://r-cas.github.io/ryacas0/.
yacas
The naming principle governing Ryacas
functions is as follows:
yac_*(x)
functions evaluate/run yacas
command x
; the result varies depending on which of the functions used (see below)y_*(x)
various utility functions (not involving calls to yacas
)There are two interfaces to yacas
: a low-level (see the "The low-level interface" vignette) and a high-level (see the "The high-level (symbol) interface" vignette).
The low-level is highly customisable, but also requires more work with text strings.
The high-level is easier dealing with vectors and matrices, but may also be less computationally efficient and less flexible.
Below, we will demonstrate both interfaces and refer to the other vignettes for more information.
A short summary of often-used yacas
commands are found
at the end of this vignette.
A short summary of often-used low-level Ryacas
functions are found
at the end of the "The low-level interface" vignette,
and a short summary of often-used high-level Ryacas
functions are found
at the end of the "The high-level (symbol) interface" vignette.
The low-level interface consists of these two main functions:
yac_str(x)
: Evaluate yacas
command x
(a string) and get result as string/character.yac_expr(x)
: Evaluate yacas
command x
(a string) and get result as an R
expression.Note, that the yacas
command x
is a string and must often be built op using paste()
/paste0()
.
Examples of this will be shown in multiple examples below.
eq <- "x^2 + 4 + 2*x + 2*x" yac_str(eq) # No task was given to yacas, so we simply get the same returned yac_str(paste0("Simplify(", eq, ")")) yac_str(paste0("Factor(", eq, ")")) yac_expr(paste0("Factor(", eq, ")")) yac_str(paste0("TeXForm(Factor(", eq, "))"))
Instead of the pattern paste0("Simplify(", eq, ")")
etc., there exists a
helper function y_fn()
that does this:
y_fn(eq, "Factor") yac_str(y_fn(eq, "Factor")) yac_str(y_fn(y_fn(eq, "Factor"), "TeXForm"))
As you see, there are a lot of nested function calls. That can be avoided by using magrittr
's pipe %>%
(automatically available with Ryacas
) together with the helper function y_fn()
:
eq %>% y_fn("Factor") eq %>% y_fn("Factor") %>% yac_str() eq %>% y_fn("Factor") %>% y_fn("TeXForm") %>% yac_str()
The polynomial can be evaluated for a value of $x$ by calling yac_expr()
instead of yac_str()
:
yac_str(paste0("Factor(", eq, ")")) expr <- yac_expr(paste0("Factor(", eq, ")")) expr eval(expr, list(x = 2))
The high-level interface consists of the main function ysym()
and
often the helper function as_r()
will be used to get back an R
object (expression, matrix, vector, ...).
Before we had eq
as a text string. We now make a ysym()
from that:
eqy <- ysym(eq) eqy as_r(eqy) eqy %>% y_fn("Factor") # Notice how we do not need to call yac_str()/yac_expr()
Notice how the printing is different from before.
We start with a small matrix example:
A <- outer(0:3, 1:4, "-") + diag(2:5) a <- 1:4 B <- ysym(A) B b <- ysym(a) b
Notice how they are printed using yacas
's syntax.
We can apply yacas
functions using y_fn()
:
y_fn(B, "Transpose") y_fn(B, "Inverse") y_fn(B, "Trace")
Some standard R
commands are available (see the section "Ryacas
high-level reference" at the end of the "The high-level (symbol) interface" vignette):
A %*% a B %*% b t(A) t(B) A[, 2:3] B[, 2:3] A %*% solve(A) B %*% solve(B)
Next we will demonstrate matrix functionality using the Hilbert matrix
[
H_{{ij}}={\frac{1}{i+j-1}}
]
In R
's solve()
help file there is code for generating it:
hilbert <- function(n) { i <- 1:n H <- 1 / outer(i - 1, i, "+") return(H) }
To avoid floating-point issues (see the "Arbitrary-precision arithmetic" vignette),
we instead generate just the denominators as a stanard R
matrix:
hilbert_den <- function(n) { i <- 1:n H <- outer(i - 1, i, "+") return(H) } Hden <- hilbert_den(4) Hden H <- 1/Hden H
To use Ryacas
's high-level interface, we use the function ysym()
that converts the matrix to yacas
representation and automatically calls
yac_str()
when needed.
Furthermore, it enables standard R
functions such as subsetting with [
, diag()
, dim()
and others.
Hyden <- ysym(Hden) Hyden Hy <- 1/Hyden Hy
Notice how the printing is different from R
's printing.
We can then to a number of things with the ysym()
.
as_r(Hy) # now floating-point and the related problems diag(Hy) Hy[upper.tri(Hy)] Hy[1:2, ] dim(Hy) A <- Hy A[lower.tri(A)] <- "x" A as_r(A) eval(as_r(A), list(x = 999))
We consider the Rosenbrock function:
x <- ysym("x") y <- ysym("y") f <- (1 - x)^2 + 100*(y - x^2)^2 f tex(f)
$$\begin{align}f(x, y) = r tex(f)
\end{align}$$
We can visualise this, too.
N <- 30 x <- seq(-1, 2, length=N) y <- seq(-1, 2, length=N) f_r <- as_r(f) f_r z <- outer(x, y, function(x, y) eval(f_r, list(x = x, y = y))) levels <- c(0.001, .1, .3, 1:5, 10, 20, 30, 40, 50, 60, 80, 100, 500, 1000) cols <- rainbow(length(levels)) contour(x, y, z, levels = levels, col = cols)
Say we want to find the minimum. We do that by finding the roots of the gradient:
g <- deriv(f, c("x", "y")) g
$$\begin{align}g(x, y) = r tex(g)
\end{align}$$
crit_sol_all <- solve(g, c("x", "y")) crit_sol_all crit_sol <- crit_sol_all[1, ] %>% y_rmvars() crit_sol crit <- crit_sol %>% as_r() crit
We now verify what type of critical point we have by inspecting the Hessian at that critical point:
H <- Hessian(f, c("x", "y")) H tex(H)
$$\begin{align}
H = r tex(H)
\end{align}$$
H_crit <- eval(as_r(H), list(x = crit[1], y = crit[2])) H_crit eigen(H_crit, only.values = TRUE)$values
Because the Hessian is positive definite, the critical point $(1, 1)$ is indeed a minimum (actually, the global minimum).
yacas
referenceBelow are some yacas
functions. A more elaborate reference is available at https://yacas.readthedocs.io/:
Expand(x)
: Expand an expressionFactor(x)
: Factorise an expressionSimplify(x)
: Simplify an expressionSolve(expr, var)
solve an equation (refer to the Ryacas
function y_rmvars()
)Variables()
: List yacas
variablesD(x) expr
: Take the derivative of expr
with respect to x
HessianMatrix(function, var)
: Create the Hessian matrixJacobianMatrix(function, var)
: Create the Jacobian matrixLimit(n, a) f(n)
: Limit of f(n)
for n
going towards a
(e.g. Infinity
or 0
)Sum(k, a, b, f(k))
: Sum of f(k)
for k
from a
to b
.TeXForm(x)
: Get a $\LaTeX$ representation of an expressionPrettyForm(x)
: Print a prettier ASCII representation of an expressionInverse(A)
: Inverse of a matrixTranspose(A)
: Transpose of a matrixA * B
: Matrix multiplication (and not as R
's %*%
)Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.