Символьные вычисления в R с Ryacas - результаты становятся символами

Question

У меня есть небольшой скрипт MATLAB, в основном выполняющий производные с использованием символического набора инструментов, который я хочу переписать в R. Я выбрал пакет Ryacas, потому что нашел rSymPy слишком сложным для установки... Вот мой код R

# 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
  }
}

Он не сообщает об ошибке. Однако выходные данные показывают, что R на самом деле не выполняет символьные вычисления, а возвращает символы. Поэтому я не могу оценить результат. Я старался

> 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 ) ) ) ) ) )"

Кажется, что deriv.Sym(psi[i-1],z) делает символическую производную и получаю правильный результат. Но если результат присваивается переменной, он становится классом символов. Я смущен expression(), yacas(), Sym() и характер. Кто-нибудь может указать на мою ошибку или помочь мне прояснить эту концепцию? Огромное спасибо.

Ниже приведен соответствующий код MATLAB для справки. Код MATLAB работает просто отлично.

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])

Answer 1

Есть несколько проблем с кодом R в вопросе:

• он пытается присвоить объект S3 элементам логической матрицы:

  typeof(NA)
  ## [1] "logical"
  

  поэтому R преобразовал его в символ (так как Sym объекты внутренне характерны), насколько это возможно. f должен быть определен как список с двумя измерениями, чтобы он мог содержать такие объекты:

  f <- matrix(list(), 4, 4)
  

• поскольку f это список с 2 измерениями все ссылки на элементы f следует использовать двойные квадратные скобки, как в:

  f[[1, 1]] <- z^2 * exp(-z) / (1 - exp(-z))
  

• так же psi должен быть инициализирован как:

  psi <- list()
  

  и затем упоминается как:

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

• оценить f[[i, 1]] использование Eval:

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

  Это также работает, но перезаписывает Sym объект z:

  z <- 1
  Eval(f[[i, 1]])
  

• в общем код должен вызывать общий deriv а не путем прямого перехода к конкретному методу deriv.Sym

Пересмотренный код находится в конце раздела, в котором вносятся эти изменения, а также некоторые стилистические улучшения.

Предлагаем вам ознакомиться с виньеткой, которая идет с Ryacas. С консоли R введите:

vignette("Ryacas")

Также ознакомьтесь с демонстрациями Ryacas:

demo(package = "Ryacas")

Пересмотренный код

# 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))