R 中的显式公式与符号导数

Question

我想评估R中某些函数f的高阶导数。 我有两种可能性。

1. 要么我确定 f^(k) 的通用表达式，即 f 的 k 阶导数（我在我的具体情况下可以做），然后我评估它；
2. 或者我利用R的能力来执行符号导数（函数D()）。

1比2有什么优点？假设 f^(k) 不是递归公式。如果 f^(k) 是递归的怎么办？

任何提示将不胜感激。

Answer 1

符号微分比手工微分更不易出错。

对于低阶，我不认为符号微分会花费太多计算机时间，但您可以轻松地计算您的具体情况，以确定它正在使用proc.time, system.time 或 rbenchmark 包。另请参阅这些示例。

您可能想尝试符号微分和手微分作为检查。

关于 R 的 deriv （以及相关函数例如 D）与 Ryacas 包相比，后者能够重复执行差异化而不要求用户迭代它们自己（deriv 的第三个参数指定顺序），它有一个 Simplify 函数，R 中不存在对应的函数，尽管 Ryacas 应该仔细检查，因为 yacas 可能有点有时有越野车。

这是一个示例：

> library(Ryacas)
> x <- Sym("x")
> y <- (x^2+x)^2
> dy <- deriv(y, x, 2) # 2nd deriv
> dy
expression(2 * (2 * x + 1)^2 + 4 * (x^2 + x))
> Simplify(dy)
expression(2 * (6 * x^2 + 6 * x + 1))

编辑：添加到示例：

> Eval(dy, list(x = 3))
[1] 146
> Eval(Simplify(dy), list(x = 3))
[1] 146
>
> f <- function(x) {}
> body(f) <- yacas(Simplify(dy))[[1]]
> f
function (x) 
2 * (6 * x^2 + 6 * x + 1)
> f(3)
[1] 146
>
> # double check
> w <- 3
> eval(D( D(expression((w^2+w)^2), "w"), "w"))
[1] 146

另请尝试 demo("Ryacas-Function") 。

编辑 2：修复了一个错误并向示例添加了更多内容。