R中是否有一个函数/包,它接受一个函数 f 和一个参数 k,然后返回 k 次 f 的泰勒近似值?
utdiscant
问问题
10763 次
5 回答
16
您可以使用Ryacas
yacas计算机代数系统(您也需要安装)
使用小插图中的示例
library(Ryacas)
# run yacasInstall() if prompted to install yacas
#
yacas("texp := Taylor(x,0,3) Exp(x)")
## expression(x + x^2/2 + x^3/6 + 1)
# or
现在,如果你想把它变成一个你可以给出值的函数x
myTaylor <- function(f, k, var,...){
.call <- sprintf('texp := Taylor( %s, 0, %s) %s', var,k,f)
result <- yacas(.call)
foo <- function(..., print = FALSE){
if(print){print(result)}
Eval(result, list(...))}
return(foo)
}
# create the function
foo <- myTaylor('Exp(x)', 3, 'x')
foo(x=1:5)
## [1] 2.666667 6.333333 13.000000 23.666667 39.333333
foo(x=1:5, print = TRUE)
## expression(x + x^2/2 + x^3/6 + 1)
## [1] 2.666667 6.333333 13.000000 23.666667 39.333333
于 2013-02-08T03:26:38.660 回答
5
将上述符号解与数值泰勒近似进行比较:
library(pracma)
p <- taylor(f = exp, x0 = 0, n = 4) # Numerical coefficients
# 0.1666667 0.5000000 1.0000000 1.0000000 # x^3/6 + x^2/2 + x + 1
polyval(p, 1:5) # Evaluate the polynomial
# 2.66667 6.33333 13.00000 23.66667 39.33334 # exp(x) at x = 1:5
于 2013-02-08T08:47:17.107 回答
1
如果您在误差传播的背景下要求泰勒近似,您可以尝试我的包的“传播”函数,它以(著名的误差传播矩阵表示法)qpcR
的形式评估符号梯度向量和协方差矩阵,g * V * t(g)
相当于一阶泰勒展开。
于 2013-02-09T09:55:37.410 回答
1
作为后续,请考虑:
foo <- myTaylor('Exp(x)', 3, 'x')
sprintf('%2.15f',foo(x=1:5))
[1] "2.666666666666667" "6.333333333333333" "13.000000000000000"
[4] "23.666666666666664" "39.333333333333329"
p <- taylor(f = exp, x0 = 0, n = 3)
sprintf('%2.15f',polyval(p,1:5))
[1] "2.666666721845557" "6.333333789579300" "13.000001556539996"
[4] "23.666670376066413" "39.333340601497312"
其中哪一个更准确,我将留给读者:-)
于 2013-02-08T17:54:16.287 回答
0
下面的函数返回通过在点 a 处使用函数 f 的 n 阶的泰勒级数逼近获得的函数。
taylor <- function(f, n, a) {
ith_derivative <- as.expression(body(f))
f_temp <- function(x) x
series <- as.character(f(a))
for (i in seq_len(n)) {
ith_derivative <- body(f_temp) <- D(ith_derivative, "x")
series <- paste0(series, "+", f_temp(a) / factorial(i), "*(x - ", a, ")^", i)
}
f_output <- function(x) x
body(f_output) <- parse(text = series)
f_output
}
taylor(f = function(x) sin(x), n = 3, a = 0)
于 2021-10-21T15:55:15.410 回答