[学习笔记] - 数值积分

我们在微积分中，通过切割面积，使用无数个矩形去逼近曲线的面积，然后求极限得到曲线的面积，但是通过上面的图可以发现，这样太慢了，其他问题还有

• $f(x)$形式复杂，求原函数非常困难

$f(x) = \sqrt{ax^2 + bx + c}$

• $f(x)$的原函数不能用初等函数形式表示

$f(x) = \frac{1}{\ln x}, \quad e^{-x^2}, \quad \sin x^2, \quad \frac{\sin x}{x}$

• $f(x)$虽有初等函数表示的原函数，但其原函数表示形式相当复杂

$f(x) = \frac{1}{1 + x^4}$

• $f(x)$本身没有解析表达式，其函数关系由表格或图像给出，如实验或测量数据。

牛顿-柯斯特公式

在数值分析上，梯形法则和辛普森法则都是数值积分的方法，用于计算定积分
这两种方法都属于牛顿-柯斯特公式，使函数于等距$n + 1$点的值，取得一个$n$次的多项式来近似原来的函数，再行求积

积分中值定理

对于连续函数$f(x)$，有

$I = \int^b_a f(x)dx = (b - a)f(c), \quad c \in [a, b]$

若能对$f(c)$提供一种近似，就能够得到对应的数值积分公式

$\begin{aligned} I &= \int^b_a f(x)dx \approx (b -a)f(a) & \text{左矩形公式} \\ I &= \int^b_a f(x)dx \approx (b -a)f(b) & \text{右矩形公式} \\ I &= \int^b_a f(x)dx \approx (b -a)f(\frac{a+b}{2}) & \text{中矩形公式} \end{aligned}$

左矩形公式	右矩形公式	中矩形公式

设$f(x)$在$[a, b]$内$n + 1$个节点$x_i$处的高度为$f(x_i)$，可通过加权平均的方法近似地得到平均高度$f(c)$，从而得到求积公式

$I = \int^b_a f(x)dx \approx \sum^n_{i=0}A_i f(x_i)$

• $x_i$为求积节点
• $A_i$为求积系数

称

$R_n(f) = \int^b_a f(x) dx - \sum^n_{i=0}A_i f(x_i)$

为求积公式的截断误差
以中矩形公式为例，设$F(x)$为$f(x)$的原函数

$\int^x_a f(s) ds = F(x)$

由泰勒展开可知

$\begin{aligned} \int^{a+h}_a f(s) ds &= F(a + h) \\ &= F(a) + F'(a)h + F''(a)\frac{h^2}{2} + F'''(a)\frac{h^3}{6} + O(h^4), \quad h = b - a \end{aligned}$

由$F(x)$的定义可知

