R 语言案例

函数知识

seq 函数

用于生成一段步长相等的序列，看例子即可理解

> seq(5)
[1] 1 2 3 4 5

> seq(2,5)
[1] 2 3 4 5

> seq(2,10,2)
[1] 2 4 6 8 10

sapply 函数

将列表或向量作为输入，并以向量或矩阵形式输出

# 先创建一个函数，这个函数的作用是让 x 除以 2
func1 <- function(x) x / 2

# 用 c() 函数创建向量，并赋值给 nums
nums <- c(10,8,6,4,2)

# sappky(x,y) 的 x 代表向量，y 代表函数
result <- sapply(nums,func1)

# 结果就是生成一个全部除以 2 了的序列
[1] 5 4 3 2 1

length 函数

# 第一个是一个普通的例子
> temp <- c(1,2,3,4,5)
> length(temp)
[1] 5

# 下面的例子比较特殊
> temp <- c(1,2,3,4,5)
> temp[2]
[1] 2
> temp[length(temp)]
[1] 5
# 这里相当于向量 temp 除去了 5，并打印剩下的序列
> temp[-length(temp)]
[1] 1 2 3 4

梯形积分法 与 Simpson

公式

复化梯形公式

若将区间 ([a, b] n) 等分, 有 (n+1) 个等距结点 (x_{k}), 在每一个小区间采用梯形公式可以得到：

[int_{a}^{b} f(x) mathrm{d} x=sum_{k=0}^{n-1} int_{x_{k}}^{x_{k+1}} f(x) mathrm{d} x approx frac{h}{2}[f(a)+f(b)]+h sum_{k=1}^{n-1} fleft(x_{k} ight)=T_{n} ]

复化辛卜生公式

若将区间 ([a, b] n) 等分, 有 (n+1) 个等距结点 (x), 在每一个小区间采用辛卜生公式可以得到：

[int_{a}^{b} f(x) mathrm{d} x=sum_{k=0}^{n-1} int_{x_{k}}^{x_{k+1}} f(x) mathrm{d} x approx frac{h}{6}left[f(a)+f(b)+2 sum_{k=1}^{n-1} fleft(x_{k} ight)+4 sum_{k=0}^{n-1} fleft(x_{k}+frac{h}{2} ight) ight] ]

全部代码

func1 <- function(x) sin(x) / x

# 梯形积分法
Trapezoid <- function(func, a, b, n = 100) {
    h <- (b - a) / n
    add_by <- seq(a, b, by = h)
    f_x <- sapply(add_by, func)
    x <- h * sum(f_x[1] / 2, f_x[2:n], f_x[n + 1] / 2)
    return(x)
}

# Simpson
Simpson <- function(func, a, b, n = 100) {
    h <- (b - a) / n
    # 奇数项
    add_by_1 <- seq(a + h, b - h, by = 2 * h)
    # 偶数项
    add_by_2 <- add_by_1 + h
    add_by_2 <- add_by_2[-length(add_by_2)]
    x <- h / 3 * sum(func(a), 4 * sapply(add_by_1, func), 2 * sapply(add_by_2, func), func(b))
    return(x)
}

代码逐步分析 -- 梯型积分法

Trapezoid <- function(func, a, b, n = 100)

Trapezoid 是梯型（函数名字），传入的函数 func 是 (frac{sin left( x ight)}{x})（由 func1 <- function(x) sin(x) / x 可知）

h <- (b - a) / n

[h=frac{left( b-a ight)}{n} ]

add_by <- seq(a, b, by = h)

[add\_by = aquad a+hquad a+2hquad ...quad b ]

f_x <- sapply(add_by, func)

[f\_x = frac{sin left( a ight)}{a}quad frac{sin left( a+h ight)}{a+h}quad frac{sin left( a+2h ight)}{a+2h}quad ...quad frac{sin left( b ight)}{b} ]

x <- h * sum(f_x[1] / 2, f_x[2:n], f_x[n + 1] / 2)

[x=h*left( frac{sin left( a ight)}{2a}+frac{sin left( a+h ight)}{a+h}+frac{sin left( a+2h ight)}{a+2h}+...+frac{sin left( n ight)}{n}+frac{sin left( n+1 ight)}{2left( n+1 ight)} ight) ]

return(x)

返回结果 x

代码逐步分析 -- Simpson

Simpson <- function(func, a, b, n = 100)

Simpson 是辛卜生（函数名字），传入的函数 func 是 (frac{sin left( x ight)}{x})（由 func1 <- function(x) sin(x) / x 可知）

h <- (b - a) / n

[h=frac{left( b-a ight)}{n} ]

# 奇数项
add_by_1 <- seq(a + h, b - h, by = 2 * h)

[add\_by\_1 = a+hquad a+3hquad a+5hquad ...quad b-h ]

# 偶数项
add_by_2 <- add_by_1 + h

[add\_by\_2 = a+2hquad a+4hquad a+6hquad ...quad b ]

add_by_2 <- add_by_2[-length(add_by_2)]

[length(add\_by\_2)quad代表quad add\_by\_2quad的长度，你可以思考为什么等于quadfrac{left( a+b ight)}{4h} ]

[add\_by\_2quad序列去掉quadfrac{left( a+b ight)}{4h} ]

x <- h / 3 * sum(func(a), 4 * sapply(add_by_1, func), 2 * sapply(add_by_2, func), func(b))

[x=frac{h}{3}*left( frac{sin left( a ight)}{a}+4*left( frac{sin left( a+h ight)}{a+h}+frac{sin left( a+3h ight)}{a+3h}+...+frac{sin left( b-h ight)}{b-h} ight) +2*left( frac{sin left( a+2h ight)}{a+2h}+frac{sin left( a+4h ight)}{a+4h}+...+frac{sin left( b ight)}{b} ight) +frac{sin left( b ight)}{b} ight) ]

最后的结果和公式不太一致，你可以思考为什么可以这样做？【提示：(frac{h}{2})】

return(x)

返回结果 x