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\) ，即证.
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}
\]
从而证明了余项估计.
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\) ，即证.
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(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}}
\]
这只是一个猜测，暂时不进行验证.
浙公网安备 33010602011771号