Symbolic computation in R with Ryacas - results become character

I have a small MATLAB script mainly doing derivatives using symbolic toolbox that I want to rewrite into R. I chose Ryacas package because I found rSymPy too tricky to install... Here is my R code

# install.packages('Ryacas')
library(Ryacas)    
z <- Sym("z")
psi=c()
psi[1]=z^2*exp(-z)/(1-exp(-z))
psi[2]=z^2*exp(-z)/(1-exp(-z))*log(z)
psi[3]=z^2*exp(-z)/(1-exp(-z))*log(z)^2
f=matrix(NA,4,4)
f[1,1]=z^2*exp(-z)/(1-exp(-z))
for(i in 2:4){
  f[i,1]=deriv.Sym(psi[i-1],z)
  j=2
  while(j<=i){
    f[i,j]=deriv.Sym(expression(f[i,j-1]/f[j-1,j-1]),z)
    j=j+1
  }
}

It does not report any error. However, the output shows that R isn't actually doing symbolic computation but returns characters. So I cannot evaluate the result. I tried

> i=2
> deriv.Sym(psi[i-1],z)
expression(((1 - exp(-z)) * (2 * (z * exp(-z)) - z^2 * exp(-z)) - 
    z^2 * exp(-z)^2)/(1 - exp(-z))^2)
> f[i,1]
[1] "( D( z , 1 ) ( ( ( z ^ 2 ) * ( Exp ( ( - z ) ) ) ) / ( 1 - ( Exp ( ( - z ) ) ) ) ) )"

It seems that deriv.Sym(psi[i-1],z) does the symbolic derivative and get the correct result. But if the result is assigned to a variable, it becomes character class. I feel confused about expression(), yacas(), Sym() and character. Anyone can point out my mistake or help me clarify these concept? Thank you so much.

Below corresponding MATLAB code for reference. The MATLAB code works just fine.

syms c;
psi(1)=c^2*exp(-c)/(1-exp(-c));
psi(2)=c^2*exp(-c)/(1-exp(-c))*log(c);
psi(3)=c^2*exp(-c)/(1-exp(-c))*log(c)^2;

f(1,1)=c^2*exp(-c)/(1-exp(-c));
for i=2:4
    f(i,1)=diff(psi(i-1),c);
    j=2;
    while j<=i
        f(i,j)=diff(f(i,j-1)/f(j-1,j-1),c);
        j=j+1;
    end
end

g11=matlabFunction(f(1,1));
fplot(g11,[0,10])
figure
g22=matlabFunction(f(2,2));
fplot(g22,[0,10])
figure
g33=matlabFunction(f(3,3));
fplot(g33,[0,10])
figure
g44=matlabFunction(f(4,4));
fplot(g44,[0,10])

Solution

There are several problems with the R code in the question:

it is attempting to assign an S3 object to elements of a logical matrix:
```
typeof(NA)
## [1] "logical"
```
so R has converted it to character (since Sym objects are internally character) which is as far as it can go. f needs to be defined as a list with 2 dimensions so that it can hold such objects:
```
f <- matrix(list(), 4, 4)
```
since f is a list with 2 dimensions all references to elements of f should use double square brackets as in:
```
f[[1, 1]] <- z^2 * exp(-z) / (1 - exp(-z))
```
similarly psi should be initialized as:
```
psi <- list()
```
and then referenced as:
```
psi[[1]] <- z^2 * exp(-z) / (1 - exp(-z))
```
to evaluate f[[i, 1]] use Eval:
```
Eval(f[[i, 1]], list(z = 1))
## [1] 0.2432798
```
This also works but overwrites the Sym object z:
```
z <- 1
Eval(f[[i, 1]])
```
in general code should be calling the generic deriv and not by directly going to the specific method deriv.Sym

The revised code is in the section at the end which makes these changes as well as some stylistic improvements.

Suggest you review the vignette that comes with Ryacas. From the R console enter:

vignette("Ryacas")

Also review the Ryacas demos:

demo(package = "Ryacas")

Revised code

# install.packages('Ryacas')
library(Ryacas)    

z <- Sym("z")

psi <- list()
psi[[1]] <- z^2 * exp(-z) / (1 - exp(-z))
psi[[2]] <- z^2 * exp(-z) / (1 - exp(-z)) * log(z)
psi[[3]] <- z^2 * exp(-z) / (1 - exp(-z)) * log(z)^2

f <- matrix(list(), 4, 4)
f[[1,1]] <- z^2 * exp(-z) / (1 - exp(-z))
for(i in 2:4) {
  f[[i, 1]] <- deriv(psi[[i-1]], z)
  j <- 2
  while(j <= i) {
    f[[i, j]] <- deriv(f[[i, j-1]] / f[[j-1, j-1]], z)
    j <- j + 1
  }
}

i <- 2
deriv(psi[[i-1]], z)
f[[i, 1]]

Eval(f[[i, 1]], list(z = 1))