数值积分和导数

目录

• Accuracy and convergence
• Newton-Cotes formulas
• Composite formulas
• Gauss formulas
• Numerical differentiation
• Problem

 

Definition 6.1. 加权求积公式 weighted quadrature formula $I_n(f)$ 是线性泛函

$$I_n(f) = \sum_{k=1}^nw_kf(x_k) \tag{6.1}$$

它逼近函数 $f\in\mathcal{C}[a,b]$ 的积分

$$I(f) = \int_a^bf(x)\rho(x)dx \tag{6.2}$$

其中权函数 $\rho\in\mathcal{C}[a,b]$ 满足 $\forall x\in[a,b],\ \rho(x)>0$ ，点 $x_k$ 称为节点 node 或横坐标 abscissas ，乘数 $\omega_k$ 称为权 weight 或系数.

 

Accuracy and convergence

Definition 6.3. $I_n(f)$ 的余项 remainder 或误差 error 为

$$E_n(f) = I(f) - I_n(f) \tag{6.3}$$

称 $I_n(f)$ 在 $\mathcal{C}[a,b]$ 收敛当且仅当

$$\forall f\in\mathcal{C}[a,b],\quad \lim_{n\to\infty} I_n(f) = I(f) \tag{6.4}$$

 
Definition 6.4. 子集 $\mathbb{V}\subset\mathcal{C}[a,b]$ 在 $\mathcal{C}[a,b]$ 上稠密 dense 当且仅当

$$\forall f\in\mathcal{C}[a,b],\ \forall\epsilon>0,\ \exists f_{\epsilon}\in\mathbb{V},\quad \mathrm{s.t.}\quad \max_{x\in[a,b]}|f(x)-f_{\epsilon}(x)| \le \epsilon \tag{6.5}$$

也就是说，可以从 $\mathbb{V}$ 无限逼近 $f\in\mathcal{C}[a,b]$ .

 

Theorem 6.5. 令 $\{I_n(f):n\in\mathbb{N}^+\}$ 为逼近 $I(f)$ 的一列求积公式，令子集 $\mathbb{V}$ 为 $\mathcal{C}[a,b]$ 的稠密子集，则 $I_n(f)$ 在 $\mathcal{C}[a,b]$ 收敛当且仅当

$$\forall f\in\mathbb{V},\quad \lim_{n\to\infty} I_n(f) = I(f)\\ \exist B\in\mathbb{R},\quad \mathrm{s.t.}\quad \forall n\in\mathbb{N}^+,\ W_n = \sum_{k=1}^n|\omega_k| < B\tag{6.6}$$

$Proof.$ 注意第一个条件是在稠密子集中收敛.

 

Definition 6.6. 加权求积公式有精确次数 $d_E$ 当且仅当

