函数知识
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