$\begin{aligned} F(a) &= 0 \\ F'(a) &= f(a) \\ F''(a) &= f'(a) \\ F'''(a) &= f''(a) \end{aligned}$

从而

$\int^{a+h}_a f(s) ds = 0 + f(a)h + f'(a)\frac{h^2}{2} + f''(a)\frac{h^3}{6} + O(h^4)$

而

$h \cdot f(a+ \frac{h}{2}) = h \cdot \Big[f(a) + f'(a)\frac{h}{2} + f''(a)\frac{h^2}{8} + O(h^3)\Big]$

于是

$\Big| \int^{a+h}_a f(s) ds - h \cdot f(a+\frac{h}{2}) \Big| = \frac{h^3}{24} f''(a) + O(h^4) = O(h^3)$

插值求积公式

值是建立逼近函数的手段，用以研究原函数的性质。因此，可以用插值函数的导数近似为原函数的导数

$f^{(k)}(x) \approx L^{(k)}_n{x}$

设$[a, b]$上的节点为$a = x_0 < x_1 < \cdot < x_n < b$，则$f(x)$的拉格朗日插值多项式为

$L_n(x) = \sum^n_{i=0}l_i(x)f(x_i), \quad l_i(x) = \prod^n_{j=0,j\neq i}\frac{x - x_j}{x_i - x_j}$

用$L_n(x)$作为$f(x)$的近似函数有

$\begin{aligned} I &= \int^b_a f(x) dx \approx \int^b_a L_n(x) dx \\ &= \int^b_a \sum^n_{i=0} l_i(x) f(x_i) dx \\ &= \sum^n_{i=0} (\int^b_a l_i(x)dx) f(x_i) \end{aligned}$

记$A_i = \int^b_a l_i(x)dx$，则有插值性求积公式

$\color{red} I = \int^b_a f(x) dx \approx \sum^n_{i = 0} A_i f(x_i)$

$A_i$只与插值节点$x_i$有关，与被积函数$f(x)$无关

上述公式的截断误差为

$\begin{aligned} R_n(f) &= \int^b_a \big[f(x) - L_n(x)\big]dx \\ &= \frac{1}{(n + 1)!}\int^b_a f^{(n+1)}(c) w_{n+1}(x)dx \\ \\ w_{n+1} &= \prod^n_{i=0}(x - x_i), \quad c \in (a, b) \end{aligned}$

Cotes系数

在插值型求积公式中，取等距节点，即将$[a,b]$作$n$等分，记节点为$x_i = a + ih(i = 1, 2, \cdots, n), h = \frac{b - a}{n}$。记$x = a+ th$，则

$\begin{aligned} A_i &= \int^b_a \prod^n_{j = 0, j \neq i}\frac{x - x_j}{x_i - x_j}dx \\ &= h\int^n_0 \prod^n_{j = 0, j \neq i}\frac{t - j}{i - j}dt \\ &= \frac{b - a}{n} \cdot \frac{(-1)^{n-i}}{i!(n - i)!} \int^n_0\prod^n_{j = 0, j \neq i}(t - j)dt \end{aligned}$

记

$C^{(n)}_i = \frac{1}{n}\frac{(-1)^{n-i}}{i!(n - i)!} \int^n_0\prod^n_{j = 0, j \neq i}(t - j)dt \Rightarrow \color{red} A_i = (b -a)C^{(n)}_i$

于是得到牛顿-柯斯特公式

$\color{red} I = \int^b_a f(x) dx \approx (b - a)\sum^n_{i=0} C^{(n)}_i f(x_i)$

$C^{(n)}_i$成为Cotes系数

• 对称性：$C^{(n)}_i = C^{(n)}_{n-i}$
• $\sum^n_{i=1}C^{(n)}_i = 1$

梯形公式

$\begin{aligned} C^{(1)}_0 &= C^{(1)}_1 = \int^1_0 t dt = \frac{1}{2} \\ \color{red}I &\color{red}= \int^b_a f(x) dx \approx \frac{b - a}{2}\Big [f(a) + f(b) \Big] \end{aligned}$

若$f''(x)$在$[a, b]$上连续，则梯形公式的截断误差为

$R_1(f) = -\frac{(b - a)^3}{12}f''(c), \quad c \in [a,b]$

辛普森公式

$n = 2, Simpson rule$

$\begin{aligned} C^{(2)}_0 &= C^{(2)}_2 = \frac{1}{4}\int^2_0 t(t - 1)dt = \frac{1}{6} \\ C^{(2)}_1 &= \frac{1}{4} \int^2_0 t(t - 2)dt = \frac{4}{6} \\ \\ \color{red}I &\color{red}= \int^b_a f(x) dx \approx \frac{b - a}{6}\Big[f(a) + 4f(\frac{a+b}{2}) + f(b)\Big] \end{aligned}$

若$f^{(4)}(x)$在$[a, b]$上连续，则辛普森公式的截断误差为

$\begin{aligned} R_2(f) &= -\frac{(b - a)^5}{2880}f^{(4)}(c) \\ &= -\frac{1}{90}(\frac{b - a}{2})^5 f^{(4)}(c), \quad c \in [a, b] \end{aligned}$

辛普森3/8公式

$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{8}\Big[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)\Big]$

柯斯特公式

$n = 4, Cotes rule$

$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{90}\Big[7f(x_0) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(x_4)\Big]$

若$f^{(6)}(x)$在$[a, b]$上连续，则柯特斯公式的截断误差为

$\begin{aligned} R_4(f) &= -\frac{(b - a)^7}{1013760}f^{(4)}(c) \\ &= -\frac{8}{495}(\frac{b - a}{2})^7 f^{(6)}(c), \quad c \in [a, b] \end{aligned}$

总结

方法	公式
梯形公式	$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{2}\Big [f(a) + f(b) \Big]$
辛普森公式	$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{6}\Big[f(a) + 4f(\frac{a+b}{2}) + f(b)\Big]$
辛普森3/8公式	$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{8}\Big[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)\Big]$
柯斯特公式	$\color{red}I = \int^b_a f(x) dx \approx \frac{b - a}{90}\Big[7f(x_0) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(x_4)\Big]$

参考

• 数值积分(Quadrature)