$$\left\{ \begin{aligned} \forall f\in\mathbb{P}_{d_E},\quad E_n(f) = 0\\ \exists g\in\mathbb{P}_{d_E+1},\quad \mathrm{s.t.}\quad E_n(g) \neq 0 \end{aligned} \right. \tag{6.7}$$

 
Lemma 6.8. 令 $x_1,\cdots,x_n$ 为 $I_n(f)$ 的互异节点，若希望有 $d_E\ge n-1$ ，则权可以由

$$\forall k=1,\cdots,n,\quad \omega_k = \int_a^b\rho(x)l_k(x)dx \tag{6.8}$$

推导，其中 $l_k(x)$ 为逐点插值的初等多项式

$$l_k(x) = \prod_{i\neq k;i=1}^n\dfrac{x-x_i}{x_k-x_i} \tag{6.9}$$

$Proof.$ 我们考虑在这 $n$ 个点的插值多项式 $p(f;x)$ ，从而有

$$\sum_{k=1}^{n}\omega_kf(x_k) = \sum_{k=1}^{n}f(x_k)\int_a^b\rho(x)l_k(x)dx = \int_a^b\rho(x)\sum_{k=1}^{n}l_k(x)f(x_k)dx = \int_a^b\rho(x)p(f;x)dx$$

当 $f\in\mathbb{P}_{n-1}$ 时，有 $f(x) = p(f;x)$ ，从而 $d_E\ge n-1$ ，即证.

 

Newton-Cotes formulas

Definition 6.9. 牛顿-科茨公式 Newton-Cotes formulas 是通过在均匀节点 $x_1,\cdots,x_n\in[a,b]$ 插值逼近 $I(f)$ 的公式，形如 $(6.1)$ .

 

Definition 6.10. 梯形法 trapezoidal rule 是通过直接连接点 $(a,f(a))^T,\ (b,f(b))^T$ 来逼近 $I(f)$ 的公式

$$I^T(f) = \omega_1f(a) + \omega_2f(b) \tag{6.10}$$

其中权 $\omega_1,\ \omega_2$ 满足

$$\omega_1 = \int_a^b\rho(x)l_1(x)dx = \int_a^b\rho(x)\dfrac{x-b}{a-b}dx\\ \omega_2 = \int_a^b\rho(x)l_2(x)dx = \int_a^b\rho(x)\dfrac{x-a}{b-a}dx\\ \tag{6.11}$$

特别的，当 $\rho(x)\equiv 1$ ，有

$$I^T(f) = \dfrac{b-a}{2}[f(a)+f(b)] \tag{6.12}$$

此时逼近结果就是它们围成的梯形面积.

 

Theorem 6.12. 对带有权函数 $\rho(x)\equiv 1$ 的 $f\in\mathcal{C}^2[a,b]$ ，梯形法的余项满足

$$\exists\zeta\in[a,b],\quad \mathrm{s.t.}\quad E^T(f) = -\dfrac{(b-a)^3}{12}f^{\prime\prime}(\zeta) \tag{6.13}$$

$Proof.$ 注意到

$$I^T(f) = \int_a^b[\rho(x)f(a)l_1(x) + \rho(x)f(b)l_2(x)]dx = \int_a^b[f(a)l_1(x) + f(b)l_2(x)]dx$$

也就是说，被积项实际上是两点间的插值多项式，因此取柯西余项，有

$$\begin{aligned} E^T(f) &= -\int_a^b\dfrac{f^{\prime\prime}(\xi(x))}{2}(x-a)(b-x)dx\\ &= -\dfrac{f^{\prime\prime}(\zeta)}{2}\int_a^b(x-a)(b-x)dx\\ &= -\dfrac{(b-a)^3}{12}f^{\prime\prime}(\zeta) \end{aligned}$$

其中使用了积分中值定理，这是因为 $(x-a)(b-x)$ 在 $(a,b)$ 上不变号.

 

Definition 6.13. 辛普森法 Simpson's rule 是通过经过点 $(a,f(a))^T,\ (b,f(b))^T,\ \left(\frac{a+b}{2},\frac{f(a)+f(b)}{2}\right)^T$ 来逼近 $I(f)$ 的公式，特别的，当权函数 $\rho(x)\equiv 1$ 时，有

$$I^S(f) = \dfrac{b-a}{6}\left[f(a)+4f\left(\dfrac{a+b}{2}\right)+f(b)\right] \tag{6.14}$$

 
Theorem 6.14. 对带有权函数 $\rho(x)\equiv 1$ 的 $f\in\mathcal{C}^4[a,b]$ ，辛普森法的余项满足

$$\exist \zeta\in[a,b],\quad \mathrm{s.t.}\quad E^S(f) = -\dfrac{(b-a)^5}{2880}f^{(4)}(\zeta) \tag{6.15}$$

此定理的证明需要一些技巧，因为 $(x-a)(x-b)(x-\frac{a+b}{2})$ 的符号在 $(a,b)$ 上变化.

$Proof.$ 我们将证明分为两个部分

首先考虑在 $[-1,1]$ 上的辛普森公式

$$\int_{-1}^1y(t)dt = \int_{-1}^1p_3(y;-1,0,0,1;t)dt + E^S(t)$$

其中 $y\in\mathcal{C}^4[-1,1]$ ，而 $p_3(y;-1,0,0,1;t)$ 是在满足 $p_3(-1)=y(-1),\ p_3(0)=y(0),\ p_3^{\prime}(0)=y^{\prime}(0),\ p_3(1)=y(1)$ 的插值多项式，计算差分表得到

$$\begin{aligned} p_3(t) &= y(-1)+[y(0)-y(-1)](t+1)\\ &+[y^{\prime}(0)-y(0)-y(-1)]t(t+1)\\ &+\frac{1}{2}\left[y(1)-2y^{\prime}(0)-y(-1)\right]t^{2}(t+1) \end{aligned}$$

从而有积分值为辛普森公式

$$\begin{aligned} \int_{-1}^{1}p_3(t)dt &= 2y(-1)+2[y(0)-y(-1)]+\frac{2}{3}[y^{\prime}(0)-y(0)-y(-1)]\\ &+\frac{1}{3}\left[y(1)-2y^{\prime}(0)-y(-1)\right]\\ &= \frac{1}{3}[y(-1)+4y(0)+y(1)]\\ &= I_S(y) \end{aligned}$$

因此可以通过这种方式得到辛普森公式的余项

$$\begin{aligned} E_S(y) &= -\int_{-1}^{1}\dfrac{f^{4}(\xi(x))}{4!}t^{2}(t+1)(1-t)dt\\ &= -\dfrac{f^{(4)}(\zeta)}{4!}\int_{-1}^{1}t^{2}(t+1)(1-t)dt\\ &= -\dfrac{f^{(4)}(\zeta)}{90} \end{aligned}$$

这里被积项是埃尔米特插值的柯西余项，于是我们证明了在 $[-1,1]$ 上的余项；

考虑变量代换 $g:[-1,1]\to[a,b]$ ，定义

$$g(t) = \dfrac{b-a}{2}t + \dfrac{b+a}{2}$$

我们仍然令 $p_3(t)$ 为 $y(t)$ 的插值多项式，但是插值条件改为 $p_3(a)=y(a),\ p_3(\frac{a+b}{2})=y(\frac{a+b}{2}),\ p_3^{\prime}(\frac{a+b}{2})=y^{\prime}(\frac{a+b}{2}),\ p_3(b)=y(b)$ ，则有

$$\int_{a}^{b}p_3(t)dt = \int_{-1}^{1}p_3(g(t))g^{\prime}(t)dt = \int_{-1}^{1}\dfrac{b-a}{2}p_3\left(\dfrac{b-a}{2}t + \dfrac{b+a}{2}\right)dt$$

对应的，我们改写插值多项式和插值函数，令

$$h(t) = \dfrac{b-a}{2}p_3\left(\dfrac{b-a}{2}t + \dfrac{b+a}{2}\right)\\ z(t) = \dfrac{b-a}{2}y\left(\dfrac{b-a}{2}t + \dfrac{b+a}{2}\right)$$

应用原先的插值条件，我们有

$$\begin{aligned} &h(-1) = \dfrac{b-a}{2}y(a) = z(-1),\\ &h(0) = \dfrac{b-a}{2}y\left(\dfrac{b+a}{2}\right) = z(0),\\ &h^{\prime}(0) = \left(\dfrac{b-a}{2}\right)^{2}y^{\prime}\left(\dfrac{b+a}{2}\right) = z^{\prime}(0),\\ &h(1) = \dfrac{b-a}{2}y(b) = z(1) \end{aligned}$$

也就是说，我们回到了 $[-1,1]$ 上的插值，应用先前的结论，有

$$\int_{a}^{b}p_3(t)dt = \int_{-1}^{1}h(t)dt = \dfrac{1}{3}[z(-1)+4z(0)+z(1)] = \dfrac{b-a}{6}\left[y(a)+4y\left(\dfrac{b+a}{2}\right)+y(b)\right]$$

这恰好就是辛普森公式，因此我们可以推导埃尔米特插值的柯西余项

$$\begin{aligned} E_S(f) &= -\int_{a}^{b}\dfrac{f^{4}(\xi(x))}{4!}\left(t-\dfrac{b+a}{2}\right)^{2}(t-a)(b-t)dt\\ &= -\dfrac{f^{(4)}(\zeta)}{4!}\int_{a}^{b}\left(t-\dfrac{b+a}{2}\right)^{2}(t-a)(b-t)dt \end{aligned}$$

可以发现，通过在中点进行两次插值，我们得到了在 $[a,b]$ 上不变号的多项式，作代换 $d = \frac{b-a}{2}$ 有

$$\begin{aligned} E_S(f) &= -\dfrac{f^{(4)}(\zeta)}{4!}\int_{-d}^{d}t^{2}(t+d)(d-t)dt\\ &= -\dfrac{f^{(4)}(\zeta)}{90}d^{5}\\ &= -\dfrac{(b-a)^{5}}{2880}f^{(4)}(\zeta) \end{aligned}$$

从而证明了余项估计.

 

Composite formulas

Definition 6.16. 带有权函数 $\rho\equiv 1$ 逼近 $f(x)$ 的复合梯形法 composite trapezoidal rule 公式为

$$I_n^T(f) = \dfrac{h}{2}f(x_0) + h\sum_{k=1}^{n-1}f(x_k) + \dfrac{h}{2}f(x_n) \tag{6.16}$$

其中 $h=\frac{b-a}{n},\ x_k=a+kh$ ；它是梯形法在每一段长度为 $h$ 的 $[x_i,x_{i+1}]$ 上的简单相加.

 

Theorem 6.17. 对 $f\in\mathcal{C}^2[a,b]$ ，复合梯形法的余项满足

$$\exists\xi\in(a,b),\quad \mathrm{s.t.}\quad E_n^T(f) = -\dfrac{b-a}{12}h^2f^{\prime\prime}(\xi) \tag{6.17}$$

$Proof.$ 在每一段区间上的余项相加即证，注意由于每段的 $\xi_k$ 不同，需要应用连续函数的介值定理得到一个共同的 $\xi$ .

 

Definition 6.18. 带有权函数 $\rho\equiv 1$ 逼近 $I(f)$ 的复合辛普森法 composite Simpson's rule 公式为

$$I_n^S(f) = \dfrac{h}{3}\sum_{k=1}^{n/2}[f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})] \tag{6.18}$$

