是否可以在 R 中求解代数方程?

Is it possible to solve an algebraic equation in R?

我想找到解决方案:

-x^3+6*x^2+51*x+44=0

但是有了 R。这可能吗?

我找到了 Ryacas 包,但似乎没有人能够让它工作。

听起来可能微不足道,但我无法找到一种简单的方法...

你有其他选择吗?

谢谢大家!

您可以使用 polynom 包:

library(polynom)
p <- polynomial(c(44,51,6,-1))
# 44 + 51*x + 6*x^2 - x^3 
solve(p)
# [1] -4 -1 11

但是您可以简单地使用 base 包中的函数 polyroot

polyroot(c(44,51,6,-1))
# [1] -1+0i -4+0i 11+0i

如果保留实部 Re:

Re(polyroot(c(44,51,6,-1)))
# [1] -1 -4 11

这里我们使用矩阵与其特征多项式之间的关系求解根。

给定多项式a0 + a1*x^1 + a2*x^2 + x^3,定义矩阵:

0   0  -a0
1   0  -a1
0   1  -a2

这个矩阵的特征值是多项式的根。

在多项式方程中代入 y = -x 得到

y^3 + 6*y^2 - 51*y + 44=0

并举出这个例子

> z <- matrix(c(0,1,0,0,0,1,-44,51,-6),3,3)
> z
     [,1] [,2] [,3]
[1,]    0    0  -44
[2,]    1    0   51
[3,]    0    1   -6
> eigen(z)
$values
[1] -11   4   1

$vectors
           [,1]        [,2]        [,3]
[1,]  0.6172134  0.73827166  0.98733164
[2,] -0.7715167 -0.67115606 -0.15707549
[3,]  0.1543033 -0.06711561 -0.02243936

或者,因为我们用 -y 代替了 x

> eigen(-z)$values
[1] 11 -4 -1

参见:http://www-math.mit.edu/~edelman/publications/polynomial_roots.pdf

我只是偶然发现了这个问题,我不确定 Ryacas 包是否有任何内在的变化,但它似乎在 2020 年运行良好,这里有一个有用的小插图:https://cran.r-project.org/web/packages/Ryacas/vignettes/getting-started.html

按照小插图,当我 运行 代码时,事情似乎按预期工作:

library(Ryacas)
# initialize equation:
eq <- "-x^3+6*x^2+51*x+44"
# simplify the equation:
library(glue)
yac_str(glue("Simplify({eq})"))

[1] "6*x^2-x^3+51*x+44"

# factor:
yac_str(glue("Factor({eq})"))

[1] "(-1)*(x-11)*(x+4)*(x+1)"

您可以像这样计算表达式,插入 x 的任何值:

# evaluate
evaluate(eq,list(x=c(0,1,10,100,-100)))
[[1]]
$src
[1] "-x^3+6*x^2+51*x+44"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

在这里您可以看到结果,其中 x=0 产生了 44 的答案,x=1 产生了 100 的答案,等等...

如果您评估新的简化版本或分解版本并评估它们,您当然会得到完全相同的结果:

evaluate(yac_str(glue("Simplify({eq})")),list(x=c(0,1,10,100,-100)))

[[1]]
$src
[1] "6*x^2-x^3+51*x+44"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

请注意 $src 输出中的公式已更改,但我们得到相同的结果。

这也是因式分解: evaluate(yac_str(glue("Factor({eq})")),list(x=c(0,1,10,100,-100)))

[[1]]
$src
[1] "(-1)*(x-11)*(x+4)*(x+1)"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

我在此处概述的内容与小插图中概述的内容之间唯一真正的区别是实际公式,以及我使用 library(glue) 而不是 paste0 的事实,这也是一个合理的选择.