使用 Ryacas 在 R 中进行符号计算 - 结果变成字符

Symbolic computation in R with Ryacas - results become character

我有一个小的 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 实际上不是在进行符号计算,而是 returns 个字符。所以我无法评估结果。我试过了

> 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)做了符号求导得到了正确的结果。但是如果把结果赋值给变量,就变成了字符class。我对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])

题中的R代码有几个问题:

  • 它正在尝试将 S3 对象分配给逻辑矩阵的元素:

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

    因此 R 已将其转换为字符(因为 Sym 对象在内部 字符)这是它可以去的。 f需要定义为列表 具有 2 个维度,以便它可以容纳这样的对象:

    f <- matrix(list(), 4, 4)
    
  • 因为 f 是一个二维列表,所有对 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))