其中 $h=\frac{b-a}{n},\ x_k=a+kh$ ，$n$ 为偶数；它是辛普森法在每一段长度为 $2h$ 的区间 $[x_i,x_{i+2}]$ 上的简单相加.

 

Theorem 6.19. 对 $f\in\mathcal{C}^4[a,b]$ 及 $n\in2\mathbb{N}^+$ ，复合辛普森法的余项满足

$$\exists\xi\in(a,b),\quad \mathrm{s.t.}\quad E_n^S(f) = -\dfrac{b-a}{180}h^4f^{(4)}(\xi) \tag{6.19}$$

$Proof.$ 在每一段区间上的余项相加即证.

 

Lemma 6.20. 梯形法满足

$$\forall f\in\mathbb{T}_{n-1}[0,2\pi],\quad E_n^T(f) = 0 \tag{6.20}$$

其中 $\mathbb{T}_{n}[0,2\pi]$ 是至多 $n$ 次的三角函数的多项式类

$$\mathbb{T}_{n}[0,2\pi] = \mathrm{span}\{1,\cos x,\sin x,\cdots,\cos nx,\sin nx,\} \tag{6.21}$$

$Proof.$ 我们只需要考虑 $e_m(x) = e^{imx} = \cos mx + i\sin mx,\ m\in\mathbb{N}$ ，则

