CS205A Lecture 15 Numerical integration and differentiation

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第十五讲，这一讲介绍了数值积分和微分。

Chapter 12 数值积分和微分

积分

首先讨论一元函数求积分的问题，即

$$f:\mathbb R\to \mathbb R$$

我们的目标是求

$$I(f)=\int_a^b f(x)dx$$

这里有两个假设：

1. $a,b$固定。
2. 可以计算$f(x)$任意一点的值。

插值求积

因为定积分是由黎曼和定义的，即

$$\int_{a}^{b} f(x) dx=\lim _{\Delta x_{k} \rightarrow 0} \sum_{k} f\left(\tilde{x}_{k}\right)\left(x_{k+1}-x_{k}\right)$$

其中

$$\begin{aligned} &1.a=x_{1}<x_{2}<\cdots<x_{n}=b\\ &2.\Delta x_k =x_{k+1}-x_k\\ &3.\tilde{x}_{k}\in [x_k,x_{k+1}] \end{aligned}$$

所以一种常见的思路是运用下式近似计算$\int_{a}^{b} f(x)$

$$Q_n[f] \triangleq\sum_{i=1}^n w_{i} f\left(x_{i}\right)$$

其中$w_i$为权重。类似插值的思想，我们希望对于多项式基$f_k(x)= x^k,k=0,\ldots n-1 $，有如下事实成立

$$\int_a^b x^kdx=\sum_{i=1}^n w_{i} x_i^k$$

即

$$\begin{array}{l} {w_{1} \cdot 1+w_{2} \cdot 1+\cdots+w_{n} \cdot 1=\int_{a}^{b} 1 {d} x=b-a} \\ {w_{1} \cdot x_{1}+w_{2} \cdot x_{2}+\cdots+w_{n} \cdot x_{n}=\int_{a}^{b} x {d} x=\left(b^{2}-a^{2}\right) / 2}\\ \vdots\\ {w_{1} \cdot x_{1}^{n-1}+w_{2} \cdot x_{2}^{n-1}+\cdots+w_{n} \cdot x_{n}^{n-1}=\int_{a}^{b} x^{n-1} {d} x=\left(b^{n}-a^{n}\right) / n} \end{array}$$ $$\left( \begin{array}{ccccc} {1} & {1} & {\cdots} & 1\\ {x_{1}}& {x_{2}} &\cdots & x_n\\ \vdots &\vdots & & \vdots \\ {x_{1}^{n-1}} & {x_{2}^{n-1}} &\ldots & {x_{n}^{n-1}} \end{array}\right) \left( \begin{array}{c} {w_{1}} \\ {w_{2}} \\ {\vdots} \\ {w_{n}} \end{array}\right)= \left( \begin{array}{c}{b-a} \\ {\left(b^{2}-a^{2}\right) / 2}\\ {\vdots} \\ {\left(b^{n}-a^{n}\right) / n} \end{array}\right)$$

系数矩阵为范德蒙矩阵，求解该方程即可，来看一个具体例子

例1

考虑$n=3$的情形：

$$Q_{3}(f)=w_{1} f\left(x_{1}\right)+w_{2} f\left(x_{2}\right)+w_{3} f\left(x_{3}\right)$$

其中

$$\begin{aligned} x_{1} &=a \\ x_{2} &=\frac{a+b}{2} \\ x_{3} &=b \end{aligned}$$

线性方程组为

$$\left(\begin{array}{ccc}{1} & {1} & {1} \\ {a} & {(a+b) / 2} & {b} \\ {a^{2}} & {((a+b) / 2)^{2}} & {b^{2}}\end{array}\right) \left( \begin{array}{c}{w_{1}} \\ {w_{2}} \\ {w_{3}}\end{array}\right)=\left( \begin{array}{c}{b-a} \\ {\left(b^{2}-a^{2}\right) / 2} \\ {\left(b^{3}-a^{3}\right) / 3}\end{array}\right)$$

解得

