R计算微积分

r
符号计算
微积分
作者

zsc

发布日期

2018年2月10日

1.1差分

Show the code
x=1:12
diff(x)  #向量差分  后面一个数减去前面一个数
 [1] 1 1 1 1 1 1 1 1 1 1 1
Show the code
z=matrix(x,3,4)
z
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
Show the code
diff(z) #矩阵差分 前行减去后行
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    1    1

1.2 符号计算–微分

1.2.1求一元函数导数— \(\sin{x}\) 的一阶导数为： \(\cos{x}\)

在R里，声明表达式对象使用 expression() 函数, 计算一阶导数用D()函数,格式：D(表达式,对谁求导)

Show the code
fun=expression(sin(x))# 声明表达式
D(fun,"x")#---方法1
cos(x)
Show the code
deriv(fun,"x")#---方法2  其中.grad[, "x"]为求x的导数表达式
expression({
    .value <- sin(x)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- cos(x)
    attr(.value, "gradient") <- .grad
    .value
})
Show the code
deriv3(fun,"x")#---方法3 其中.grad[, "x"]为求x的导数表达式
expression({
    .expr1 <- sin(x)
    .value <- .expr1
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .hessian <- array(0, c(length(.value), 1L, 1L), list(NULL, 
        c("x"), c("x")))
    .grad[, "x"] <- cos(x)
    .hessian[, "x", "x"] <- -.expr1
    attr(.value, "gradient") <- .grad
    attr(.value, "hessian") <- .hessian
    .value
})
1.2.2 计算n阶导数

  1. 计算n阶导数—方法一: 结合一阶导数写递归函数

    函数： \(\sin{x}+\cos{2x}+x^2+xy+y^2+2x^3+y^3\) 的3阶导数为：\(12 + 8\sin{2x} -\cos{x}\)

Show the code
   fun=expression(sin(x)+cos(2*x)+x^2+x*y+y^2+2*x^3+y^3)
   DD <- function(expr, name, order = 1) {
     if (order < 1)
       stop("'order' must be >= 1")
     if (order == 1){
       D(expr, name)
     }else{
       DD(D(expr, name), name, order - 1)
     }
   }
   DD(fun,"x",3)
2 * (3 * 2) - (cos(x) - sin(2 * x) * 2 * 2 * 2)
  1. 计算n阶导数—方法二: Deriv 包中Simplify()化简表达式
Show the code
   library(Deriv)
   DD(fun,"x",3)
2 * (3 * 2) - (cos(x) - sin(2 * x) * 2 * 2 * 2)
Show the code
   Simplify(DD(fun, "x", 3))
12 + 8 * sin(2 * x) - cos(x)
1.2.3 通过函数计算导数

有时候我们有的就是函数，这怎么计算导数呢？—还是用上面的函数

Show the code
f=function(x,y)sin(x)+cos(2*x)+x^2+x*y+y^2+2*x^3+y^3 #这里是函数，而不是表达式
body(f)
sin(x) + cos(2 * x) + x^2 + x * y + y^2 + 2 * x^3 + y^3
Show the code
Simplify(D(body(f), "x"))# 注意：函数体有花括号{}会出错
cos(x) + x * (2 + 6 * x) + y - 2 * sin(2 * x)
1.2.4 求二元函数偏导数及梯度
Show the code
D(expression(x^2+x*y+y^2),"x")# x偏导数
2 * x + y
Show the code
D(expression(x^2+x*y+y^2),"y")# y偏导数
x + 2 * y
1.2.5 符号计算扩展包 Ryacas

想要做更多的符号计算内容，如解方程，泰勒展开等，可以借助第三方 R 扩展包 Ryacas

解方程: \(\frac{x}{1+x}=a\) 求解\(x=\frac{a}{1-a}\)

Show the code
library(Ryacas)
ysym("Solve(x/(1+x) == a, x)") 
{x==a/(1-a)}

多项式展开：如\((1+x)^3\) 展开

Show the code
ysym(expression(Expand((1 + x)^3)))# 把(1+x)^3展开
y: x^3+3*x^2+3*x+1

求解常微分方程：\(y''=4y\)

Show the code
ysym("OdeSolve(y''== 4 * y)")
y: C105*Exp(2*x)+C109*Exp((-2)*x)

泰勒展开：\(\cos{x}=1-\frac{1}{2!}x^2+\frac{1}{4!}x^4+o(x^4)\)

Show the code
ysym("Taylor(x, a, 3) cos(x)") # cos(x)函数在a点的3阶泰勒展开
y: cos(a)+(Deriv(a)cos(a))*(x-a)+((x-a)^2*(Deriv(a)Deriv(a)cos(a)))/2+((x-a)^3*(Deriv(a)Deriv(a)Deriv(a)cos(a)))/6
1.3 表达式转为函数值

很多时候我们使用 R 目的是计算，符号计算后希望可以直接代入计算，那么只需要在 deriv 中指定 function.arg 参数为 TRUE。 \[ \sin{x}+\cos{2x}+x^2+xy+y^2对x求偏导为：\cos{x}-2sin{2x}+2x+y \]

Show the code
fun=expression(sin(x)+cos(2*x)+x^2+x*y+y^2)
dx=D(fun,"x") #用D()函数得到符号运算结果，然后代入数值即可得到最后结果
dx
cos(x) - sin(2 * x) * 2 + 2 * x + y
Show the code
x=0;y=pi# 对x、y赋值
eval(dx)#求出数值解
[1] 4.141593
Show the code
Dfun=deriv(fun,c("x","y"),function.arg = TRUE)# 同时对x、y求偏导--D()函数不可以同时求偏导
Dfun(x=0,y=pi/2) # 代值计算，其中attr(,"gradient")的值为导数值 ，另一个为原函数在该处的函数值
[1] 3.467401
attr(,"gradient")
            x        y
[1,] 2.570796 3.141593
Show the code
#我们可以作如下简单验证：
fun=function(x,y){sin(x)+cos(2*x)+x^2+x*y+y^2}
fun(x=0,y=pi/2)
[1] 3.467401

1.3、求积分—暂时只找到数值计算的–没找到符号计算的

积分函数： integrate(fun,a,b) fun被积函数，不需要表达式，因为这是数值计算, a,b为上下限

Show the code
f <- function (x) sin(x)
integrate(f,0,pi/2)
1 with absolute error < 1.1e-14

有些时候只想要值输入 integrate(f,0,1)$value

Show the code
integrate(f,0,pi/2)$value
[1] 1
Show the code 
R version 4.2.1 (2022-06-23)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.5.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Ryacas_1.1.3.1 Deriv_4.1.3   

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9        digest_0.6.29     jsonlite_1.8.0    magrittr_2.0.3   
 [5] evaluate_0.16     rlang_1.0.4       stringi_1.7.8     cli_3.3.0        
 [9] rstudioapi_0.14   rmarkdown_2.16.1  tools_4.2.1       stringr_1.4.1    
[13] htmlwidgets_1.5.4 pkgload_1.3.0     xfun_0.32         yaml_2.3.5       
[17] fastmap_1.1.0     compiler_4.2.1    htmltools_0.5.3   knitr_1.40