$$\begin{aligned} E_n^T(e_m) &= \int_0^{2\pi} e_m(x)dx - \dfrac{2\pi}{n}\left[\dfrac{e_m(0)+e_m(2\pi)}{2}+\sum_{k=1}^{n-1}e_m\left(\dfrac{2k\pi}{n}\right)\right]\\ &= \int_0^{2\pi} e^{imx}dx - \dfrac{2\pi}{n}\sum_{k=0}^{n-1}e^{imk\cdot 2\pi/n} \end{aligned}$$

容易证明 $E_n^T(e_m) = 0,\ m = 0,\cdots,n-1$ ，即证.

 

Gauss formulas

Lemma 6.21. 令 $n,m\in\mathbb{N}^+,\ m\le n$ ，给定多项式 $p=\sum_{i=0}^{n+m}p_ix^i\in\mathbb{P}_{n+m}$ 以及 $s=\sum_{i=0}^ns_ix^i\in\mathbb{P}_n$ 满足 $p_{n+m}\neq 0,\ s_n\neq 0$ ，则存在唯一多项式 $q\in\mathbb{P}_m,\ r\in\mathbb{P}_{n-1}$ 使得

$$p = qs + r \tag{6.22}$$

$Proof.$ 此引理本质上是带余除法得到的，容易证明.

 

Definition 6.22. 在节点 $x_k$ 上的节点多项式 node polynomial 为

$$v_n(x) = \prod_{k=1}^n(x-x_k) \tag{6.23}$$

 
Theorem 6.23. 设求积公式有 $d_E\ge n-1$ ，则能且只能通过添加关于它的节点多项式条件

$$\forall p\in\mathbb{P}_{j-1},\quad \int_a^bv_n(x)p(x)\rho(x)dx = 0 \tag{6.24}$$

改善精确次数为 $d_E\ge n+j-1,\ j\in(0,n]$ .

$Proof.$ 若有精确次数 $d_E\ge n+j-1,\ j\in(0,n]$ ，由 $v_n(x)p(x)\in\mathbb{P}_{n+j-1}$ 显然有

$$\int_a^bv_n(x)p(x)\rho(x)dx = \sum_{k=1}^{n}\omega_kv_n(x_k)p(x_k) = 0$$

必要性得证；对于充分性，我们取 $p(x)\in\mathbb{P}_{n+j-1}$ ，则可分解为 $p = qv_n + r$ ，有

$$\int_a^bp(x)\rho(x)dx = \int_a^bv_n(x)q(x)\rho(x)dx + \int_a^br(x)\rho(x)dx = \int_a^br(x)\rho(x)dx$$

由于余项 $r\in\mathbb{P}_{n-1}$ ，则有

$$\int_a^br(x)\rho(x)dx = \sum_{k=1}^n\omega_kr(x_k) = \sum_{k=1}^n\omega_k[p(x_k)-q(x_k)v_n(x_k)] = \sum_{k=1}^n\omega_kp(x_k)$$

即证.

 

Note. 注意到对于节点多项式 $v_n(x)$ 的所有零点的插值多项式恰为 $n-1$ 次，这也就意味着，对任意 $p\in\mathbb{P}_{n-1}$ ，以这些零点为节点推导的加权求积公式都是完全精确的，由 Lemma 6.8 可知，天然地有 $d_E \ge n-1$ ，因此这一条件其实是自然的，但是并不显然.

 

Definition 6.24. 高斯求积公式 Gaussian quadrature formula 以节点多项式 $v_n(x)$ 的零点为节点，在条件 $(6.24)$ 中 $j=n$ .

 

Corollary 6.25. 高斯公式有 $d_E = 2n-1$ .

$Proof.$ 在 Theorem 6.23 中，我们可将节点条件看做

$$\forall p\in\mathbb{P}_{n-1},\quad \left<v_n,p\right> = 0$$

即 $v_n$ 与所有不大于 $n-1$ 次多项式正交，它们共有 $n$ 个方向，那么显然最多改善 $n$ ，因为 $v_n\in\mathbb{P}_n$ ，它不可能与自己正交；总的来看，我们可以这样理解：在加权求积公式中， $\omega_k,\ x_k$ 为 $2n$ 个变量，因此可以通过 $2n$ 阶线性系统确定至多 $2n-1$ 次多项式.

 

Corollary 6.26. 高斯公式 $I_n(f)$ 的权为

$$\forall k=1,\cdots,n,\quad \omega_k=\int_a^b\dfrac{v_n(x)}{(x-x_k)v_n^{\prime}(x_k)}\rho(x)dx \tag{6.25}$$

$Proof.$ 这是由 Lemma 6.8 直接得到的，是在 $v_n(x)$ 节点处的 $l_k(x)$ .

 

Note. 推导高斯公式

对于给定的 $n$ 和权函数 $\rho(x)$ ，要求节点多项式 $v_n(x)$ 来确定插值点，注意它是首一多项式；

我们首先根据 Corollary 6.25 要求 $\pi_n(x)$ 应当与 $x^i,\ i=0,1,\cdots,n-1$ 正交，即有

$$\int_a^b\rho(x)\pi_n(x)x^idx = \int_a^b\rho(x)(c_0+c_1x+\cdots+x^n)x^idx = 0$$

