2

我有一个小的 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])
4

1 回答 1

3

问题中的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))
于 2018-01-01T13:36:31.517 回答