函数数值计算(04)：Clenshaw 逆向迭代, Padé 近似
Clenshaw 逆向迭代
真正用类似的展开来算目标函数值的时候, 我们可以采用 Clenshaw 提出的逆向迭代的算法. 先回顾对一个函数展开的近似方案.
首先有一个可展开的函数\(f(x)\), 以及一个基\(\{P_i(x)\}\)(不一定要求正交归一), 对于无穷维的基, 在\(i=N\)处截断, 用部分的基函数来近似描述\(f(x)\): \[ f(x)\approx S_N(x)=\sum_{k=0}^Nc_kP_k(x) \] 基函数具有一个至多二阶的递推关系: \[ P_{k+1}(x)=\alpha_k(x)P_k(x)+\beta_k(x)P_{k-1}(x) \] 现在, 假设我们已经得知\(\{c_i\}\), 目的是计算部分和函数\(S_N(x)\), 这个计算可以通过构造新的函数进行.
构造新函数
我们引入一类新的函数\(\{b_i(x)\}\), 并且规定 \[ b_{N+1}=b_{N+2}=0 \] 递推关系规定为 \[ b_k(x)=c_k+\alpha_k(x)b_{k+1}(x)+\beta_{k+1}(x)b_{k+2}(x) \]
重新计算部分和函数
我们顺势将部分和函数中的\(\{c_i\}\)全部替换, 分裂为三个和: \[ S_N(x)=\sum_{k=0}^Nb_k(x)P_k(x)-\sum_{k=0}^N\alpha_k(x)b_{k+1}(x)P_k(x)-\sum_{k=0}^N\beta_{k+1}(x)b_{k+2}(x)P_k(x) \] 基函数\(\{P_i(x)\}\)同样具有递推关系: \[ \begin{aligned} \beta_{k+1}(x)P_k(x)&=P_{k+2}(x)-\alpha_{k+1}(x)P_{k+1}(x) \newline&\Downarrow\newline \sum_{k=0}^N\beta_{k+1}(x)b_{k+2}(x)P_k(x)&=\sum_{k=0}^Nb_{k+2}(x)P_{k+2}(x)-\sum_{k=0}^N\alpha_{k+1}(x)b_{k+2}(x)P_{k+1}(x) \newline&= \sum_{k=2}^Nb_{k}(x)P_{k}(x)-\sum_{k=1}^N\alpha_{k}(x)b_{k+1}(x)P_{k}(x) \end{aligned} \] 代回部分和函数, 得到 \[ S_N(x)=b_0(x)P_0(x)+b_1(x)P_1(x)-\alpha_0(x)b_1(x)P_0(x) \] 我们通常不使用\(b_0(x)\)(这样就可以少迭代一阶了), 而是考虑\(b_0(x)=c_0+\alpha_0(x)b_1(x)+\beta_1(x)b_2(x)\), 得到 \[ \boxed{S_N(x)=c_0P_0(x)+b_1(x)P_1(x)+\beta_1(x)b_2(x)P_2(x)}\label{eq1}\tag{eq1} \]
方法总结
所谓 Clenshaw 逆向迭代法, 就是已知:
一个基\(\{P_i(x)\}\), 以及它的递推关系 \[ P_{k+1}(x)=\alpha_k(x)P_k(x)+\beta_k(x)P_{k-1}(x) \]
一个原函数\(f(x)\), 以及它基于某种方案在基\(\{P_i(x)\}\)上的有限维近似展开(前文多次称为 部分和函数)的系数\(\{c_i\}\) \[ f(x)\approx\sum_{k=0}^Nc_kP_k(x) \]
于是构造新函数 \[ b_k(x)=c_k+\alpha_k(x)b_{k+1}(x)+\beta_{k+1}(x)b_{k+2}(x) \newline b_{N+1}=b_{N+2}=0 \] 从而把\(N\)项的求和化简为少数几项和, 如\(\eqref{eq1}\): \[ S_N(x)=c_0P_0(x)+b_1(x)P_1(x)+\beta_1(x)b_2(x)P_2(x) \] 因为计算新函数\(\{b_i(x)\}\)时, 指标从大到小递推, 因此得名为"逆向迭代法".
经过反向递推, 我们把一个多项的求和式通过递推关系变成了少量几项的求和, 这一优点有些时候可以改进求和计算的收敛性和稳定性.
Padé 近似
还有一个常用的近似计算函数的方法是 Padé 近似, 也是当作级数展开的补充方法来用. 因为计算很多高阶的级数展开并不容易, Padé 近似相当于算了一些阶的级数展开以后来获得更精确结果的方法.
Padé 近似函数是一个分式函数的形式 \[ R(x)=\frac{\displaystyle\sum_{k=0}^Ma_kx^k}{1+\displaystyle\sum_{k=1}^Nb_kx^k} \] 我们希望它在展开点(例如\(x=0\))处与函数相等至\(M+N\)阶导数: \[ R^{(k)}(0)=c_k\newline (k=0,1,\cdots,M+N) \] 并称它为函数的 \((M,N)\)阶 Padé 近似. 其中, \(c_k=f^{(k)}(0)\), 或者说\(f(x)\)的幂级数为\(\displaystyle\sum_{k=0}^\infty c_kx^k\).
分式函数求导或许比较繁琐, 我们退而求其次, 令\(R(x)=\displaystyle\sum_{k=0}^\infty c_kx^k\), 并把分母乘到等式另一边, 再考虑各阶导数. 特别是\(M=N\)的情形, 得到 \[ \begin{aligned} \sum_{m=1}^Nc_{N-m+k}b_m&=-c_{N+k} \newline a_k&=\sum_{m=0}^kc_{k-m}b_m \newline (k&=1,\cdots,N) \end{aligned}\label{eq2}\tag{eq2} \] \(\eqref{eq2}\)上式是一个\(N\)元线性方程组, 可解出\(\{b_i\}\), 再代入\(\eqref{eq2}\)下式, 计算出\(\{a_i\}\). 显然, 我们需要在知道\(\{c_i\}_{0\sim2N}\)的前提下进行求解.