这样得到关于系数 $c_i,\ i=0,1,\cdots,n-1$ 的 $n$ 个线性方程，从而可以解出 $\pi_n(x)$ ；

通过 $\pi_n(x)$ 的零点确定 $n$ 个插值点，然后求取对应各点的权 $\omega_i$ ，当然可以通过 Corollary 6.26 来计算权，但是这样每次需要计算一个复杂的高次积分，因此我们考虑利用插值多项式的特性。由于在 $\pi_n(x)$ 零点的求积公式对于不高于 $n-1$ 次的多项式是精确的，则

$$\begin{aligned} \omega_1+\omega_2+\cdots+\omega_n &= \int_a^b\rho(x)x^0dx\\ x_1\omega_1+x_2\omega_2+\cdots+x_n\omega_n &= \int_a^b\rho(x)x^1dx\\ &\vdots\\ x_1^{n-1}\omega_1+x_2^{n-1}\omega_2+\cdots+x_n^{n-1}\omega_n &= \int_a^b\rho(x)x^{n-1}dx \end{aligned}$$

其中 $x_1,x_2,\cdots,x_n$ 是 $\pi_n(x)$ 的零点，我们利用对单项式的求积公式得到 $n$ 个线性方程，从而可以解出 $\omega_i,\ i=1,2,\cdots,n$ ；

最终有高斯公式

$$I_n^G(f) = \omega_1f(x_1) + \omega_2f(x_2) + \cdots + \omega_nf(x_n)$$

它对于不高于 $2n-1$ 次的多项式是精确的.

 

Definition 6.28. 正交多项式集 orthogonal polynomials 是满足

$$\forall p_i,p_j\in P,\quad i\neq j \Rightarrow \left<p_i,p_j\right> = 0 \tag{6.26}$$

的集合 $P = \{p_i:\deg(p_i) = i\}$ .

 

Theorem 6.30. 在 $[a,b]$ 上实正交多项式的每一个零点都是实值单零点，并且在 $(a,b)$ 中.

$Proof.$ 注意到 $p_0$ 是常多项式，从 $n\ge 1$ 开始，有

$$\int_a^b\rho(x)p_n(x)p_0(x)dx = \left<p_n,p_0\right> = 0$$

若 $p_n(x)$ 在 $[a,b]$ 上没有零点，则它不变号，因此

$$\int_a^b\rho(x)p_n(x)p_0(x)dx = p_0\int_a^b\rho(x)p_n(x)dx \neq 0$$

矛盾，故至少有根 $x_1\in[a,b],\ p_n(x_1)=0$ ；若有重根，则 $\frac{p_n(x)}{(x-x_1)^2}$ 是一个 $n-2$ 次多项式，从而

$$0 = \left<p_n(x),\dfrac{p_n(x)}{(x-x_1)^2}\right> = \left<1,\dfrac{p_n^2(x)}{(x-x_1)^2}\right> > 0$$

矛盾，则 $x_1$ 只能是单根；

设只有 $j<n$ 个零点 $x_1,\cdots,x_j$ 在 $(a,b)$ 中，令 $v_j = \prod_{i=1}^j(x-x_i)\in\mathbb{P}_{j}$ ，显然 $v_j$ 是 $p_n$ 的因式，则 $p_nv_j = P_{n-1}v_j^2$ ，其中 $P_{n-j}$ 是在 $[a,b]$ 上不变号的 $n-j$ 次多项式，则有

$$0 = \left<p_n,v_j\right> = \left<P_{n-1},v_j^2\right> \neq 0$$

矛盾，因此所有零点在 $(a,b)$ 中.

 

Corollary 6.31. 高斯公式的所有节点都互异实值，并且在 $(a,b)$ 中.

 

Lemma 6.32. 高斯公式有正权.

$Proof.$ 由于初等插值多项式 $l_k(x)$ 满足 $l_k(x_j) = l_k^2(x_j),\ j=1,\cdots,n$ ，并且 $l_k^2(x)\in\mathbb{P}_{2n-1}$ ，则有

$$\omega_k = \sum_{j=1}^nw_jl_k^2(x_j) = \int_a^b\rho(x)l_k^2(x)dx > 0$$

即证.

 

Lemma 6.33. 高斯公式满足

$$\sum_{k=1}^n\omega_k = \mu_0\in(0,+\infty) \tag{6.27}$$

$Proof.$ 这里 $\mu_0$ 是由权函数的重量

$$\mu_j = \int_a^bx^j\rho(x)dx,\quad j\in\mathbb{N}$$

定义的；只需令 $j=0$ ，即有对 $I(f),\ f(x)\equiv 1$ 进行逼近，由定义即证.

 

Theorem 6.34. 高斯公式在 $\mathcal{C}[a,b]$ 上收敛.

$Proof.$ 不妨设 $\mathbb{P}$ 为实多项式集，则由魏尔斯特拉斯逼近定理， $\mathbb{P}$ 在 $\mathcal{C}[a,b]$ 上稠密，并且由于

$$\sum_{k=1}^n|\omega_k| = \sum_{k=1}^n\omega_k = \mu_0 < +\infty$$