$$\begin{aligned} w_{1}&=\frac{b-a} 6\\ w_{2}&=\frac{2(b-a) } 3\\ w_{3}&=\frac{b-a} 6 \end{aligned}$$

即

$$Q_{3}(f)=\frac{b-a}{6}\left(f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right)$$

这个公式被称为辛普森公式。

Newton-Cotes求积

这部分讨论节点$\{x_1,\dots ,x_n\}$的选择问题，最常见的节点选择方式为选择等距节点，Newton-Cotes求积就有该特点，分别有如下两种：

1. Closed：包括端点$a,b​$

   $$x_{k} = a+\frac{(k-1)(b-a)}{n-1},k=1,\ldots, n$$

2. Open：不包括端点$a,b$

   $$x_{k} = a+\frac{k(b-a)}{n+1},k=1,\ldots, n$$

对应图片如下：

根据这两种节点选择的方式，产生了许多常用公式。

Open：

取$n=1$，于是产生了中点公式：

$$\int_{a}^{b} f(x) d x \approx(b-a) f\left(\frac{a+b}{2}\right)$$

Closed：

取$n=2$，那么

$$\begin{aligned} x_1 &=a\\ x_2 &=b \end{aligned}$$

于是产生梯形公式：

$$\int_{a}^{b} f(x) d x \approx(b-a) \frac{f(a)+f(b)}{2}$$

取$n=3$，那么

$$\begin{aligned} x_{1} &=a \\ x_{2} &=\frac{a+b}{2} \\ x_{3} &=b \end{aligned}$$

此即为例1的情形，所以

$$\int_{a}^{b} f(x) d x \approx\frac{b-a}{6}\left(f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right)$$

这三种方法的图示如下：

复合求积

实际中我们肯定希望使用多个节点求数值积分，那么相邻两个节点之间可以使用之前介绍的方法：

假设

$$\Delta x = \frac{b-a}{k}, x_{i}= a+i \Delta x$$

那么复合中点公式为

$$\int_{a}^{b} f(x) d x \approx \sum_{i=0}^{k-1} f\left(\frac{x_{i+1}+x_{i}}{2}\right) \Delta x$$

复合梯形公式为

$$\begin{aligned} \int_{a}^{b} f(x) d x &\approx \sum_{i=0}^{k-1}\left(\frac{f\left(x_{i}\right)+f\left(x_{i+1}\right)}{2}\right) \Delta x \\ &=\Delta x\left(\frac{1}{2} f(a)+f\left(x_{1}\right)+\cdots+f\left(x_{k-1}\right)+\frac{1}{2} f(b)\right) \end{aligned}$$

在$n$为偶数的情形下，复合辛普森公式为

$$\begin{aligned} \int_{a}^{b} f(x) d x & \approx \frac{\Delta x}{3}\left[f(a)+2 \sum_{i=1}^{n-2-1} f\left(x_{2 i}\right)+4 \sum_{i=1}^{n / 2} f\left(x_{2 i-1}\right)+f(b)\right] \\ &=\frac{\Delta x}{3}\left[f(a)+4 f\left(x_{1}\right)+2 f\left(x_{2}\right)+\cdots+4 f\left(x_{n-1}\right)+f(b)\right] \end{aligned}$$

精确度

这一部分讨论三种积分公式的精确度，首先定义

$$c=\frac 1 2(a+b)$$

那么$f$在$c$处的泰勒展开为

$$f(x)=f(c)+f^{\prime}(c)(x-c)+\frac{1}{2} f^{\prime \prime}(c)(x-c)^{2}+\frac{1}{6} f^{\prime \prime \prime}(c)(x-c)^{3}+\frac{1}{24} f^{\prime \prime \prime \prime}(c)(x-c)^{4}+\cdots \tag 1$$

对上式积分可得

$$\int_{a}^{b} f(x) d x=(b-a) f(c)+\frac{1}{24} f^{\prime \prime}(c)(b-a)^{3}+\frac{1}{1920} f^{\prime \prime \prime \prime}(c)(b-a)^{5}+\cdots$$

