R语言学习之科学计算——求导与积分

说明

R语言的版本为4.0.2，IDE为Rstudio，版本为1.3.959。学习过程中参考了以下文章：
R笔记（1）：formula和Formula（CSDN）
R语言笔记.formula（知乎）
R语言进阶之五：表达式、数学公式与特殊符号（CSDN）
R语言的导数计算（CSDN）
积分_R（CSDN）

1 R中的表达式

要进行求导或者积分运算，首先需要有一个表达式，注意，是表达式，而不是实现了这个表达式的函数，这就要用到 expression 对象。
使用expression() 函数可以创建expression 对象，expression 对象实际上是以列表的形式储存表达式的，例如创建以下表达式：

$$y=x^3+sin(x)+e^x$$

当然除了使用 expression() 函数创建表达式之外，还可以使用formula() 函数创建公式，二者都可以当成参数进行求导和积分运算。

# f1,f2,f3分别创建了表达式、公式和函数
# 区别是函数可以传参并运行，前面两种不行
f1  <-  expression(y=x^3+sin(x)+exp(x))
f2 <- formula(y~x^3+sin(x)+exp(x))
f3 <- function(x) x^3 + sin(x) + exp(x)
mode(f1)
mode(f2)
mode(f3)

# 输出
# > mode(f1)
# [1] "expression"
# > mode(f2)
# [1] "call"
# > mode(f3)
# [1] "function"

2 求导

求导可以使用函数deriv() 来实现，可以用它来求

1. 一阶导数
2. 高阶导数
3. 偏导数等
4. ……

2.1 求一阶导数

用于求导的函数表达式为：

$$y=x^3+sin(x)+e^x$$

求一阶导数，可以直接使用deriv() 函数，或者使用D()函数，只是后者只支持传入expression不能传入formula，也不能设置让其返回函数用于计算，以下例子里面仅仅使用了deriv() 函数。

# f1,f2,f3分别创建了表达式、公式和函数
# 区别是函数可以传参并运行，前面两种不行
f1  <-  expression(y=x^3+sin(x)+exp(x))
f2 <- formula(y~x^3+sin(x)+exp(x))
f3 <- function(x) x^3 + sin(x) + exp(x)
mode(f1)
mode(f2)
mode(f3)

# 分别使用表达式、公式和手动编求导后的函数
# 验证求导结果，发现结果一致
df1 <- deriv(f1,"x",function.arg = TRUE)
df2 <- deriv(f2,"x",function.arg = TRUE)
df3 <- function(x) 3*x^2 + cos(x) + exp(x)
df1(1)
df2(1)
df3(1)

# 输出
# > df1(1)
# [1] 4.559753
# attr(,"gradient")
# x
# [1,] 6.258584
# > df2(1)
# [1] 4.559753
# attr(,"gradient")
# x
# [1,] 6.258584
# > df3(1)
# [1] 6.25858

2.2求高阶导数

用于求导的函数表达式依旧为：

$$y=x^3+sin(x)+e^x$$

求二阶导数可直接使用deriv3() 函数。

# f1,f2,f3分别创建了表达式、公式和函数
# 区别是函数可以传参并运行，前面两种不行
f1  <-  expression(y=x^3+sin(x)+exp(x))
f2 <- formula(y~x^3+sin(x)+exp(x))
f3 <- function(x) x^3 + sin(x) + exp(x)
mode(f1)
mode(f2)
mode(f3)

# 分别使用表达式、公式和手动编求二阶导数后的函数
# 验证求导结果，发现结果一致
df1 <- deriv3(f1,"x",function.arg = TRUE)
df2 <- deriv3(f2,"x",function.arg = TRUE)
df3 <- function(x) 6*x - sin(x) + exp(x)
df1(1)
df2(1)
df3(1)

# 输出
# > df1(1)
# [1] 4.559753
# attr(,"gradient")
# x
# [1,] 6.258584
# attr(,"hessian")
# , , x
# 
# x
# [1,] 7.876811
# 
# > df2(1)
# [1] 4.559753
# attr(,"gradient")
# x
# [1,] 6.258584
# attr(,"hessian")
# , , x
# 
# x
# [1,] 7.876811
# 
# > df3(1)
# [1] 7.876811

求更高阶导数则需要自己编一个函数，方便高阶求导

# f1,f2,f3分别创建了表达式、公式和函数
# 区别是函数可以传参并运行，前面两种不行
f1  <-  expression(y=x^3+sin(x)+exp(x))
f2 <- formula(y~x^3+sin(x)+exp(x))
f3 <- function(x) x^3 + sin(x) + exp(x)
mode(f1)
mode(f2)
mode(f3)

# 自定义的求高阶导数的函数，设置order值以确定求导阶数
DD <- function(expr, name, order = 1,function.arg = FALSE) {
  if(order < 1) stop("'order' must be >= 1")
  if(order == 1) deriv(expr, name, function.arg = function.arg)
  else DD(D(expr, name), name, order - 1, function.arg)
}

# 分别使用表达式、公式和手动编求3阶导数后的函数
# 验证求导结果，发现结果一致
df1 <- DD(f1,"x",3,function.arg = TRUE)
df3 <- function(x) 6 - cos(x) + exp(x)
df1(1)
df3(1)

# 输出
# [1] 7.876811
# attr(,"gradient")
# x
# [1,] 8.17798
# > df3(1)
# [1] 8.17798

2.3 求偏导数

对以下函数求偏导数：
$$f(x,y)=\frac{x^3}{y}+sin(x{y}^3)+e^{xy}$$

f1  <-  expression(f=x^3/y+sin(x*y^3)+exp(x*y))
f2 <- function(x,y) x^3/y+sin(x*y^3)+exp(x*y)

# 分别使用表达式和手动编求偏导数后的函数
# 验证求导结果，发现结果一致
df1xy <- deriv(f1,c("x","y"),function.arg = TRUE)
df2x <- function(x,y) 3*x^2/y+cos(x*y^3)*y^3+y*exp(x*y)
df2y <- function(x,y) -x^3/y^2+3*cos(x*y^3)*y^2+x*exp(x*y)
df1xy(1,1)
df2x(1,1)
df2y(1,1)

# 输出
# [1] 4.559753
# attr(,"gradient")
# x        y
# [1,] 6.258584 3.339189
# > df2x(1,1)
# [1] 6.258584
# > df2y(1,1)
# [1] 3.339189

3 积分

对于积分而言，有以下几种形式：

1. 定积分
2. 不定积分
3. 多元函数积分
4. ……

3.1 定积分

求定积分可以使用integrate() 函数，该函数要求传入的是仅有单变量的函数而非表达式，因此只能求一元函数定积分，对以下函数求定积分：

$$y=x^3+sin(x)+e^x$$

# 定义一个用于积分的单变量函数，和积分之后的函数
f  <-  function(x) x^3 + sin(x) + exp(x)
If <- function(x) x^4/4-cos(x) +exp(x)

# 计算区间（1，2）上的积分
# 并手动验证，发现一致，误差很小
integrate(f, 1,2)
If(2)-If(1)

# 输出
# > integrate(f, 1,2)
# 9.377223 with absolute error < 1e-13
# > If(2)-If(1)
# [1] 9.377223

3.2 不定积分

暂时还不知道和不定积分相关的包有哪些，日后找到再补上。