由 Theorem 6.5 即证.

 

 
Theorem 6.35. 对 $f\in\mathcal{C}^{2n}[a,b]$ ，高斯公式 $I_n(f)$ 的余项满足

$$\exists \xi\in[a,b],\quad \mathrm{s.t.}\quad E_n^G(f) = \dfrac{f^{(2n)}(\xi)}{(2n)!}\int_a^b\rho(x)v_n^2(x)dx \tag{6.28}$$

其中 $v_n$ 是确定 $I_n$ 的节点多项式.

 

Numerical differentiation

Formula 6.36 (The method of undetermined coefficients). 推导 FD 公式来逼近 $u^{(k)}(\overline{x})$ 的一般方法基于任意互异点 $x_1,x_2,\cdots,x_n,\ n>k$ 构成的模板 stencil ；在模板上逐点 $x_i$ 泰勒展开得到

$$u(x_i) = u(\overline{x}) + (x_i-\overline{x})u^{\prime}(\overline{x}) + \cdots + \dfrac{1}{k!}(x_i-\overline{x})^ku^{(k)}(\overline{x}) + \cdots \tag{6.29}$$

这样可以把 $u^{(k)}(\overline{x})$ 表示为 $u(x_i)$ 的线性组合

$$u^{(k)}(\overline{x}) = c_1u(x_1) + c_2u(x_2) + \cdots + c_nu(x_n) + O(h^p) \tag{6.30}$$

其中 $c_i$ 是为使 $p$ 足够大而特别选取的，有