注意中点公式为

$$(b-a) f(c)$$

因此中点公式的精确度为$O\left(\Delta x^{3}\right)​$。

对(1)式令$x=a,b$可得

$$\begin{aligned} f(a) &=f(c)+f^{\prime}(c)(a-c)+\frac{1}{2} f^{\prime \prime}(c)(a-c)^{2}+\frac{1}{6} f^{\prime \prime \prime}(c)(a-c)^{3}+\cdots \\ f(b) &=f(c)+f^{\prime}(c)(b-c)+\frac{1}{2} f^{\prime \prime}(c)(b-c)^{2}+\frac{1}{6} f^{\prime \prime \prime}(c)(b-c)^{3}+\cdots \end{aligned}$$

两式相加并乘以$\frac {b-a}{2}​$得到

$$(b-a) \frac{f(a)+f(b)}{2}=f(c)(b-a)+\frac{1}{4} f^{\prime \prime}(c)(b-a)\left((a-c)^{2}+(b-c)^{2}\right)+\cdots$$

因此梯形公式的精确度同样为$O\left(\Delta x^{3}\right)$。

运用类似的方法可以证明辛普森公式的精确度为$O\left(\Delta x^{5}\right)​$.

注意泰勒展开只有在$\Delta x$充分小的时候成立，所以实际中先要将区间划分为$\frac{b-a}{\Delta x}$个部分，因此实际的误差还要乘以$\frac 1{\Delta x}$，所以误差分别为：

• 中点公式：$O(\Delta x^2)$
• 梯形公式：$O(\Delta x^2)​$
• 辛普森公式：$O(\Delta x^4)$

还有两种求积分的方法为高斯求积以及自适应求积，这部分讲义介绍的比较简单，这里略过。

多变量积分

由于维度灾难问题，如果函数为

$$f:\mathbb R^k \to \mathbb R$$

那么按照之前介绍的方法，我们取点的数量为$O(n^k)​$，所以之前的方法很难推广到多变量的情形，此时我们一般采用蒙特卡洛方法。

条件数

之前讨论的都是精确度问题，这部分条件该方法的条件数。假设$\hat f$是$f$的扰动，那么

$$\begin{aligned} \left|Q_{n}(\hat{f})-Q_{n}(f)\right| &=\left|Q_{n}(\hat{f}-f)\right| \\ &=\left|\sum_{i=1}^{n} w_{i}\left(\hat{f}\left(x_{i}\right)-f\left(x_{i}\right)\right)\right| \\ & \le \sum_{i=1}^{n}\left(\left|w_{i}\right| \cdot\left|\hat{f}\left(x_{i}\right)-f\left(x_{i}\right)\right|\right) \\ &\le \left(\sum_{i=1}^{n}\left|w_{i}\right|\right)\|\hat{f}-f\|_{\infty} \end{aligned}$$

回顾插值求积部分的公式，我们有

$$\sum_{i=1}^n w_i =b-a$$

如果所有权重都非负，那么条件数为

$$\sum_{i=1}^n |w_i| =b-a$$

如果某些权重是负数，那么某些$|w_i|$可能非常大，求积公式就会不太稳定。

微分

这部分讨论如何计算数值微分。

有限差分

回顾泰勒展开：

$$f(x+h)=f(x)+f^{\prime}(x) h+\frac{1}{2} f^{\prime \prime}(x) h^{2}+\cdots$$

调整该等式，我们有

$$f^{\prime}(x)=\frac{f(x+h)-f(x)}{h}+O(h)$$

由此得到前向差分公式：

$$f^{\prime}(x)\approx \frac{f(x+h)-f(x)}{h}$$

同理可得后向差分公式：

$$f^{\prime}(x)\approx \frac{f(x)-f(x-h)}{h}$$

注意这两种方法的精确度都为$O(h)$。

现在考虑如下两个展开：

