函数知识

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} f\left(x_{k}\right)=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} f\left(x_{k}\right)+4 \sum_{k=0}^{n-1} f\left(x_{k}+\frac{h}{2}\right)\right] \]

全部代码

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 \right)}{x}\)（由 func1 <- function(x) sin(x) / x 可知）

h <- (b - a) / n

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

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

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

f_x <- sapply(add_by, func)

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

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

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

return(x)

返回结果 x

代码逐步分析 -- Simpson

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

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

h <- (b - a) / n

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

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

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

# 偶数项
add_by_2 <- add_by_1 + h

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

add_by_2 <- add_by_2[-length(add_by_2)]

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

\[add\_by\_2\quad序列去掉\quad\frac{\left( a+b \right)}{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 \right)}{a}+4*\left( \frac{\sin \left( a+h \right)}{a+h}+\frac{\sin \left( a+3h \right)}{a+3h}+...+\frac{\sin \left( b-h \right)}{b-h} \right) +2*\left( \frac{\sin \left( a+2h \right)}{a+2h}+\frac{\sin \left( a+4h \right)}{a+4h}+...+\frac{\sin \left( b \right)}{b} \right) +\frac{\sin \left( b \right)}{b} \right) \]

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

return(x)

返回结果 x