$$\forall i=0,\cdots,p-1,\quad \dfrac{1}{i!}\sum_{j=1}^nc_j(x_j-\overline{x})^i = \left\{ \begin{aligned} &1,\quad i=k\\ &0,\quad \mathrm{otherwise} \end{aligned} \right. \tag{6.31}$$

即将右式泰勒展开后，只有 $k$ 阶导数的项的系数为 $1$ ，其它项系数都为 $0$ .

 

Definition 6.38. 逼近导数的 FD 公式称为是 $p$ 阶精确，如果误差 $E$ 有形式

$$E(h) = \Theta(h^p) \tag{6.32}$$

其中 $h$ 为模板中相邻点的最大距离.

 

Lemma 6.42. 逼近 $u\in\mathcal{C}^4(\mathbb{R})$ 的二阶导数公式

$$D^2u(\overline{x}) = \dfrac{u(\overline{x}-h)-2u(\overline{x})+u(\overline{x}+h)}{h^2} \tag{6.33}$$

是二阶精确，进一步若输入函数 $u(\overline{x}-h),\ u(\overline{x}),\ u(\overline{x}+h)$ 有误差范围 $\epsilon\in[-E,E]$ ，则存在 $\xi\in[\overline{x}-h,\overline{x}+h]$ 使得

$$\left|u^{\prime\prime}(\overline{x})-D^2u(\overline{x})\right| \le \dfrac{h^2}{12}\left|u^{(4)}(\xi)\right| + \dfrac{4E}{h^2} \tag{6.34}$$

$Proof.$ 构造 FD 公式

$$D^{2}u(\overline{x}) = au(\overline{x}-h) + bu(\overline{x}) + cu(\overline{x}+h)$$

在 $x = \overline{x}$ 处 Taylor 展开得

$$\begin{aligned} D^{2}u(\overline{x}) &= (a + b + c)u(\overline{x}) + (-a+c)hu^{\prime}(\overline{x}) \\ &+ \frac{1}{2}(a+c)h^{2}u^{\prime\prime}(\overline{x}) + \frac{1}{6}(-a+c)h^{3}u^{\prime\prime\prime}(\overline{x}) + O(h^{4}) \end{aligned}$$

求解线性方程组

$$\left\{ \begin{aligned} &a+b+c = 0\\ &-a+c = 0\\ &a+c = \frac{2}{h^{2}} \end{aligned} \right.$$

得到 $a = c = \frac{1}{h^{2}},\ b = -\frac{2}{h^{2}}$ ，从而有

$$D^{2}u(\overline{x}) = \dfrac{u(\overline{x}-h) - 2u(\overline{x}) + u(\overline{x}+h)}{h^{2}}$$

注意到 $-a + c = 0$ ，则 $3$ 阶项为 $0$ ，而由 $4$ 阶展开的 Lagrange 余项

$$\dfrac{1}{24}(a+c)h^{4}f^{(4)}(\xi) = \dfrac{h^{2}}{12}f^{(4)}(\xi),\ \xi\in[\overline{x}-h,\overline{x}+h]$$

再由每一项的扰动 $\epsilon\in[-E,E]$ 有

$$\left|\Delta D^{2}u(\overline{x})\right|\le \dfrac{\left|\Delta u(\overline{x}-h)\right| + 2\left|\Delta u(\overline{x})\right| + \left|\Delta u(\overline{x}+h)\right|}{h^{2}} = \dfrac{4E}{h^{2}}$$

从而有

$$\left|u^{\prime\prime}(\overline{x})-D^{2}u(\overline{x})\right|\le\dfrac{h^{2}}{12}\left|f^{(4)}(\xi)\right| + \dfrac{4E}{h^{2}}$$

应用均值不等式

$$\dfrac{h^{2}}{12}\left|f^{(4)}(\xi)\right| + \dfrac{4E}{h^{2}}\ge \sqrt{\dfrac{4}{3}E\left|f^{(4)}(\xi)\right|}$$

其中等式成立当且仅当

$$h = \sqrt[4]{\dfrac{48E}{\left|f^{(4)}(\xi)\right|}}$$

此时误差界最小，即证.

 

Problem

$\mathrm{IV}.$ 高斯公式的余项。考虑埃尔米特插值问题：寻找 $p\in\mathbb{P}_{n-1}$ 使得

$$\forall m=1,2,\cdots,n,\quad p(x_m) = f_m,\ p^{\prime}(x_m) = f^{\prime}_m$$

有初等埃尔米特插值多项式 $h_m,q_m$ 使得问题的解可以表示为

$$p(t) = \sum_{m=1}^n[h_m(t)f_m+q_m(t)f_m^{\prime}]$$

类似于拉格朗日插值公式.

• 寻找 $h_m,q_m$ 具有形式

$$h_m(t) = (a_m+b_mt)l_m^2(t),\quad q_m(t) = (c_m+d_mt)l_m^2(t)$$

其中 $l_m$ 是初等拉格朗日多项式.

首先，根据这种形式，我们希望 $h_m(t)$ 只在 $x_m$ 非零且为 $1$ ，并且导数在节点恒为 $0$ ；令 $h_m(t)$ 在 $x_1,\cdots, x_n$ 处满足

$$h_m(x_i) = \left\{ \begin{aligned} &1,\quad i = m\\ &0,\quad i \neq m \end{aligned} \right.\quad\quad h_m^{\prime}(x_i) = 0,\ i = 1,\cdots,n$$

观察 $l_m$ 的形式

$$l_m(t) = \prod_{i = 1\atop i \neq m}^{n}\dfrac{t-x_i}{x_m-x_i}\quad\quad l_m^{\prime}(t) \sum_{i = 1\atop i \neq m}^{n}\left(\dfrac{1}{x_m-x_i}\prod_{j = 1\atop j \neq m,i}^{n}\dfrac{t-x_j}{x_m-x_j}\right)$$

显然有

$$l_m(x_i) = \left\{ \begin{aligned} &1,\quad i = m\\ &0,\quad i \neq m \end{aligned} \right.\quad\quad l_m^{\prime}(x_m) = \sum_{i = 1\atop i \neq m}^{n}\dfrac{1}{x_m-x_i}$$

根据 $h_m(t)$ 的形式有

$$h_m(t) = (a_m+b_mt)l_m^{2}(t),\quad h_m^{\prime}(t) = b_ml_m^{2}(t) + 2(a_m+b_mt)l_m^{\prime}(t)l_m(t)$$

代入节点，按照我们的要求，只需求解线性方程组

$$\left\{ \begin{aligned} &a_m+b_mx_m = 1\\ &b_m+2(a_m+b_mx_m)l_m^{\prime}(x_m) = 0 \end{aligned} \right.$$

得到

$$a_m = 1 + \sum_{i = 1\atop i \neq m}^{n}\dfrac{2x_m}{x_m-x_i},\quad b_m = -\sum_{i = 1\atop i \neq m}^{n}\dfrac{2}{x_m-x_i}$$

类似的，要求

$$q_m(x_i) = 0,\ i = 1,\cdots,n\quad\quad q_m^{\prime}(x_i) = \left\{ \begin{aligned} &1,\quad i = m\\ &0,\quad i \neq m \end{aligned} \right.$$

得到线性方程组

$$\left\{ \begin{aligned} &c_m+d_mx_m = 0\\ &d_m+2(c_m+d_mx_m)l_m^{\prime}(x_m) = 1 \end{aligned} \right.$$

解得

$$c_m = -x_m,\quad d_m = 1$$

代回我们便得到需要的多项式.

 

• 推导求积公式

$$I_n(f) = \sum_{k=1}^n[\omega_kf(x_k)+\mu_kf^{\prime}(x_k)]$$

满足 $E_n(p)=0,\ \forall p\in\mathbb{P}_{2n-1}$ .

直接对先前所得 $p(x)$ 积分得到

$$I_n(f) = \sum_{k=1}^{n}\left[\omega_kf(x_k)+\mu_kf^{\prime}(x_k)\right]$$

其中

$$\omega_k = \int_{a}^{b}\rho(t)h_m(t)dt,\quad \mu_k = \int_{a}^{b}\rho(t)q_m(t)dt$$

由于 $p(t)$ 实际上是 $f(x)$ 的 Hermite 插值多项式，根据 Hermite 插值的余项，对于 $f\in\mathbb{P}_{2n-1}$ ，有 $f(t) = p(t)$ ，从而 $E_n(f) = 0$ .

 

• 在什么条件下 $\mu_k=0,\ k=1,2,\cdots,n$ .

根据要求有

$$\int_{a}^{b}\rho(t)q_k(t)dt = 0,\ k=1,\cdots,n$$

节点多项式为 $v_n(x) = \prod_{k=1}^{n}(x-x_k)$ ，则

$$q_k(t) = (t-x_k)\left(\prod_{i = 1\atop i \neq k}^{n}\dfrac{t-x_i}{x_k-x_i}\right)^{2} = v_n(t)l_k(t)\prod_{i = 1\atop i \neq k}^{n}\dfrac{1}{x_k-x_i}$$

代入得

$$\int_{a}^{b}\rho(t)v_n(t)l_k(t)dt = 0,\ k=1,\cdots,n$$

则节点多项式满足 $\left<v_n, l_k\right> = 0,\ k=1,\cdots,n$ ，即节点多项式与任意初等 Lagrange 多项式正交.

 

$\mathrm{V}.$ 在对称模板上构造 $4$ 阶精确的公式.

设有 FD 公式

$$D^{4}u(\overline{x}) = au(\overline{x}-2h) + bu(\overline{x}-h) + cu(\overline{x}) + du(\overline{x}+h) + eu(\overline{x}+2h)$$

在 $x = \overline{x}$ 处 Taylor 展开得

$$\begin{aligned} D^{4}u(\overline{x}) &= (a+b+c+d+e)u(\overline{x}) + (-2a-b+d+2e)hu^{\prime}(\overline{x})\\ &+ \dfrac{1}{2}(4a+b+d+4e)h^{2}u^{\prime\prime}(\overline{x}) + \dfrac{1}{6}(-8a-b+d+8e)h^{3}u^{\prime\prime\prime}(\overline{x})\\ &+ \dfrac{1}{24}(16a+b+d+16e)h^{4}u^{(4)}(\overline{x}) + O(h^{5}) \end{aligned}$$

从而有线性方程组

$$\left\{ \begin{aligned} &a+b+c+d+e = 0\\ &-2a-b+d+2e = 0\\ &4a+b+d+4e = \frac{2}{h^{2}}\\ &-8a-b+d+8e = 0\\ &16a+b+d+16e = 0 \end{aligned} \right.$$

解得

$$a = -\dfrac{1}{12h^{2}},\ b = \dfrac{4}{3h^{2}},\ c = -\dfrac{5}{2h^{2}},\ d = \dfrac{4}{3h^{2}},\ e = -\dfrac{1}{12h^{2}}$$

从而有

$$D^{4}u(\overline{x}) = \dfrac{-u(\overline{x}-2h)+16u(\overline{x}-h)-30u(\overline{x})+16u(\overline{x}+h)-u(\overline{x}+2h)}{12h^{2}}$$

注意到 $5$ 阶余项为 $0$ ，而 $6$ 阶余项为

$$\dfrac{h^{6}}{6!}(64a+b+d+64e)f^{(6)}(\xi) = -\dfrac{h^{4}}{90}f^{(6)}(\xi),\ \xi\in[\overline{x}-2h,\overline{x}+2h]$$

因此它是 $4$ 阶精确的，并且

$$\left|u^{\prime\prime}(\overline{x})-D^{4}u(\overline{x})\right|\le\dfrac{h^{4}}{90}\left|f^{(6)}(\xi)\right| + \dfrac{16E}{3h^{2}}$$

由均值不等式有

$$\dfrac{h^{4}}{90}\left|f^{(6)}(\xi)\right| + \dfrac{8E}{3h^{2}} + \dfrac{8E}{3h^{2}}\ge 4\sqrt[3]{\dfrac{E^{2}\left|f^{(6)}(\xi)\right|}{30}}$$

其中等式成立当且仅当

$$h = \sqrt[6]{\dfrac{240E}{\left|f^{(6)}(\xi)\right|}}$$

此时误差界最小.

通过观察可以发现，对于 $u^{\prime\prime}$ 的逼近公式具有类似的形式

$$\begin{aligned} D^2u(\overline{x}) &= \dfrac{u(\overline{x}-h)-2u(\overline{x})+u(\overline{x}+h)}{h^2}\\ D^{4}u(\overline{x}) &= \dfrac{-u(\overline{x}-2h)+16u(\overline{x}-h)-30u(\overline{x})+16u(\overline{x}+h)-u(\overline{x}+2h)}{12h^{2}} \end{aligned}$$

进行简单的作差，我们有

$$\begin{aligned} D^2u(\overline{x}) &= 2\cdot\dfrac{u(\overline{x}-h)-2u(\overline{x})+u(\overline{x}+h)}{2!h^2}\\ D^{4}u(\overline{x}) - D^2u(\overline{x}) &= 2\cdot\dfrac{-u(\overline{x}-2h)+4u(\overline{x}-h)-6u(\overline{x})+4u(\overline{x}+h)-u(\overline{x}+2h)}{4!h^{2}} \end{aligned}$$

其中 $2$ 阶精确公式系数 $1\ 2\ 1$ ， $4$ 阶精确公式增加系数 $1\ 4\ 6\ 4\ 1$ ，我们可以合理猜测 $6$ 阶精确公式

$$D^{6}u(\overline{x}) - D^{4}u(\overline{x}) = 2\cdot\dfrac{u(\overline{x}-3h)-6u(\overline{x}-2h)+15u(\overline{x}-h)-20u(\overline{x})+15u(\overline{x}+h)-6u(\overline{x}+2h)+u(\overline{x}+3h)}{6!h^{2}}$$

这只是一个猜测，暂时不进行验证.