微积分数值计算(02)：外推积分法, Gauss 积分法
外推积分法
所谓 外推积分法, 就是将积分区间大小为\(h\)的数值积分值, 按照一定的形式外推至\(h\to0\)的情况. 比如我们如果考虑积分值满足如下表达式 \[ I(h)=\tau_0+\tau_1h^2+\cdots+\tau_mh^{2m}+\cdots \] 那么我们就可以通过递推计算获得关于\(h\)的某一阶的数值积分.
外推积分法的基础是 Euler-Maclaurin 公式 \[ T(h)-\int_a^bf(t)dt=\sum_{k=1}^m\frac{B_{2k}h^{2k}}{(2k)!}[f^{(2k-1)}(a)-f^{(2k-1)}(b)]-\frac{B_{2m+2}h^{2m+2}}{(2m+2)!}f^{(2m+2)}(\xi)\ ,\ \xi\in[a,b] \] 其中, \(T(h)\)和\(\displaystyle\int_a^bf(t)dt\)分别是微元法给出的积分值和真实的积分值. \(B_{2k}\)是所谓的 Bernoulli 数, 即 \[ \frac{z}{e^z-1}=\sum_{n=0}^\infty\frac{B_n}{n!}z^n \] 根据这个公式, \(T(h)\)可以写为形如 \[ T(h)=\tau_0+\tau_1h^2+\cdots+\tau_mh^{2m}+o(h^{2m}) \] 的多项式. 因此取多个\(h\)并计算\(T(h)\)的值, 由此进行多项式插值, 得到的插值函数的零阶项, 便是对积分的数值计算.
Neville 形式
我们可以模拟 Neville 插值法进行外推. 首先, 选取一个整数序列\(\{N_i\}\)来调整\(h\), 即\(\{h_i\}=\{\displaystyle\frac{b-a}{N_i}\}\). 为了方便, 将\(T(h_k)\)记作\(T_{k0}\). 而内插的\(m\)阶多项式记为\(P_{mm}(h)\), 它显然满足 \[ P_{mm}(h_k)=T_{k0} \] 假设\(P_{jk}(h)\)是关于\(h^2\)的\(k\)次插值函数, 支撑点为\(\{h_{j-k},\cdots,h_{j}\}\). 模仿 Neville 迭代有 \[ P_{jk}=\frac{\displaystyle\frac{h_{j-k}^2}{h_j^2}P_{j,k-1}-P_{j-1,k-1}}{\displaystyle\frac{h_{j-k}^2}{h_j^2}-1} \] 我们可以这么表示迭代过程: \[ \begin{pmatrix} T_{00} && T_{10} && T_{20} && \cdots && T_{m0} \newline & \searrow & \downarrow & \searrow & \downarrow &\searrow & \cdots & \searrow & \downarrow \newline &&P_{11} && P_{21} && \cdots && P_{m1} \newline &&& \searrow & \downarrow &\searrow & \cdots & \searrow & \downarrow \newline &&&& P_{22} && \cdots && P_{m2} \newline &&&&&\searrow & \cdots & \searrow & \downarrow \newline &&&&&& \cdots && \cdots \newline &&&&&&& \searrow & \downarrow \newline &&&&&&&&P_{mm} \end{pmatrix} \]
Gauss 积分法
如果一个函数可以利用正交多项式的展开进行近似, 我们就可以利用这个展开计算定积分.
现在先给定积分区间\([a,b]\)以及一个内积空间\(\mathbb{V}\), 权函数为\(\omega\), 主要研究对象为函数\(f\in\mathbb{V}\). 内积空间上全体最高次系数为\(1\)的\(j\)次多项式, 构成了内积空间\(\mathbb{V}\)的真子集, 记作\(\mathbb{P}_j\).
我们不加证明地指出一个定理系统:
正交多项式与 Chebyshev 系统
正交基的存在性
存在一系列多项式\(p_k\in\mathbb{P}_k\ ,\ k=0,1,\cdots\), 满足正交关系 \[ (p_j,p_k)=0\ ,\ j\ne k \]
Gram-Schimidt 正交化方法
这个正交基可以由下列递推关系给出 \[ p_{k+1}(x)=(x-\alpha_{k+1})p_k(x)-\beta_{k+1}^2p_{k-1}(x) \newline p_0(x)\equiv1\ ,\ \text{规定 }p_{-1}(x)\equiv0 \] 其中, 参数\(\alpha\)和\(\beta\)分别为 \[ \alpha_{k+1}=\frac{(xp_{k+1},p_{k+1})}{(p_{k+1},p_{k+1})} \\ \beta_{k+1}^2=\frac{(p_k,p_k)}{(p_{k-1},p_{k-1})} \\ (\text{规定 }\beta^2_1=0) \]
零点的性质
\(p_k(x)\)的\(k\)个零点都是单重实根, 并且都在\([a,b]\)上
Haar 条件
给定\(n\)个两两不同的自变量\(\{t_i\}_{0\sim n-1}\), 下列系数矩阵是可逆的 \[ \begin{aligned} \boldsymbol{A}&=\begin{pmatrix} p_0(t_0) & \cdots & p_0(t_{n-1}) \newline \vdots & \cdots & \vdots \newline p_{n-1}(t_0) & \cdots & p_{n-1}(t_{n-1}) \end{pmatrix} \newline\newline (&\text{或者 }a_{jk}=p_{j-1}(t_{k-1})) \end{aligned} \]
Haar 条件意味着什么呢? 当我们计算一个插值问题时, 插值函数为 \[ P_{N}(x)=\sum_{k=0}^{N-1}c_kp_k(x) \] 支撑点为\(\{(t_i,f_i)\}_{0\sim N-1}\), 因此它满足一个线性方程组 \[ \sum_{k=0}^{N-1}p_k(t_{j})c_k=f_j \] 即 \[ \begin{pmatrix}c_0 & c_1 & \cdots & c_{N-1}\end{pmatrix}\cdot\boldsymbol{A}= \begin{pmatrix}f_0 & f_1 & \cdots & f_{N-1}\end{pmatrix} \] 它有唯一解的条件是\(\boldsymbol{A}\)具有full rank, 也就要求可逆, 即所谓的 Haar 条件保证了插值问题解的唯一性. 满足以上几个条件的一系列多项式, 称之为 Chebyshev 系统.
Gauss 点, 权重因子
我们把 Gram-Schimidt 正交化方法中的参数排成一个三对角矩阵 \[ \boldsymbol{J}_n=\begin{pmatrix} \alpha_1 & \beta_2 & \newline \beta_2 & \alpha_2 & \ddots \newline & \ddots & \ddots & \ddots \newline & & \ddots &\alpha_{n-1} & \beta_n \newline & & & \beta_n & \alpha_n \end{pmatrix} \] 定义函数 \[ p_n(x)=\det(x\boldsymbol{I}_n-\boldsymbol{J}_n) \] 明显有 \[ p_n(x)=(x-\alpha_n)p_{n-1}(x)-\beta_n^2p_{n-2}(x) \] 这明显就是之前定义的多项式基\(\{p_i(x)\}\), 而\(\boldsymbol{J}_n\)的本征值, 恰好是\(p_n(x)\)的\(n\)个零点.
假设我们选择了这\(n\)个零点作为积分的 Gauss 点: \[ I[f]\approx\sum_{j=1}^nw_jf(x_j) \] \(w_j\)称为 权重因子, 为了保证支撑点都有效, 并且不出现正负相消, 它们应该都为正. 考虑在正交多项式基上展开\(\displaystyle f(x)\approx\sum_{k=0}^{n-1}c_kp_k(x)\), 近似替代\(f\), 我们希望三个约等号联系起一个等号: \[ \sum_{k=0}^{n-1}c_k\int_a^b\omega(x)p_k(x)dx=\sum_{k=0}^{n-1}c_k\sum_{j=1}^nw_jp_k(x_j) \] 观察左边, 不难发现其实只有\(k=0\)项存在, 因为积分相当于\(p_0\equiv1\)和\(p_k\)的内积. 为了得到普适的式子(指对于任何\(\{c_i\}\)总能取等或者近似取等), 有 \[ \sum_{j=1}^nw_jp_k(x_j)=\begin{cases} \begin{aligned} \displaystyle(p_0,p_0)\qquad k=0& \newline\newline\displaystyle 0\qquad \text{otherwise}& \end{aligned} \end{cases} \]
Gauss 积分法总结
对于一个函数\(f\in\mathbb{C}^{2n}[a,b]\)和一个正交多项式基\(\{p_i\}\), 取\(p_n\)的零点\(\{x_i\}_{1\sim n}\)为 Gauss 点, 线性方程组 \[ \sum_{j=1}^np_k(x_j)w_j=\begin{cases} \begin{aligned} \displaystyle(p_0,p_0)\qquad k=0& \newline\newline\displaystyle 0\qquad \text{otherwise}& \end{aligned} \end{cases} \] 的解\(\vec{w}\)为权重因子, 则可以用如下公式近似替代定积分: \[ \int_a^b\omega(x)f(x)dx=\sum_{j=1}^nw_jf(x_j)+\frac{f^{(2n)}(\xi)}{(2n)!}(p_n,p_n) \] 后一项表示误差, \(\xi\)是\([a,b]\)上的某个数.