$$\begin{eqnarray*} f(x+h) &=f(x)+f^{\prime}(x) h+\frac{1}{2} f^{\prime \prime}(x) h^{2}+\frac{1}{6} f^{\prime \prime \prime}(x) h^{3}+\cdots \tag 1\\ f(x-h) &=f(x)-f^{\prime}(x) h+\frac{1}{2} f^{\prime \prime}(x) h^{2}-\frac{1}{6} f^{\prime \prime \prime}(x) h^{3}+\cdots \tag 2 \end{eqnarray*}$$

两式相减得到

$$\begin{aligned} f(x+h)-f(x-h) &=2 f^{\prime}(x) h+\frac{1}{3} f^{\prime \prime \prime}(x) h^{3}+\cdots \\ \frac{f(x+h)-f(x-h)}{2h} &=f^{\prime}(x)+O\left(h^{2}\right) \end{aligned}$$

即

$$\frac{f(x+h)-f(x-h)}{2h}\approx f^{\prime}(x)$$

该公式被称为中点差分公式，其精确度为$O(h^2)$。

(1),(2)相加得到二阶导数的中点差分公式：

$$\begin{aligned} f(x+h)+f(x-h) &=2 f(x)+f^{\prime \prime}(x) h^{2}+O\left(h^{3}\right) \\ \frac{f(x+h)-2 f(x)+f(x-h)}{h^{2}} &=f^{\prime \prime}(x)+O\left(h^{2}\right) \end{aligned}$$

最后介绍一个提升精度的方法，定义

$$D(h) = \frac{f(x+h)-f(x)}{h}$$

那么$\forall 0 <\alpha <1$，计算$D(h),D(\alpha h)$可得

$$\begin{aligned} D(h)&=f^{\prime}(x)+\frac{1}{2} f^{\prime \prime}(x) h+O\left(h^{2}\right)\\ D(\alpha h)&=f^{\prime}(x)+\frac{1}{2} f^{\prime \prime}(x) \alpha h+O\left(h^{2}\right) \end{aligned}$$

写成矩阵形式得到

$$\begin{array}{c} {\left( \begin{array}{cc}{1} & {\frac{1}{2} h} \\ {1} & {\frac{1}{2} \alpha h}\end{array}\right) \left( \begin{array}{c}{f^{\prime}(x)} \\ {f^{\prime \prime}(x)}\end{array}\right)=\left( \begin{array}{c}{D(h)} \\ {D(\alpha h)}\end{array}\right)+O\left(h^{2}\right)} \end{array}$$

求解可得

$$\begin{aligned} \left( \begin{array}{cc}{f^{\prime}(x)} \\ {f^{\prime \prime}(x)}\end{array}\right) &=\left( \begin{array}{cc}{1} & {\frac{1}{2} h} \\ {1} & {\frac{1}{2} \alpha h}\end{array}\right)^{-1}\left[\left( \begin{array}{c}{D(h)} \\ {D(\alpha h)}\end{array}\right)+O\left(h^{2}\right)\right] \\ &=\frac{1}{1-\alpha} \left( \begin{array}{cc}{-\alpha} & {1} \\ {\frac{2}{h}} & {-\frac{2}{h}}\end{array}\right) \left[\left( \begin{array}{c}{D(h)} \\ {D(\alpha h)}\end{array}\right)+O\left(h^{2}\right)\right] \\ &=\frac{1}{1-\alpha} \left( \begin{array}{cc}{-\alpha} & {1} \\ {\frac{2}{h}} & {-\frac{2}{h}}\end{array}\right) \left( \begin{array}{c}{D(h)} \\ {D(\alpha h)}\end{array}\right)+\left( \begin{array}{c}{O\left(h^{2}\right)} \\ {O(h)}\end{array}\right) \end{aligned}$$

因此

$$f'(x) =\frac 1 {1-\alpha} \left(-\alpha D(h)+D(\alpha h)\right) +O(h^2)$$

这说明我们利用$D(h)$和$D(\alpha h)$得到了对$f’(x)$的计算公式，其精度为$O(h^2)$。