3 Interpolation Methods

3.1 Lagrange Interpolation

考虑拟合一个函数, 其中 \(n+1\) 个互异节点的函数值已知:

\[y_0=f(x_0),\quad y_1=f(x_1),\quad \cdots,\quad y_n=f(x_n)\]

Lagrange interpolation 要求拟合函数 \(P(x)\) 的每个点都经过插值点:

\[P(x_k)=f(x_k),\quad k=0,1,2,\cdots,n\]

Lagrange interpolation 函数 \(P(x)\):

\[P(x)=\sum_{k=0}^n f(x_k)L_{n,k}(x)\]

插值基函数:

\[\begin{aligned} L_{n,k}(x) &= \frac{(x-x_0)(x-x_1)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)} {(x_k-x_0)(x_k-x_1)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}\\ &= \prod_{0\leq i\leq n,\, i\neq k} \frac{x-x_i}{x_k-x_i}. \end{aligned}\]

可以看到

\[L_{n,k}(x_i)= \begin{cases} 1, & k = i,\\ 0, & k \neq i. \end{cases}\]

若 \(f\in C^{n+1}\) 于包含插值节点和 \(x\) 的闭区间, 插值余项为

\[R_n(x)=f(x)-P_n(x) =\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n}(x-x_i),\]

其中 \(\xi\) 位于这些点之间. 误差估计通常写成

\[|R_n(x)| \leq \frac{\max_\xi |f^{(n+1)}(\xi)|}{(n+1)!} \left|\prod_{i=0}^{n}(x-x_i)\right|.\]

Example \(f(x)=x^{-1}\), 计算插值误差. 插值点为 \(x_0=2,x_1=2.75,x_2=4\). 相应函数值为

\[f(2)=\frac12,\qquad f(2.75)=\frac{4}{11},\qquad f(4)=\frac14.\]

Solution 三个 Lagrange 基函数:

\[\begin{aligned} L_0(x) &=\frac{(x-2.75)(x-4)}{(2-2.75)(2-4)} =\frac{(x-\frac{11}{4})(x-4)}{\frac32},\\ L_1(x) &=\frac{(x-2)(x-4)}{(2.75-2)(2.75-4)} =-\frac{16}{15}(x-2)(x-4),\\ L_2(x) &=\frac{(x-2)(x-2.75)}{(4-2)(4-2.75)} =\frac{2}{5}(x-2)(x-\frac{11}{4}). \end{aligned}\]

因此

\[\begin{aligned} P(x) &=\frac12 L_0(x)+\frac{4}{11}L_1(x)+\frac14 L_2(x)\\ &=\frac{1}{3}(x-\frac{11}{4})(x-4) -\frac{64}{165}(x-2)(x-4) +\frac{1}{10}(x-2)(x-\frac{11}{4})\\ &=\frac{1}{22}x^2-\frac{35}{88}x+\frac{49}{44} \end{aligned}\]

这里有三个插值点, 所以 \(n=2\), 需要三阶导数. 因为 \(f^{(3)}(x)=-6x^{-4}\), 二次插值的余项为

\[R(x) = \frac{f^{(3)}(\xi)}{3!}(x-x_0)(x-x_1)(x-x_2) = -\xi^{-4}(x-2)(x-2.75)(x-4).\]

在 \(x,\xi\in [2,4]\) 上, 对右侧两项取极值得到

\[|R(x)|\leq \frac{1}{16}\cdot\frac{9}{16} =\frac{9}{256}\approx 0.03516.\]正文倒数函数例题的二次 Lagrange 插值曲线和实际误差

图 3.1 采用本节倒数函数例题的三个节点 2、2.75、4. 左图叠加原函数与二次插值多项式, 黑点表示插值节点; 右图单独显示原函数减插值多项式的误差, 三个节点处误差为零. 右图纵轴是实际误差, 不是前文给出的统一误差上界. 由本节例题直接计算绘制.

Example 证明

\[\sum_{k=0}^{n}\frac{x_k^n}{(x_k-x_0)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}=1\]

Solution 对 \(p(x)=x^n\) 作 \(n\) 次 Lagrange 插值. 因 \(p\) 为 \(n\) 次多项式, 插值多项式即为自身

\[x^n = \sum_{k=0}^{n} x_k^n L_k(x), \qquad L_k(x)=\frac{\prod_{j\ne k}(x-x_j)}{\prod_{j\ne k}(x_k-x_j)}\]

比较两边 \(x^n\) 的最高次项系数. 左边为 1, \(L_k(x)\) 的最高次项系数为

\[\frac{1}{\prod_{j\ne k}(x_k-x_j)}\]

因此

\[\sum_{k=0}^{n} \frac{x_k^n}{\prod_{j\ne k}(x_k-x_j)} = 1\]

3.2 Hermite Interpolation *

Hermite interpolation 是 Lagrange interpolation 的升级版, 不仅要求拟合函数和原函数在插值点的函数值相同, 还要求在插值点的斜率相同. 插值函数为

\[H_{2n+1}(x)=\sum_{i=0}^n f(x_i)H_{n,i}(x)+\sum_{i=0}^n f'(x_i)\hat H_{n,i}(x)\]

\(H_{n,i}(x)\) 为了处理函数值:

\[\begin{aligned} H_{n,i}(x) &= \left[1-2(x-x_i)L'_{n,i}(x_i)\right]L_{n,i}^2(x),\\ H_{n,i}(x_i) &= 1,\qquad H'_{n,i}(x_i)=0. \end{aligned}\]

\(\hat{H}_{n,i}\left(x\right)\) 为了处理导数:

\[\begin{aligned} \hat H_{n,i}(x) &= (x-x_i)L_{n,i}^2(x),\\ \hat H_{n,i}(x_i) &= 0,\qquad \hat H'_{n,i}(x_i)=1. \end{aligned}\]

若节点互异且 \(f\in C^{2n+2}\) 于包含节点和 \(x\) 的闭区间, Hermite 插值余项为:

\[f(x)=H_{2n+1}(x)+\frac{(x-x_0)^2\cdots(x-x_n)^2}{(2n+2)!}f^{(2n+2)}(\xi)\]

3.3 Cubic Spline Interpolation

Runge phenomenon 是数值分析中的一个经典问题, 指的是在等距节点下, 使用高次多项式插值时, 插值多项式在区间两端会出现剧烈振荡, 导致插值误差反而比低次多项式更大, 甚至完全不收敛于原函数.

想要减小 Runge phenomenon, 可以使用低次分段插值方法, 比较经典的方法是三次样条插值.

函数有 \(n+1\) 个数据点, 节点满足 \(x_0\lt x_1\lt\cdots\lt x_n\), 记 \(h_i=x_{i+1}-x_i>0\), 分为 \(n\) 个区间. 考虑第 \(i\) 个区间, 用三次多项式函数拟合:

\[S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3,\quad i=0,1,\cdots,n-1\]

一共有 \(4n\) 个系数.

  • 插值条件

    \[\begin{aligned} S_i(x_i)=f_i &\implies a_i=f_i,\\ S_i(x_{i+1})=f_{i+1} &\implies a_i+b_i h_i+c_i h_i^2+d_i h_i^3=f_{i+1}. \end{aligned}\]

    给出了 \(2n\) 个方程, 若消去 \(f\) 可以得到

    \[a_{i+1} = a_i+b_i h_i+c_i h_i^2+d_i h_i^3\]

  • 一阶导数连续

    \[S_i'(x_{i+1}) = S_{i+1}'(x_{i+1}) \implies b_i+2c_i h_i+3d_i h_i^2 = b_{i+1}.\]

    给出 \(n-1\) 个方程.

  • 二阶导数连续

    \[S_i''(x_{i+1}) = S_{i+1}''(x_{i+1}) \implies c_i+3d_i h_i = c_{i+1}.\]

    给出 \(n-1\) 个方程.

  • 边界条件: 要根据实际情况来确定. 常用的是 (1) 自然样条边界条件: 两端二阶导数为 \(0\).

    \[\begin{aligned} S_0''(x_0)=0 &\implies c_0=0,\\ S_{n-1}''(x_n)=0 &\implies c_{n-1}+3d_{n-1}h_{n-1}=0. \end{aligned}\]

    (2) 固定边界条件:

    \[\begin{aligned} S_0'(x_0)=f_0' &\implies b_0=f_0',\\ S_{n-1}'(x_n)=f_n' &\implies b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=f_n'. \end{aligned}\]

3.3.1 Coefficient Relations

已知方程:

\[\begin{aligned} a_{i+1} &= a_i+b_i h_i+c_i h_i^2+d_i h_i^3,\\ b_{i+1} &= b_i+2c_i h_i+3d_i h_i^2,\\ c_{i+1} &= c_i+3d_i h_i. \end{aligned}\]

消去 \(d_{i}\) 和 \(b_{i}\), 过程中得到

\[d_i=\frac{c_{i+1}-c_i}{3h_i},\qquad b_i=\frac{a_{i+1}-a_i}{h_i}-\frac{h_i}{3}(2c_i+c_{i+1})\]

结果为

\[h_{i-1}c_{i-1}+2(h_{i-1}+h_i)c_i+h_i c_{i+1} =\frac{3}{h_i}(a_{i+1}-a_i)-\frac{3}{h_{i-1}}(a_i-a_{i-1})\]

上述内部节点方程适用于 \(i=1,\cdots,n-1\). 这里把 \(a_n=f_n\) 和 \(c_n=S_{n-1}^{\prime\prime}(x_n)/2\) 作为右端点的记号, 并非再增加一个区间多项式. 已知 \(a_i=f_i\), 加上两个边界条件即可求解 \(c_i\).

3.3.2 Boundary Conditions

  1. 自然样条: \(S''(x_0)=S''(x_n)=0\), 令 \(\mathbf{c}=(c_0,c_1,\cdots,c_n)^{\mathsf{T}}\). 自然样条对应的三对角线性方程组为

    \[\begin{cases} c_0=0,\\ h_{i-1}c_{i-1}+2(h_{i-1}+h_i)c_i+h_i c_{i+1} =\dfrac{3}{h_i}(a_{i+1}-a_i)-\dfrac{3}{h_{i-1}}(a_i-a_{i-1})\\ c_n=0. \end{cases}\]

  2. 固定边界条件: 这里 \(a=x_0\), \(b=x_n\) 表示区间端点, 并给定端点一阶导数:

    \[S'(x_0)=f'(a),\qquad S'(x_n)=f'(b).\]

    左端点 \(x_0\): 利用 \(S_0(x_1) = a_1\) 且 \(S_0^{\prime}(x_0) = f^{\prime}(a)\)

    \[b_0 =\frac{a_1-a_0}{h_0}-\frac{h_0}{3}(2c_0+c_1), \qquad b_0=f'(a),\]

    消去 \(b_0\) 得到

    \[2h_0c_0+h_0c_1 = \frac{3}{h_0}(a_1-a_0)-3f'(a).\]

    右端点 \(x_n\): 利用 \(S_{n-1}(x_n)=a_n\) 且 \(S_{n-1}'(x_n)=f'(b)\). 由

    \[a_n=a_{n-1}+b_{n-1}h_{n-1}+c_{n-1}h_{n-1}^2+d_{n-1}h_{n-1}^3\]

    以及 \(d_{n-1}=(c_n-c_{n-1})/(3h_{n-1})\), 得

    \[b_{n-1} =\frac{a_n-a_{n-1}}{h_{n-1}} -\frac{h_{n-1}}{3}(2c_{n-1}+c_n), \qquad b_{n-1}+h_{n-1}(c_{n-1}+c_n)=f'(b).\]

    消去 \(b_{n-1}\) 得到

    \[h_{n-1}c_{n-1}+2h_{n-1}c_n = 3f'(b)-\frac{3}{h_{n-1}}(a_n-a_{n-1}).\]

    于是固定边界条件对应的三对角线性方程组为

    \[\begin{cases} 2h_0c_0+h_0c_1 =\dfrac{3}{h_0}(a_1-a_0)-3f'(a),\\ h_{i-1}c_{i-1}+2(h_{i-1}+h_i)c_i+h_i c_{i+1} =\dfrac{3}{h_i}(a_{i+1}-a_i)-\dfrac{3}{h_{i-1}}(a_i-a_{i-1}), \quad i=1,\cdots,n-1,\\ h_{n-1}c_{n-1}+2h_{n-1}c_n =3f'(b)-\dfrac{3}{h_{n-1}}(a_n-a_{n-1}). \end{cases}\]

解出 \(c_i\) 后, 由

\[d_i=\frac{c_{i+1}-c_i}{3h_i},\qquad b_i=\frac{a_{i+1}-a_i}{h_i}-\frac{h_i}{3}(2c_i+c_{i+1})\]

恢复每个区间上的 \(S_i(x)\).

3.4 Akima Interpolation *

Akima 插值用局部数据确定节点斜率, 再构造分段三次 Hermite 多项式. 它保证函数及一阶导数连续, 一般不保证二阶导数连续, 因而与上一节的三次样条条件不同.1

\[S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3\]

需要四个方程来确定四个系数:

\[\begin{aligned} S_i(x_i)&=f_i,\qquad S_i(x_{i+1})=f_{i+1},\\ S_i'(x_i)&=t_i,\qquad S_i'(x_{i+1})=t_{i+1}. \end{aligned}\]

定义割线斜率和权重:

\[\begin{aligned} m_i &= \frac{f_{i+1}-f_i}{x_{i+1}-x_i},\\ w_{i-1} &= |m_{i+1}-m_i|,\\ w_i &= |m_{i-1}-m_{i-2}|. \end{aligned}\]

\(w_{i-1}\) 是右侧差分, \(w_i\) 是左侧差分. 节点一阶导数 \(t_i=S'(x_i)\) 取为

\[t_i=\frac{w_{i-1}m_{i-1}+w_i m_i}{w_{i-1}+w_i}\]

一些其他的代替方案:

  • \(m_{i-2}=m_{i-1},m_i\neq m_{i+1}\) 时, 取左边的稳定坡度值

    \[t_i=m_{i-1}\]

  • \(m_i=m_{i+1},m_{i-2}\neq m_{i-1}\) 时, 取右侧的稳定坡度值

    \[t_i=m_i\]

  • \(m_{i-1}=m_i\) 时, 取共同的坡度值

    \[t_i=m_{i-1}=m_i\]

  • \(m_{i-2}=m_{i-1}\neq m_i=m_{i+1}\) 时, 取平均值

    \[t_i=\frac{m_{i-1}+m_i}{2}\]

节点斜率 \(t_i\) 用到四个割线斜率 \(m_{i-2},m_{i-1},m_i,m_{i+1}\), 涉及五个节点. 当两个权重均为零时, 采用上面的共同坡度或平均值规则. 对至少三个节点, 边界缺少的割线斜率可按下式外推; 只有两个节点时可直接采用线性插值:

\[\begin{aligned} m_{-1} &= 2m_0-m_1,\\ m_{-2} &= 2m_{-1}-m_0,\\ m_n &= 2m_{n-1}-m_{n-2},\\ m_{n+1} &= 2m_n-m_{n-1}. \end{aligned}\]

Example 拟合函数

\[y=a+b(x-x_1)+c(x-x_1)^2+d(x-x_1)^3,\quad x\in[x_1,x_2]\]

边界条件:

\[\begin{aligned} y(x_1)&=y_1,\qquad y(x_2)=y_2,\\ y'(x_1)&=t_1,\qquad y'(x_2)=t_2. \end{aligned}\]

Solution 求解出系数:

\[\begin{aligned} a &= y_1,\\ b &= t_1,\\ c &= \frac{1}{x_2-x_1}\left[\frac{3(y_2-y_1)}{x_2-x_1}-2t_1-t_2\right],\\ d &= \frac{1}{(x_2-x_1)^2}\left[t_1+t_2-\frac{2(y_2-y_1)}{x_2-x_1}\right]. \end{aligned}\]

3.5 Rational Function Interpolation *

分数多项式插值可通过分母描述某些极点结构, 在合适的函数及节点上可能比多项式插值更准确. 但它也可能引入伪极点, 不能普遍保证避免 Runge phenomenon 或优于多项式插值. 一种方法是经过给定数据点, 即有理函数插值. 已知函数的 \(m+1\) 个点的数值 \((x_i, y_i)\), \(i=1,2,\cdots,m+1\). 插值函数:

\[R_{i(i+1)\cdots(i+m)}(x)=\frac{P_\mu(x)}{Q_\nu(x)} =\frac{p_0+p_1x+\cdots+p_\mu x^\mu}{1+q_1x+\cdots+q_\nu x^\nu}\]

其中 \(m=\mu+\nu\), 其数值应根据实际情况选取, \(P_{\mu}(x)\), \(Q_{\nu}(x)\) 为两个多项式. 确定系数:

\[\begin{aligned} R(x_i)&=\frac{P_\mu(x_i)}{Q_\nu(x_i)}=y_i,\\ P_\mu(x_i)-y_i Q_\nu(x_i)&=0,\\ \sum_{j=0}^{\mu}p_j x_i^j-y_i\left(1+\sum_{k=1}^{\nu}q_k x_i^k\right)&=0. \end{aligned}\]

这是关于未知系数的线性方程组. 需要检查该归一化下方程组是否有解, 并确保 \(Q_\nu(x_i)\ne0\); 交叉相乘所得等式本身不排除节点处 \(P_\mu=Q_\nu=0\) 的未定义情形:

\[\begin{aligned} p_0+x_1p_1+\cdots+x_1^\mu p_\mu-y_1x_1q_1-\cdots-y_1x_1^\nu q_\nu &= y_1,\\ p_0+x_2p_1+\cdots+x_2^\mu p_\mu-y_2x_2q_1-\cdots-y_2x_2^\nu q_\nu &= y_2,\\ &\cdots\\ p_0+x_{m+1}p_1+\cdots+x_{m+1}^\mu p_\mu-y_{m+1}x_{m+1}q_1-\cdots-y_{m+1}x_{m+1}^\nu q_\nu &= y_{m+1}. \end{aligned}\]

3.6 Padé Approximation *

Padé 近似是另一种分式多项式插值方法, 已知函数 \(f(x)\) 的 Taylor expansion, 我们要构造有理函数, 使其 Taylor expansion 在 \(x=0\) 处与 \(f(x)\) 匹配到尽可能高的阶数. 拟合函数:

\[R(x)=\frac{a_0+a_1x+\cdots+a_m x^m}{1+b_1x+\cdots+b_n x^n}\]

\(f(x)\) 在原点的 Taylor 系数已知, 或给定相应形式幂级数. 约定 \(c_k=0\) 当 \(k\lt0\), 并令 \(b_0=1\). 下述系数匹配以这个归一化的线性方程组存在解为前提:

\[f(x)=\sum_{k=0}^{\infty}c_k x^k\]

要求两个展开的前 \(m+n+1\) 项完全相同, 即

\[\begin{aligned} \frac{a_0+a_1x+\cdots+a_m x^m}{1+b_1x+\cdots+b_nx^n} &=c_0+c_1x+\cdots+c_{m+n}x^{m+n}+\cdots,\\ (1+b_1x+\cdots+b_nx^n)(c_0+c_1x+c_2x^2+\cdots) &=a_0+a_1x+\cdots+a_m x^m+O(x^{m+n+1}). \end{aligned}\]

各阶系数一般满足

\[a_k=\sum_{i=0}^{\min(k,n)}b_i c_{k-i},\qquad 0\leq k\leq m.\]

下面的一阶、二阶示例只在相应次数的系数存在时使用; 若某个 \(b_i\) 超过分母次数, 则取其为零.

零阶项:

\[a_0=c_0\]

一阶项:

\[a_1=c_1+b_1c_0\]

二阶项:

\[\begin{aligned} a_2 &= c_2+b_1c_1+b_2c_0,\\ &\cdots\\ a_m &= c_m+\sum_{i=1}^{\min(m,n)}b_i c_{m-i}. \end{aligned}\]

\(x^{m+1}\) 到 \(x^{m+n}\) 的系数都为零, 因为右边没有 \(m+1\) 阶以上的多项式项.

\[\begin{pmatrix} c_{m-n+1} & c_{m-n+2} & \cdots & c_m\\ c_{m-n+2} & c_{m-n+3} & \cdots & c_{m+1}\\ \vdots & \vdots & \ddots & \vdots\\ c_m & c_{m+1} & \cdots & c_{m+n-1} \end{pmatrix} \begin{pmatrix} b_n\\ b_{n-1}\\ \vdots\\ b_1 \end{pmatrix} =- \begin{pmatrix} c_{m+1}\\ c_{m+2}\\ \vdots\\ c_{m+n} \end{pmatrix}\]

4 Numerical Differentiation and Integration

4.1 First-Derivative Formulas

若 \(h>0\) 且 \(f\in C^2[x_0,x_0+h]\), 最简单的一阶微分公式是向前差分公式, 其中 \(\xi\in(x_0,x_0+h)\):

\[f'(x_0)=\frac{f(x_0+h)-f(x_0)}{h}-\frac{h}{2}f''(\xi)\]

截断误差精度为一阶, 误差相对较大. 考虑其他微分公式来提高精度.

取三个互异节点 \(x_0,x_1,x_2\), 令

\[P_2(x)=\sum_{k=0}^{2}f(x_k)L_k(x),\qquad \omega_3(x)=\prod_{k=0}^{2}(x-x_k).\]

在节点 \(x_i\) 处对插值多项式求导, 得到三点微分近似. 若 \(f\in C^3\) 于包含节点的区间, 则完整公式为

\[\begin{aligned} f'(x_i) &=f(x_0)\frac{2x_i-x_1-x_2}{(x_0-x_1)(x_0-x_2)} +f(x_1)\frac{2x_i-x_0-x_2}{(x_1-x_0)(x_1-x_2)}\\ &\quad+f(x_2)\frac{2x_i-x_0-x_1}{(x_2-x_0)(x_2-x_1)}+R(x_i),\\ R(x_i)&=\frac{f^{(3)}(\xi_i)}{6}\omega_3'(x_i) =\frac{f^{(3)}(\xi_i)}{6}\prod_{\substack{k=0\\k\ne i}}^2(x_i-x_k). \end{aligned}\]

这里不能把插值余项中的存在性记号 \(\xi\) 直接当作可微函数求导. 余项可独立验证: 取

\[K=\frac{f'(x_i)-P_2'(x_i)}{\omega_3'(x_i)},\qquad \phi(x)=f(x)-P_2(x)-K\omega_3(x).\]

\(\phi\) 在三个节点为零且 \(\phi'(x_i)=0\). 把 \(x_i\) 的重零点计入后, 反复使用 Rolle 定理可得某个 \(\xi_i\) 满足 \(\phi^{(3)}(\xi_i)=0\). 因 \(P_2^{(3)}=0\)、\(\omega_3^{(3)}=6\), 所以 \(K=f^{(3)}(\xi_i)/6\).2

若等间距取点, 即 \(x_1=x_0+h,x_2=x_0+2h\), 得到三个点的微分公式:

\[\begin{aligned} f'(x_0)&=\frac{1}{h}\left[-\frac{3}{2}f(x_0)+2f(x_0+h)-\frac{1}{2}f(x_0+2h)\right]+\frac{h^2}{3}f^{(3)}(\xi_0),\\ f'(x_0+h)&=\frac{1}{h}\left[-\frac{1}{2}f(x_0)+\frac{1}{2}f(x_0+2h)\right]-\frac{h^2}{6}f^{(3)}(\xi_1),\\ f'(x_0+2h)&=\frac{1}{h}\left[\frac{1}{2}f(x_0)-2f(x_0+h)+\frac{3}{2}f(x_0+2h)\right]+\frac{h^2}{3}f^{(3)}(\xi_2). \end{aligned}\]

平移得到

  • 端点三点微分公式:

    \[f'(x_0)=\frac{1}{2h}\left[-3f(x_0)+4f(x_0+h)-f(x_0+2h)\right]+\frac{h^2}{3}f^{(3)}(\xi_0)\]

  • 中心差分公式:

    \[f'(x_0)=\frac{1}{2h}\left[-f(x_0-h)+f(x_0+h)\right]-\frac{h^2}{6}f^{(3)}(\xi_1)\]

  • 端点三点公式:

    \[f'(x_0)=\frac{1}{2h}\left[f(x_0-2h)-4f(x_0-h)+3f(x_0)\right]+\frac{h^2}{3}f^{(3)}(\xi_2)\]

4.2 Second-Derivative Formulas

若 \(f\in C^4\) 于包含 \(x_0-h,x_0+h\) 的区间, 在左右对称点展开

\[\begin{aligned} f(x_0+h)&=f(x_0)+f'(x_0)h+\frac{1}{2}f''(x_0)h^2+\frac{1}{6}f'''(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_1)h^4,\\ f(x_0-h)&=f(x_0)-f'(x_0)h+\frac{1}{2}f''(x_0)h^2-\frac{1}{6}f'''(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_{-1})h^4 \end{aligned}\]

相加消去 Taylor 展开中的奇次导数项

\[f''(x_0)=\frac{1}{h^2}\left[f(x_0-h)-2f(x_0)+f(x_0+h)\right]-\frac{h^2}{24}\left[f^{(4)}(\xi_1)+f^{(4)}(\xi_{-1})\right].\]

因此记二阶中心差分公式为

\[f''(x_0) = \frac{f(x_0-h)-2f(x_0)+f(x_0+h)}{h^2} +O(h^2).\]

4.3 Simpson Formula

由 Lagrange interpolation 公式:

\[\begin{aligned} P_n(x)&=\sum_{i=0}^{n}f(x_i)L_i(x),\\ \int_a^b f(x)\,\mathrm{d}x &=\int_a^b\sum_{i=0}^{n}f(x_i)L_i(x)\,\mathrm{d}x +\frac{1}{(n+1)!}\int_a^b\prod_{i=0}^{n}(x-x_i)f^{(n+1)}(\xi)\,\mathrm{d}x\\ &=\sum_{i=0}^{n}a_i f(x_i)+\frac{1}{(n+1)!}\int_a^b\prod_{i=0}^{n}(x-x_i)f^{(n+1)}(\xi)\,\mathrm{d}x\\ &\approx \sum_{i=0}^{n}a_i f(x_i). \end{aligned}\]

其中

\[a_i=\int_a^b L_i(x)\,\mathrm{d}x\]

选三个等距点:

\[x_0=a,\qquad x_1=\frac{a+b}{2},\qquad x_2=b\]

插值多项式为

\[P_2(x)=L_0(x)f(x_0)+L_1(x)f(x_1)+L_2(x)f(x_2)\]

因此积分可写作:

\[\begin{aligned} \int_a^b P_2(x)\,\mathrm{d}x &=\sum_{k=0}^{2}a_k f(x_k)\\ &=a_0f(x_0)+a_1f(x_1)+a_2f(x_2),\\ L_k(x)&=\prod_{i=0,i\neq k}^{2}\frac{x-x_i}{x_k-x_i},\\ a_0&=\int_a^b L_0(x)\,\mathrm{d}x=\frac{1}{6}(b-a). \end{aligned}\]\[\begin{aligned} a_1&=\int_a^b L_1(x)\,\mathrm{d}x=\frac{2}{3}(b-a),\\ a_2&=\int_a^b L_2(x)\,\mathrm{d}x=\frac{1}{6}(b-a). \end{aligned}\]

得到 Simpson 公式:

\[\int_a^b f(x)\,\mathrm{d}x\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]\]

若 \(f\in C^4[a,b]\), 其余项为

\[R=\int_a^b f(x)\,\mathrm{d}x-\int_a^b P_2(x)\,\mathrm{d}x =-\frac{(b-a)^5}{2880}f^{(4)}(\xi),\qquad \xi\in(a,b).\]

因此 Simpson 公式对三次及以下多项式精确, 比二次插值本身的次数多一阶.

4.4 Newton-Cotes Formulas

取 \(x_k=a+kh\), \(k=0,\cdots,n\), 共 \(n+1\) 个等距节点, 积分公式为

\[\begin{aligned} \int_a^b f(x)\,\mathrm{d}x&\approx \sum_{k=0}^{n}a_k f(x_k),\\ a_k&=\int_a^b L_k(x)\,\mathrm{d}x =\int_a^b\prod_{i=0,i\ne k}^{n}\frac{x-x_i}{x_k-x_i}\,\mathrm{d}x. \end{aligned}\]

定义科特斯系数:

\[c_k=\frac{a_k}{b-a}\]

步长:

\[\begin{aligned} h&=\frac{b-a}{n},\\ \int_a^b f(x)\,\mathrm{d}x &\approx (b-a)\sum_{k=0}^{n}c_k f(x_k) =(b-a)\sum_{k=0}^{n}c_k f(a+kh),\\ c_k &=\frac{1}{b-a}\int_a^b\prod_{i=0,i\ne k}^{n}\frac{x-x_i}{x_k-x_i}\,\mathrm{d}x\\ &=\frac{1}{nh}\int_0^n\prod_{i=0,i\ne k}^{n}\frac{\lambda-i}{k-i}\,h\,\mathrm{d}\lambda\\ &=\frac{(-1)^{n-k}}{n(n-k)!k!}\int_0^n\prod_{i=0,i\ne k}^{n}(\lambda-i)\,\mathrm{d}\lambda. \end{aligned}\]

科特斯系数的性质:

  • 对于确定的 \(n\), 有

    \[\sum_{k=0}^{n}c_k=1\]

  • 对称性

    \[c_k=c_{n-k}\]

  • 闭型 Newton-Cotes 公式从 \(n=8\) 开始会出现负权重. 负权重会增加抵消及数据误差放大的风险, 因而实践中通常采用低阶复合公式; 这不意味着每个高阶公式对每个函数都必然不稳定.

比较常用的是 \(n=4\) 的闭型 Newton-Cotes 公式(Boole 公式). 取 \(h=(b-a)/4\), 若 \(f\in C^6[a,b]\), 可写出求积式及余项:

\[\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &\approx (b-a)\left[\frac{7}{90}f(x_0)+\frac{16}{45}f(x_1)+\frac{2}{15}f(x_2)+\frac{16}{45}f(x_3)+\frac{7}{90}f(x_4)\right],\\ R&=-\frac{8h^7}{945}f^{(6)}(\xi). \end{aligned}\]

4.5 Composite Quadrature Formulas

复合积分公式将积分区间分为若干子区间, 分别利用低阶积分公式, 再把各段积分加起来. 考虑复合 Simpson 公式: 令 \(n\) 为正偶数, \(h=(b-a)/n\), \(x_k=a+kh\), 并假设 \(f\in C^4[a,b]\). 每个 Simpson 面板包含两个长度为 \(h\) 的子区间:

\[\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &=\sum_{k=1}^{n/2}\int_{x_{2k-2}}^{x_{2k}}f(x)\,\mathrm{d}x\\ &=\sum_{k=1}^{n/2}\left\{\frac{h}{3}\left[f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})\right]-\frac{h^5}{90}f^{(4)}(\xi_k)\right\} \end{aligned}\]

每个 \(\xi_k\) 位于对应面板内. 利用 \(f^{(4)}\) 的连续性合并余项, 同时合并重复节点的权重:

\[\int_a^b f(x)\,\mathrm{d}x =\frac{h}{3}\left[f(a)+2\sum_{k=1}^{n/2-1}f(x_{2k})+4\sum_{k=1}^{n/2}f(x_{2k-1})+f(b)\right] -\frac{b-a}{180}h^4 f^{(4)}(\xi')\]

4.6 Remainder Extrapolation

假设精确值为 \(M\), 用步长 \(h\) 计算近似值为 \(N_1(h)\). 下面以余项具有幂次渐近展开、且系数 \(k_i\) 与 \(h\) 无关为前提; 实际计算按已有光滑性截断展开, 不要求无穷级数收敛:

\[M=N_1(h)+k_1h+k_2h^2+k_3h^3+\cdots\]

换成一半步长:

\[\begin{aligned} M&=N_1\left(\frac{h}{2}\right)+k_1\frac{h}{2}+k_2\frac{h^2}{4}+k_3\frac{h^3}{8}+\cdots,\\ M&=2N_1\left(\frac{h}{2}\right)-N_1(h)-\frac{k_2}{2}h^2-\frac{3k_3}{4}h^3+\cdots. \end{aligned}\]

定义:

\[\begin{aligned} N_2(h)&=2N_1\left(\frac{h}{2}\right)-N_1(h),\\ M&=N_2(h)-\frac{k_2}{2}h^2-\frac{3k_3}{4}h^3+\cdots. \end{aligned}\]

精度从 \(O(h)\) 提高到了 \(O(h^{2})\).

下面考虑只有偶数阶余项的情形:

\[\begin{aligned} M&=N_1(h)+k_1h^2+k_2h^4+k_3h^6+\cdots,\\ M&=N_1\left(\frac{h}{2}\right)+k_1\frac{h^2}{4}+k_2\frac{h^4}{16}+k_3\frac{h^6}{64}+\cdots,\\ 3M&=4N_1\left(\frac{h}{2}\right)-N_1(h)-\frac{3k_2}{4}h^4-\frac{15k_3}{16}h^6+\cdots. \end{aligned}\]

定义:

\[\begin{aligned} N_2(h)&=\frac{1}{3}\left[4N_1\left(\frac{h}{2}\right)-N_1(h)\right],\\ M&=N_2(h)-\frac{k_2}{4}h^4-\frac{5k_3}{16}h^6+\cdots. \end{aligned}\]

下一阶:

\[\begin{aligned} N_3(h)&=\frac{1}{15}\left[16N_2\left(\frac{h}{2}\right)-N_2(h)\right],\\ M&=N_3(h)+\frac{k_3}{64}h^6+\cdots. \end{aligned}\]

考虑复合梯形公式:

\[\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &=\frac{h}{2}\left[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right]+k_1h^2+k_2h^4+k_3h^6+\cdots,\\ N_1(h)&=\frac{h}{2}\left[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right]. \end{aligned}\]

按相同方式逐阶消去偶次误差项, 统一递推为

\[N_i(h)=N_{i-1}\left(\frac{h}{2}\right)+\frac{N_{i-1}\left(\frac{h}{2}\right)-N_{i-1}(h)}{4^{i-1}-1}.\]

4.7 Romberg Integration

Romberg algorithm 是一种加速算法, 可以得到更高精度的求积公式. 考虑变步长梯形法, 将积分区间 \([a,b]\) 等分为 \(n\) 个区间, 有 \(n+1\) 个点.

\[x_k=a+kh,\quad h=\frac{b-a}{n},\quad k=0,1,2,\cdots,n\]

积分:

\[T_n=\sum_{k=0}^{n-1}\frac{h}{2}\left[f(x_k)+f(x_{k+1})\right]\]

对每一个子区间再加一个节点 \(x_{k+1/2}\), 则每一段区间的积分为

\[\frac{h}{4}\left[f(x_k)+2f(x_{k+1/2})+f(x_{k+1})\right]\]

整个区间:

\[\begin{aligned} T_{2n} &=\frac{h}{4}\sum_{k=0}^{n-1}\left[f(x_k)+f(x_{k+1})\right] +\frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+1/2})\\ &=\frac{T_n}{2}+\frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+1/2}). \end{aligned}\]

等间隔取点:

\[h=\frac{b-a}{n},\quad x_{k+1/2}=a+\left(k+\frac{1}{2}\right)h\]

截断误差:

\[\begin{aligned} I-T_n&=-\frac{b-a}{12}h^2f''(\xi_n) =-\frac{b-a}{12}\left(\frac{b-a}{n}\right)^2f''(\xi_n),\\ I-T_{2n}&=-\frac{b-a}{12}\left(\frac{b-a}{2n}\right)^2f''(\xi_{2n}). \end{aligned}\]

不同步长对应的 \(\xi_n\) 一般不同, 不能直接令两个余项严格相差四倍. 若 \(f\) 足够光滑, 用 Euler-Maclaurin 型渐近展开

\[I-T_n=k_1h^2+k_2h^4+k_3h^6+\cdots\]

消去最低阶项, 得

\[I-T_{2n}=\frac{T_{2n}-T_n}{3}+O(h^4),\qquad S_n\equiv\frac{4T_{2n}-T_n}{3},\qquad I=S_n+O(h^4).\]

仅在该渐近展开适用且主导项未因抵消消失时, 相邻近似值之差才是可靠的余项估计; 两个近似值接近本身不能保证真误差很小.

这里 \(S_n\) 是由 \(n\) 个宽度为 \(h\) 的 Simpson 面板组成的复合公式, 其节点间距为 \(h/2\), 因而若 \(f\in C^4[a,b]\),

\[\begin{aligned} I-S_n&=-\frac{b-a}{2880}h^4f^{(4)}(\xi_n),\\ I-S_{2n}&=-\frac{b-a}{2880}\left(\frac{h}{2}\right)^4f^{(4)}(\xi_{2n}). \end{aligned}\]

继续外推, 定义

\[C_n\equiv\frac{16S_{2n}-S_n}{15},\qquad I=C_n+O(h^6).\]

\(C_n\) 恰为 \(n\) 个宽度为 \(h\)、节点间距为 \(h/4\) 的 Boole 面板组成的复合公式. 若 \(f\in C^6[a,b]\),

\[\begin{aligned} I-C_n&=-\frac{2(b-a)}{945}\left(\frac{h}{4}\right)^6f^{(6)}(\xi_n),\\ I-C_{2n}&=-\frac{2(b-a)}{945}\left(\frac{h}{8}\right)^6f^{(6)}(\xi_{2n}). \end{aligned}\]

在具有相应高阶渐近展开时, 再消去六阶项得

\[I=\frac{64C_{2n}-C_n}{63}+O(h^8).\]

各阶余项估计所需的光滑性不能从低阶公式自动推出; 例如 \(f\in C^{2m+2}[a,b]\) 足以支持梯形余项展开到 \(h^{2m}\), 剩余为 \(O(h^{2m+2})\).3

4.8 Gaussian Quadrature

Definition 取 \(n+1\) 个合适的点 \(x_{k}, k=0,1,2, \cdots, n\), 使得积分公式

\[\int_a^b f(x)\,\mathrm{d}x\approx\sum_{k=0}^n a_k f(x_k)\]

对所有不超过 \(2n+1\) 次的多项式精确成立, 就称该求积公式具有 \(2n+1\) 次代数精度. 此时称为 Gaussian quadrature, 对应节点称为 Gauss points; 对一般函数仍有求积余项.

  • 然后确定权重 \(a_k\). 令 Lagrange 基函数为 \[\ell_k(x)=\prod_{\substack{i=0\\i\neq k}}^n\frac{x-x_i}{x_k-x_i}\] 任意不超过 \(n\) 次的多项式 \(\Pi_n(x)\) 都可以精确写成 \(\Pi_n(x)=\sum_{k=0}^n \Pi_n(x_k)\ell_k(x)\), 对两边积分得 \[\int_a^b \Pi_n(x)\,\mathrm{d}x = \sum_{k=0}^n \Pi_n(x_k)\int_a^b \ell_k(x)\,\mathrm{d}x\] 与上述积分公式比较, 得到 \[a_k=\int_a^b\ell_k(x)\,\mathrm{d}x\]

Theorem 节点 \(x_{k}, k=0,1,2, \cdots, n\) 是 Gauss points 的充分必要条件是, 以这些点为零点的多项式

\[\omega_{n+1}(x)=\prod_{k=0}^n(x-x_k)\]

与任意不超过 \(n\) 次的多项式 \(Q_{n}(x)\) 均正交, 即

\[\int_a^b Q_n(x)\omega_{n+1}(x)\,\mathrm{d}x=0\]

Proof 必要性: 取 \(f=Q_n\omega_{n+1}\), 次数至多 \(2n+1\), 各节点函数值为零, 精确求积即给出正交性. 充分性: 任意次数至多 \(2n+1\) 的多项式均可作带余除法 \(f=q\omega_{n+1}+r\), 其中 \(\deg q,\deg r\le n\). 第一项的积分与节点求和值均为零, 第二项由插值权重精确积分, 故整个公式精确. 对正权重的正交多项式, 节点位于开区间内且互异; 此处权重为 1.

根据这个判据, 构造 Gauss points 的关键是找到满足上述正交性的 \(\omega_{n+1}(x)\). 在 \([-1,1]\) 且权重为 1 的情形下, 可以用 Legendre 多项式来完成这个构造.

Proposition Gauss-Legendre 求积 对于 \([-1,1]\) 上的 Gaussian quadrature, 可以取 \(\omega_{n+1}(x)\) 为 \(n+1\) 次 Legendre 多项式 \(P_{n+1}(x)\) 的常数倍. 此时 \(P_{n+1}(x)\) 的零点就是 Gauss-Legendre quadrature 的 Gauss points.

Proof Legendre 多项式在 \([-1,1]\) 上满足正交关系:

\[\int_{-1}^1 P_i(x)P_j(x)\,\mathrm{d}x=\frac{2}{2i+1}\delta_{ij}\]

因此可以取

\[\omega_{n+1}(x) =\frac{2^{n+1}[(n+1)!]^2}{(2n+2)!}P_{n+1}(x) =\frac{(n+1)!}{(2n+2)!}\frac{\mathrm{d}^{n+1}}{\mathrm{d}x^{n+1}}(x^2-1)^{n+1}\]

对于不超过 \(n\) 次的多项式 \(Q_{n}(x)\), 可以用 Legendre 多项式展开:

\[Q_n(x)=\sum_{k=0}^n\lambda_k P_k(x) \implies \int_{-1}^1\omega_{n+1}(x)Q_n(x)\,\mathrm{d}x =C(n)\sum_{k=0}^n\lambda_k\int_{-1}^1P_k(x)P_{n+1}(x)\,\mathrm{d}x=0.\]

所以 \(P_{n+1}(x)\) 的零点满足 Gauss points 的正交判据.3

  • 此时, 权重可以给出另解, 即求解线性方程组 \[\sum_{k=0}^{n}a_kP_i(x_k) = \int_{-1}^{1}P_i(x)\,\mathrm{d}x\] 注意此处 \(i\) 的编号从 \(0\) 到 \(n\) 而不是 \(n+1\).

Remark 通过区间变换, 可把其他区间上的积分化为 \([-1,1]\) 上的, 从而可以用Gauss-Legendre quadrature

  • 若积分区间为 \(t \in [a,b]\), 令

    \[t=\frac{b-a}{2}x+\frac{b+a}{2}\]

    积分区间化为 \(x \in [-1,1]\), 即

    \[\int_a^b f(t)\,\mathrm{d}t=\frac{b-a}{2}\int_{-1}^1 f\left(\frac{b-a}{2}x+\frac{b+a}{2}\right)\,\mathrm{d}x\]

  • 若为无穷区间积分 \(\int_0^\infty f(t)\,\mathrm{d}t\), 令

    \[t=\tan\left[\frac{\pi}{4}(1+x)\right]\]

    积分区间化为 \(x \in [-1,1]\), 即

    \[\begin{aligned} \int_0^\infty f(t)\,\mathrm{d}t &=\int_{-1}^1 t'(x)f(t(x))\,\mathrm{d}x\\ &\approx\sum_{k=0}^n a_k t'(x_k)f(t(x_k)),\\ t'(x)&=\frac{\pi}{4}\sec^2\left[\frac{\pi}{4}(1+x)\right]. \end{aligned}\]

变换后的积分应收敛; \(t(x)\) 在 \(x=1\) 发散, 变换后的被积函数未必满足有限区间 Gauss 余项所需的光滑性.

Example 计算积分

\[\int_{-1}^1x^2\,\mathrm{d}x\]

Solution 考虑两个 Gauss points 的情形, 取二阶 Legendre 多项式

\[P_2(x)=\frac{1}{2}(3x^2-1)\]

Gauss points 为 \(x_0=-\frac{1}{\sqrt{3}},x_1=\frac{1}{\sqrt{3}}\), 对应的 Lagrange 基函数为

\[\begin{aligned} \ell_0(x) &=\frac{x-x_1}{x_0-x_1} =\frac{x-\frac{1}{\sqrt{3}}}{-\frac{2}{\sqrt{3}}} =\frac{1-\sqrt{3}x}{2},\\ \ell_1(x) &=\frac{x-x_0}{x_1-x_0} =\frac{x+\frac{1}{\sqrt{3}}}{\frac{2}{\sqrt{3}}} =\frac{1+\sqrt{3}x}{2}. \end{aligned}\]

由权重公式 \(a_k=\int_{-1}^{1}\ell_k(x)\,\mathrm{d}x\) 得

\[a_0 =\int_{-1}^{1}\frac{1-\sqrt{3}x}{2}\,\mathrm{d}x=1, \qquad a_1 =\int_{-1}^{1}\frac{1+\sqrt{3}x}{2}\,\mathrm{d}x=1.\]

代入积分公式

\[\int_{-1}^1x^2\,\mathrm{d}x=a_0x_0^2+a_1x_1^2=\frac{1}{3}+\frac{1}{3}=\frac{2}{3}\]

Remark 上述 \(a_k\) 也可直接由矩条件得到. 对 2 点 Gauss-Legendre 公式, 要求对 \(1,x\) 精确:

\[\begin{aligned} 2 &= a_0 + a_1,\\ 0 &= a_0x_0 + a_1x_1 = -\frac{a_0}{\sqrt3} + \frac{a_1}{\sqrt3} \end{aligned}\]

解得 \(a_0 = a_1 = 1\), 与 Lagrange 基积分结果一致. 在此权重下, 公式自动对 \(x^2\) 和 \(x^3\) 也精确 — 这正是 \(N\) 个节点的 Gauss 求积具有 \(2N-1\) 次代数精度(本节 \(N=n+1\)): 2 个节点精确到 3 次多项式, 多出的精度由节点位置的 Legendre 根性质保证.

4.9 Weighted Gaussian Quadrature *

很多情况下我们需要处理带权函数的 Gaussian quadrature:

\[\int_a^b \rho(x)f(x)\,\mathrm{d}x\]

以下假设:

  • \(\rho\geq0\);
  • 所需多项式矩有限;
  • 对所需次数的任意非零多项式 \(q\) 有 \(\int_a^b\rho(x)q(x)^2\,\mathrm{d}x>0\).

这些条件保证相应正交多项式存在.

此时 Gauss points 成立的条件为: 多项式

\[\omega_{n+1}(x)=\prod_{k=0}^n(x-x_k)\]

为区间 \([a,b]\) 上关于权函数 \(\rho(x)\) 的正交多项式:

\[\int_a^b \rho(x)\omega_{n+1}(x)Q_n(x)\,\mathrm{d}x=0\]

各公式使用 \(n+1\) 个节点, 取对应 \(n+1\) 次正交多项式的零点; 对次数至多 \(2n+1\) 的多项式精确, 对一般函数是近似求积. 下面列举常见公式:

  • Gauss-Laguerre quadrature formula

    \(\mathcal{L}_{n}(x)\) 为 Laguerre polynomial:

    \[\begin{aligned} \int_0^\infty \mathrm{e}^{-x}f(x)\,\mathrm{d}x&\approx\sum_{k=0}^n a_k f(x_k),\\ a_k&=-\frac{1}{(n+1)\mathcal{L}'_{n+1}(x_k)\mathcal{L}_n(x_k)}. \end{aligned}\]\[\int_0^\infty \mathrm{e}^{-x}\mathcal{L}_n(x)^2\,\mathrm{d}x=1\]

    递推关系为

    \[\begin{aligned} (n+1)\mathcal{L}_{n+1}(x)&=(2n+1-x)\mathcal{L}_n(x)-n\mathcal{L}_{n-1}(x),\\ \mathcal{L}_0(x)&=1,\\ \mathcal{L}_1(x)&=1-x. \end{aligned}\]

  • Gauss-Hermite quadrature formula

    \[\begin{aligned} \int_{-\infty}^{\infty}\mathrm{e}^{-x^2}f(x)\,\mathrm{d}x&\approx\sum_{k=0}^n a_k f(x_k),\\ a_k&=\frac{2^{n+1}n!\sqrt{\pi}}{H'_{n+1}(x_k)H_n(x_k)}. \end{aligned}\]

    \(H_{n}\left(x\right)\) Hermite polynomial:

    \[\int_{-\infty}^{\infty}\mathrm{e}^{-x^2}H_n(x)^2\,\mathrm{d}x=2^n n!\sqrt{\pi}\]

    递推关系为

    \[\begin{aligned} H_{n+1}(x)&=2xH_n(x)-2nH_{n-1}(x),\\ H_0(x)&=1,\\ H_1(x)&=2x. \end{aligned}\]

  • Gauss-Chebyshev quadrature formula

    \[\begin{aligned} \int_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}\,\mathrm{d}x&\approx\sum_{k=0}^n a_k f(x_k),\\ a_k&=\frac{\pi}{n+1},\qquad x_k=\cos\frac{(2k+1)\pi}{2(n+1)},\quad k=0,\cdots,n. \end{aligned}\]

5 Nonlinear Equations

以下讨论复杂方程的求根方法.

5.1 Bisection Method

二分法 利用介值定理. 假设 \(f\) 在 \([a,b]\) 上连续且 \(f(a)f(b)\lt 0\); 若端点已经为根则直接返回. 取区间中点 \(p_1=\frac{a+b}{2}\)

  1. 如果 \(f(p_{1})\) 与 \(f(a)\) 同号, 说明根在右半区间 \((p_{1}, b)\) 内, 更新 \((a_{2}, b_{2}) = (p_{1}, b)\);

  2. 如果 \(f(p_{1})\) 与 \(f(b)\) 同号, 说明根在左半区间 \((a, p_{1})\) 内, 更新 \((a_{2}, b_{2}) = (a, p_{1})\).

若中点恰为根则终止; 通常按区间半宽小于指定误差容限终止, 不应等待浮点函数值精确为零. 每次保留异号区间, 其长度减半, 中点到该区间中至少一个根的距离不超过半区间长.

二分法具有全局线形收敛性, 对于振荡, 不光滑函数具有较好效果.

5.2 Fixed-Point Iteration

不动点法 对待求根函数 \(f(x)\), 寻找一个不动点为其根的函数, 即

\[f(p)=0\implies p=g(p).\]

以化为寻找不动点的问题. 选取初始值 \(p_{0}\), 迭代求解:

\[p_n=g(p_{n-1})\]

若迭代收敛且 \(g\) 在极限点连续, 则极限是它的不动点. 必须在所选定义域保证不动点与原方程的根等价; 一般的 \(g\) 并不自动收敛:

\[p=\lim_{n\to\infty}p_n=\lim_{n\to\infty}g(p_n)=g(p)\]

Remark 对于同一个 \(f(x)\), 有很多不同的 \(g(x)\) 选取方法. \(g(x)\) 的选取和迭代的收敛速度有很大关系.

Theorem 压缩映射定理: 若

  1. \(g\) 将区间 \([a,b]\) 映射到自身
  2. \(g(x)\) 在 \([a,b]\) 上可导
  3. 存在常数 \(0\leq k\lt1\) 使得 \(|g'(x)| \leq k, \forall x \in (a,b)\)

则对于任意初始值 \(p_0 \in [a,b]\), 由 \(p_n = g(p_{n-1})\) 生成的序列必定收敛到 \([a,b]\) 内的唯一不动点 \(p\).

Example 对于方程

\[f(x)=x^3+4x^2-10=0\]

Solution 尝试以下五种 \(g(x)\) 的取值方法:

\[\begin{aligned} &g_1(x)=x-f(x), &&g_2(x)=\left(\frac{10}{x}-4x\right)^{1/2}, &g_3(x)=\frac{1}{2}(10-x^3)^{1/2},\\ &g_4(x)=\left(\frac{10}{4+x}\right)^{1/2}, &&g_5(x)=x-\frac{f(x)}{f^{\prime}(x)}. \end{aligned}\]

取 \(p_{0}=1.5\).

\(n\)\(g_1(x)\)\(g_2(x)\)\(g_3(x)\)\(g_4(x)\)\(g_5(x)\)
01.51.51.51.51.5
1-0.8750.81651.2869537681.3483997251.373333333
26.7322.99691.4025408041.3673763721.365262015
3-469.7\((-8.65)^{1/2}\)1.3454583741.3649570151.365230014
4\(1.03\times10^8\)1.3751702531.3652647481.365230013
51.3600941931.365225594
61.3678469681.365230576
71.3638870041.365229942
81.3659167341.365230022
91.3648782171.365230012
101.3654100621.365230014
151.3652236801.365230013
201.365230236

第一个函数会发散, 第二个函数会迭代到负数开方而在实数范围内失效. 从第三个到第五个函数, 收敛速度会不断变快.

同一三次方程的两个不动点映射从相同初值出发的蛛网迭代过程

图 5.1 采用本节三次方程例题的第三、第四种不动点映射, 均从 1.5 出发绘制七次迭代. 从对角线上的当前值竖直到蓝色映射曲线, 再水平返回对角线, 得到下一个值. 两者均交替趋近同一根, 右侧收缩更快; 接近根时的步长会在图上重叠. 红线表示算法步骤, 不表示物理运动轨迹. 由本节例题直接计算绘制.

5.3 Newton Method

5.3.1 Newton Method for Simple Roots

Newton method 是另一种求解方程的方法, 其核心思想为在当前位置用切线近似原函数, 将切线与 \(x\) 轴的交点作为下一个近似根.

对于方程 \(f(p)=0\), 函数 \(f(p)\) 可以 Taylor 展开:

\[f(p)=f(p_0)+(p-p_0)f'(p_0)+\frac{(p-p_0)^2}{2}f''(\xi(p))\]

不考虑二阶微分项:

\[p\approx p_0-\frac{f(p_0)}{f'(p_0)}\equiv p_1\]

迭代:

\[p_n=p_{n-1}-\frac{f(p_{n-1})}{f'(p_{n-1})}\]

Remark 若 \(f\in C^2\) 于单根 \(p\) 的邻域、\(f'(p)\ne0\), 且初值足够接近 \(p\), Newton method 至少二阶局部收敛. 迭代中要求 \(f'(p_{n-1})\ne0\); 初值偏离根较远时可能发散或走出定义域.

Remark Newton method 也可以看作不动点迭代. 对于方程 \(f(x)=0\), 定义 Newton map

\[G(x)=x-\frac{f(x)}{f'(x)}\]

Newton method 的迭代就是 \(p_n=G(p_{n-1})\). 上一节例子中的 \(g_5\) 正是 Newton map.

5.3.2 Newton Method for Multiple Roots

Newton method 求解多重根问题时, 如果采用前面的公式, 收敛速度将变为线性收敛. 为了恢复平方收敛速度, 可以将问题再转换为求单根问题.

例如对于 \(m\geq2\) 重根 \(f(x)=(x-p)^m q(x)\), 假设 \(q(p)\ne0\) 且 \(q\) 在根的邻域足够光滑. 普通 Newton 法的渐近误差因子为 \(1-1/m\); 若 \(m\) 已知, 可用 \(x_{n+1}=x_n-mf(x_n)/f'(x_n)\) 恢复至少二阶局部收敛. 若不使用重数, 定义

\[\mu(x)=\frac{f(x)}{f'(x)}\]

将问题转化为求解 \(\mu(x)\) 的单根, 可以计算出

\[\begin{aligned} \mu(x)&=\frac{(x-p)q(x)}{mq(x)+(x-p)q'(x)},\\ g(x)&=x-\frac{\mu(x)}{\mu'(x)} =x-\frac{f(x)f'(x)}{f'(x)^2-f(x)f''(x)}. \end{aligned}\]

在 \(p\) 处对 \(\mu\) 作连续延拓, 有 \(\mu(p)=0\)、\(\mu'(p)=1/m\), 因而它是单根. 对 \(\mu\) 使用 Newton 法, 迭代公式为

\[p_{n+1}=p_n-\frac{f(p_n)f'(p_n)}{f'(p_n)^2-f(p_n)f''(p_n)}\]

计算时要求分母非零; 达到根或残差、步长容限后应停止, 避免在多重根处计算 \(0/0\).

5.4 Secant Method

割线法可以看作 Newton method 的离散化版本: 它避免直接计算导数, 而是用两点连线的斜率来近似切线斜率.

割线法中斜率近似为

\[f'(p_{n-1}) =\lim_{x\to p_{n-1}}\frac{f(x)-f(p_{n-1})}{x-p_{n-1}} \approx\frac{f(p_{n-1})-f(p_{n-2})}{p_{n-1}-p_{n-2}}\]

若 \(f(p_{n-1})\ne f(p_{n-2})\), 代入原迭代式

\[p_n=p_{n-1}-\frac{f(p_{n-1})(p_{n-1}-p_{n-2})}{f(p_{n-1})-f(p_{n-2})}\]

图 5.2 采用 \(f(x)=x^3+4x^2-10\), 左右两栏使用同一函数.

同一三次函数上 Newton 切线与割线分别给出下一次求根近似的几何过程

图 5.2 比较同一三次函数的求根步骤. 左: Newton 法在初值 1.5 处作切线, 切线与横轴的交点给出下一次近似. 右: 割线法取两个初值 1.1、1.5, 两个函数值的连线与横轴交于新的近似. 每栏只显示一步, 旨在比较所需信息与几何构造, 不用于比较不同初值下的总收敛速度. 由该函数直接计算绘制.

5.5 Aitken Acceleration *

利用 Aitken 加速算法可以加速收敛. 考虑序列 \(p_0, p_1, p_2, \cdots\), 该序列将收敛为 \(p\). 可以构造一个新的序列, 使其收敛速度变快. 假设对于足够大的 \(n\), 有

\[\begin{aligned} \frac{p_{n+1}-p}{p_n-p} &\approx\frac{p_{n+2}-p}{p_{n+1}-p},\\ p &\approx\frac{p_{n+2}p_n-p_{n+1}^2}{p_{n+2}-2p_{n+1}+p_n}\\ &=p_n-\frac{(p_{n+1}-p_n)^2}{p_{n+2}-2p_{n+1}+p_n}. \end{aligned}\]

当 \(\Delta^2p_n\ne0\) 时可构造新的序列 \(\hat{p}_n\). 对渐近线性误差且极限比值 \(\lambda\ne1\) 的序列, 它可消去主导的几何误差项; 并非任意收敛序列都必然加速:

\[\hat{p}_n=p_n-\frac{(p_{n+1}-p_n)^2}{p_{n+2}-2p_{n+1}+p_n}\]

如果定义

\[\begin{aligned} \Delta p_n&=p_{n+1}-p_n,\\ \Delta^2p_n&=p_{n+2}-2p_{n+1}+p_n. \end{aligned}\]

可以写成差分形式:

\[\hat{p}_n=p_n-\frac{(\Delta p_n)^2}{\Delta^2p_n}\]

5.6 Steffensen Method *

Steffensen 方法把 Aitken 加速和不动点迭代结合起来, 每轮都对两步迭代作加速. 若 \(g\) 在不动点 \(p\) 邻域足够光滑、\(g'(p)\ne1\), 初值足够接近且下面的分母非零, 通常可获得至少二阶局部收敛; 它不需要显式导数. 若当前点已为不动点则先终止.

从当前近似 \(p_{0}^{(k)}\) 出发, 先做两步不动点迭代:

\[\begin{aligned} p_1^{(k)}&=g(p_0^{(k)}),\\ p_2^{(k)}&=g(p_1^{(k)}). \end{aligned}\]

这三个点利用 Aitken 公式, 计算出加速度的值作为下一轮的起点:

\[p_0^{(k+1)}=\hat{p}_0^{(k)} =p_0^{(k)}-\frac{(p_1^{(k)}-p_0^{(k)})^2}{p_2^{(k)}-2p_1^{(k)}+p_0^{(k)}}\]

再重复上述步骤即可.

5.7 Qin Jiushao Algorithm *

考虑代数多项式:

\[P(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\]

可以将其写作如下形式:

\[P(x)=(\cdots((a_nx+a_{n-1})x+a_{n-2})x+\cdots+a_1)x+a_0\]

令

\[\begin{aligned} b_n&=a_n,\\ b_k&=a_k+b_{k+1}x_0,\\ b_0&=P(x_0). \end{aligned}\]

考虑多项式:

\[Q(x)=b_nx^{n-1}+b_{n-1}x^{n-2}+\cdots+b_2x+b_1\]

这个多项式是 \(P(x)\) 除以 \((x-x_0)\) 的商, 余数为 \(b_0\), 即

\[P(x)=(x-x_0)Q(x)+b_0\]

求导:

\[\begin{aligned} P'(x)&=Q(x)+(x-x_0)Q'(x),\\ P'(x_0)&=Q(x_0). \end{aligned}\]

再利用 Newton method:

\[x_1=x_0-\frac{P(x_0)}{P'(x_0)}=x_0-\frac{P(x_0)}{Q(x_0)}\]

收敛到根 \(t\) 时, \(P(t)=0\), 此时对应的 \(b_{0}=0\). 于是

\[P(x)=(x-t)Q(x)\]

\(Q(x)\) 在精确代数下保留除去一次重数后的其他根. 可以继续对 \(Q(x)\) 求根以逐次降阶, 但 Newton 法仍需要合适初值, 也可能遇到复根; 用近似根作降阶会传播舍入误差, 宜回到原多项式检查最终根的残差.

5.8 Brent Method *

Brent 算法是一种混合型求根方法. 假设 \(f\) 在初始区间连续, 两端函数值异号(端点为根则先返回). 它保持括根区间, 在割线、逆二次插值和二分法之间选择, 用步长保护获得可靠的区间收缩.

下面采用 \(a,b\) 为括根端点的记号, 始终保持 \(f(a)f(b)\lt0\)、\(|f(b)|\leq|f(a)|\), 而 \(c\) 保留上一轮的 \(b\). 三个点的函数值互异时, 可通过逆二次插值构造候选点(不同实现对 \(a,b,c\) 的命名可能不同):

\[\begin{aligned} x &=\frac{[y-f(a)][y-f(b)]c}{[f(c)-f(a)][f(c)-f(b)]}\\ &\quad+\frac{[y-f(b)][y-f(c)]a}{[f(a)-f(b)][f(a)-f(c)]}\\ &\quad+\frac{[y-f(c)][y-f(a)]b}{[f(b)-f(c)][f(b)-f(a)]}. \end{aligned}\]

\(y=0\) 时可以求解出对应的 \(x\) 值, 解得

\[\begin{aligned} x&=b+\frac{P}{Q},\\ R&=\frac{f(b)}{f(c)},\qquad S=\frac{f(b)}{f(a)},\qquad T=\frac{f(a)}{f(c)},\\ P&=S[T(R-T)(c-b)-(1-R)(b-a)],\\ Q&=(T-1)(R-1)(S-1). \end{aligned}\]

下面要考虑这个计算出的结果 \(x\) 要不要接受, 一般来说, 接受新值的具体条件:

  • 候选点位于当前括根区间 \((\min(a,b),\max(a,b))\) 内; 端点重排后并不保证 \(a\lt b\);

  • 与上次迭代相比, 区间收缩不会太慢.

若三个函数值不互异, 优先退化为割线候选 \(x=b-f(b)(b-a)/(f(b)-f(a))\); 若候选不安全或分母过小, 使用二分法. 上述区间与步长条件只是保护思路, 完整实现还必须保存前两步长度、检查插值步是否足够缩短以及绝对/相对容限, 不能仅凭区间内条件保证 Brent 的收敛界. 可参照 Netlib 的成熟实现.4 放弃插值时使用二分法:

\[x=\frac{a+b}{2}\]

随后再迭代, 更新三个点: 计算

\[y=f(x)\]

若 \(y=0\) 则返回 \(x\); 否则按异号性更新:

  1. 如果 \(f(b)\) 与 \(y\) 同号, 新括根端点为 \((a',b')=(a,x)\).
  2. 否则新括根端点为 \((a',b')=(x,b)\).

这里括根区间始终按两端点的最小值与最大值理解; 达到区间宽度容限时也应终止.

如果 \(|f(a')| \lt |f(b')|\), 则需要更换 \(a',b'\), 因为要保持 \(f(b')\) 的绝对值最小. 最后再取 \(c' = b\), 重复上面的过程.

5.9 Newton Method for Nonlinear Systems

考虑非线性方程组 \(\mathbf{F}(\mathbf{x})=\mathbf{0}\) 的求解. 每一步要求当前 Jacobian 可逆; 若解处 Jacobian 可逆、\(\mathbf{F}\) 在其邻域二阶连续可微且初值足够接近, Newton 法至少二阶局部收敛, 并不保证任意初值都成功:

  1. 构造 Jacobian matrix:

    \[\mathbf{J}(\mathbf{x})= \left(\frac{\partial f_i}{\partial x_j}\right)\]

  2. 令

    \[G(\mathbf{x})=\mathbf{x}-\mathbf{J}(\mathbf{x})^{-1}\mathbf{F}(\mathbf{x})\]

    迭代公式:

    \[\mathbf{x}^{(k)} =G(\mathbf{x}^{(k-1)})=\mathbf{x}^{(k-1)}-\mathbf{J}(\mathbf{x}^{(k-1)})^{-1}\mathbf{F}(\mathbf{x}^{(k-1)}).\]

  3. 在实际的求解过程中, 一般不会直接求解 Jacobian matrix 的逆矩阵. 定义 \(\mathbf{y}=\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\), 上式化为

    \[\mathbf{J}(\mathbf{x}^{(k-1)})\mathbf{y}=-\mathbf{F}(\mathbf{x}^{(k-1)}).\]

    求解该线性方程组即可.

Example 求解非线性方程组

\[\begin{aligned} f_1(x_1,x_2,x_3)&=3x_1-\cos(x_2x_3)-\frac{1}{2},\\ f_2(x_1,x_2,x_3)&=x_1^2-81(x_2+0.1)^2+\sin x_3+1.06,\\ f_3(x_1,x_2,x_3)&=\mathrm{e}^{-x_1x_2}+20x_3+\frac{10\pi-3}{3}. \end{aligned}\]

Solution Jacobian matrix 为:

\[\mathbf{J}(x_1,x_2,x_3)= \begin{pmatrix} 3&x_3\sin(x_2x_3)&x_2\sin(x_2x_3)\\ 2x_1&-162(x_2+0.1)&\cos x_3\\ -x_2\mathrm{e}^{-x_1x_2}&-x_1\mathrm{e}^{-x_1x_2}&20 \end{pmatrix}\]

取定初始值 \(\mathbf{x}^{(0)}=(0.1,0.1,-0.1)^{\mathsf{T}}\), 计算得

\[\mathbf{F}(\mathbf{x}^{(0)})=(-1.19995000042,-2.26983341665,8.46202534572)^{\mathsf{T}}\]

求解线性方程组:

\[\mathbf{J}(\mathbf{x}^{(0)})\mathbf{y}^{(0)}=-\mathbf{F}(\mathbf{x}^{(0)})\]

求解得

\[\mathbf{y}^{(0)} =\begin{pmatrix} 0.39986967293\\ -0.08053315146\\ -0.42152047194 \end{pmatrix}, \qquad \mathbf{x}^{(1)} =\mathbf{x}^{(0)}+\mathbf{y}^{(0)} =\begin{pmatrix} 0.49986967293\\ 0.01946684854\\ -0.52152047194 \end{pmatrix}.\]

继续迭代并检验原方程残差:

\(k\)\(x_1^{(k)}\)\(x_2^{(k)}\)\(x_3^{(k)}\)\(\Vert \mathbf{F}(\mathbf{x}^{(k)})\Vert _\infty\)
20.50001424020.0015885914-0.52355696430.0258891
30.50000011350.0000124448-0.52359845010.000201223
40.5000000000\(7.7579\times10^{-10}\)-0.5235987756\(1.2543\times10^{-8}\)
50.5000000000\(3.36\times10^{-18}\)-0.5235987756\(1.78\times10^{-15}\)

表中残差用未舍入的迭代值计算. 对应精确解为

\[\boxed{\mathbf x_*=(1/2,0,-\pi/6)^{\mathsf T}},\qquad\mathbf F(\mathbf x_*)=\mathbf0,\]

可直接代入三式验证. 这是该初值收敛到的根, 未据此声称全部根已枚举.

5.10 Quasi-Newton Method: Broyden Method *

如果 Jacobian matrix 很难计算, 则需要用近似矩阵 \(\mathbf{A}\) 来代替. 近似矩阵 \(\mathbf{A}\) 需要满足:

  • 割线条件

    对于标量函数 \(f(x)\), 割线近似为:

    \[f'(x_1)\approx\frac{f(x_1)-f(x_0)}{x_1-x_0}\]

    对于向量函数, 自然推广为

    \[\mathbf{A}_1(\mathbf{x}^{(1)}-\mathbf{x}^{(0)}) =\mathbf{F}(\mathbf{x}^{(1)})-\mathbf{F}(\mathbf{x}^{(0)})\]

    记

    \[\begin{aligned} \mathbf{s}_1&=\mathbf{x}^{(1)}-\mathbf{x}^{(0)},\\ \mathbf{y}_1&=\mathbf{F}(\mathbf{x}^{(1)})-\mathbf{F}(\mathbf{x}^{(0)}). \end{aligned}\]

    割线条件可以写作

    \[\mathbf{A}_1\mathbf{s}_1=\mathbf{y}_1\]

  • 正交条件

    对于任意与 \(\mathbf{s}_1\) 正交的向量 \(\mathbf{z}\).

    \[\begin{aligned} \mathbf{s}_1^{\mathsf{T}}\mathbf{z}&=0,\\ \mathbf{A}_1\mathbf{z}&=\mathbf{A}_0\mathbf{z}. \end{aligned}\]

\(\mathbf{A}_1\) 需要满足秩-1 更新:

\[\mathbf{A}_1=\mathbf{A}_0+\frac{\mathbf{y}_1-\mathbf{A}_0\mathbf{s}_1}{\lVert\mathbf{s}_1\rVert^2}\mathbf{s}_1^{\mathsf{T}}\]

验证正交条件:

\[\mathbf{A}_1\mathbf{z} =\mathbf{A}_0\mathbf{z} +\frac{\mathbf{y}_1-\mathbf{A}_0\mathbf{s}_1}{\lVert\mathbf{s}_1\rVert^2}(\mathbf{s}_1^{\mathsf{T}}\mathbf{z}) =\mathbf{A}_0\mathbf{z}\]

于是 \(\mathbf{A}_{i}\) 的迭代公式:

\[\mathbf{A}_i=\mathbf{A}_{i-1}+\frac{\mathbf{y}_i-\mathbf{A}_{i-1}\mathbf{s}_i}{\lVert\mathbf{s}_i\rVert^2}\mathbf{s}_i^{\mathsf{T}}\]

\(\mathbf{x}_{i}\) 的迭代公式:

\[\mathbf{x}^{(i+1)}=\mathbf{x}^{(i)}-\mathbf{A}_i^{-1}\mathbf{F}(\mathbf{x}^{(i)})\]

\(\mathbf{x}_{i}\) 的迭代中需要求解 \(\mathbf{A}_{i}\) 的逆矩阵, 根据 Sherman-Morrison 公式:

\[\begin{aligned} (\mathbf{A}+\mathbf{u}\mathbf{v}^{\mathsf{T}})^{-1} &=\mathbf{A}^{-1}-\frac{\mathbf{A}^{-1}\mathbf{u}\mathbf{v}^{\mathsf{T}}\mathbf{A}^{-1}}{1+\mathbf{v}^{\mathsf{T}}\mathbf{A}^{-1}\mathbf{u}},\\ \mathbf{A}_i^{-1} &=\left(\mathbf{A}_{i-1}+\frac{\mathbf{y}_i-\mathbf{A}_{i-1}\mathbf{s}_i}{\lVert\mathbf{s}_i\rVert^2}\mathbf{s}_i^{\mathsf{T}}\right)^{-1}\\ &=\mathbf{A}_{i-1}^{-1} +\frac{(\mathbf{s}_i-\mathbf{A}_{i-1}^{-1}\mathbf{y}_i)\mathbf{s}_i^{\mathsf{T}}\mathbf{A}_{i-1}^{-1}}{\mathbf{s}_i^{\mathsf{T}}\mathbf{A}_{i-1}^{-1}\mathbf{y}_i}. \end{aligned}\]

上述公式用于更新近似逆矩阵, 可节省重复求逆的计算量, 本身不改变算法的收敛阶. 它要求 \(\mathbf{A}_{i-1}\) 可逆且 \(\mathbf{s}_i^{\mathsf{T}}\mathbf{A}_{i-1}^{-1}\mathbf{y}_i\ne0\); 秩-1 更新还要求 \(\mathbf{s}_i\ne0\). 若分母过小应重建近似矩阵或跳过更新. 也可直接求解 \(\mathbf{A}_i\Delta\mathbf{x}=-\mathbf{F}\) 而不显式求逆.

Example 求解非线性方程组

\[\mathbf{g}(x,y,z)= \begin{pmatrix} x^2+y^2+z^2-3\\ x^2+y^2-z-1\\ x+y+z-3 \end{pmatrix}\]

初始值:

\[\mathbf{x}^{(0)}=(x^{(0)},y^{(0)},z^{(0)})^{\mathsf{T}}=(1,0,1)^{\mathsf{T}}\]

Solution 取 \(\mathbf{A}_0=\mathbf{J}(\mathbf{x}^{(0)})\). Jacobian matrix 为:

\[\mathbf{J}(x,y,z)= \begin{pmatrix} 2x&2y&2z\\ 2x&2y&-1\\ 1&1&1 \end{pmatrix}\]\[\begin{aligned} \mathbf{J}(\mathbf{x}^{(0)})(\mathbf{x}^{(1)}-\mathbf{x}^{(0)})&=-\mathbf{g}(\mathbf{x}^{(0)}),\\ \mathbf{x}^{(1)}&=\left(\frac{3}{2},\frac{1}{2},1\right)^{\mathsf{T}}. \end{aligned}\]\[\mathbf{A}_1=\mathbf{A}_0+\frac{\mathbf{y}_1-\mathbf{A}_0\mathbf{s}_1}{\lVert\mathbf{s}_1\rVert^2}\mathbf{s}_1^{\mathsf{T}} = \begin{pmatrix} \frac{5}{2}&\frac{1}{2}&2\\ \frac{5}{2}&\frac{1}{2}&-1\\ 1&1&1 \end{pmatrix}\]

再迭代:

\[\begin{aligned} \mathbf{A}_1(\mathbf{x}^{(2)}-\mathbf{x}^{(1)})&=-\mathbf{g}(\mathbf{x}^{(1)}),\\ \mathbf{x}^{(2)}&=\left(\frac{5}{4},\frac{3}{4},1\right)^{\mathsf{T}}. \end{aligned}\]

这两步是迭代示范, 尚未达到最终解. 从前两式相减得 \(z^2+z-2=0\), 故 \(z=1\) 或 \(-2\); 后者使 \(x^2+y^2=-1\), 无实解. 当 \(z=1\) 时, \(x+y=2\)、\(x^2+y^2=2\), 从而

\[\boxed{(x,y,z)=(1,1,1)}\]

是唯一实根. 该根处 Jacobian 的前两列相同, 矩阵奇异; 因而不能直接将要求根处 Jacobian 非奇异的局部超线性收敛结论套用于本例.

6 Optimization Problems *

优化问题是在指定约束或定义域内寻找目标函数的极小或极大值. 本章的梯度法与 Nelder-Mead 法主要用于局部无约束极小化, 一般不保证全局最优; 线性规划在可行且有界时可求全局最优.

6.1 Fletcher-Reeves Algorithm *

FR(Fletcher-Reeves) 算法是非线性共轭梯度法, 常用于大规模无约束优化. 它只使用函数值及一阶梯度, 不需要显式计算或存储 Hessian matrix. 对正定二次目标配合精确线搜索时, 它与线性共轭梯度法一致, 搜索方向关于固定 Hessian 共轭; 对一般非线性函数不能直接断言各方向严格共轭.

对于一般可微函数 \(f(x)\), 求解其极小值.

初始化:

\[p_0=-\nabla f(x_0)\]

现在要沿着 \(p_k\) 方向寻找合适的步长 \(t_k\), 使得 \(f(x)\) 下降最多. 迭代公式:

\[x_{k+1}=x_k+t_kp_k\]

定义

\[\beta_{k+1}^{\mathrm{FR}}=\frac{\lVert\nabla f(x_{k+1})\rVert^2}{\lVert\nabla f(x_k)\rVert^2}\]

新方向:

\[p_{k+1}=-\nabla f(x_{k+1})+\beta_{k+1}^{\mathrm{FR}}p_k\]

\(t_k>0\) 一般采用一维线搜索来确定, 不必求整条射线上的精确最小值. 为保持 FR 的下降方向, 常用强 Wolfe 条件:

\[\begin{aligned} f(x_k+t_kp_k)&\leq f(x_k)+c_1t_k\nabla f(x_k)^{\mathsf{T}}p_k,\\ |\nabla f(x_k+t_kp_k)^{\mathsf{T}}p_k|&\leq c_2|\nabla f(x_k)^{\mathsf{T}}p_k|,\qquad 0\lt c_1\lt c_2\lt\tfrac12. \end{aligned}\]

若方向不再满足 \(\nabla f(x_k)^{\mathsf{T}}p_k\lt0\), 可重启为 \(p_k=-\nabla f(x_k)\). 梯度已经为零时应停止, 避免计算 \(\beta\) 的零分母.

继续迭代条件:

\[\lVert\nabla f(x_k)\rVert>\varepsilon\]

6.2 Linear Conjugate Gradient Method *

线性共轭梯度法适用于实对称正定矩阵 \(\mathbf{A}\) 所定义的二次函数:

\[f(\mathbf{x})=\frac{1}{2}\mathbf{x}^{\mathsf{T}}\mathbf{A}\mathbf{x}-\mathbf{b}^{\mathsf{T}}\mathbf{x}+c\]

于是

\[\nabla f=0\implies\mathbf{A}\mathbf{x}-\mathbf{b}=0\]

问题转化为求解线性方程组.

初始化:

\[\begin{aligned} \mathbf{r}_0&=\mathbf{A}\mathbf{x}_0-\mathbf{b},\\ \mathbf{p}_0&=-\mathbf{r}_0. \end{aligned}\]

沿当前方向作精确线搜索. 由

\[\frac{\mathrm{d}}{\mathrm{d}t}f(\mathbf{x}_k+t\mathbf{p}_k) =\langle\mathbf{r}_k,\mathbf{p}_k\rangle+t\langle\mathbf{p}_k,\mathbf{A}\mathbf{p}_k\rangle=0,\]

得到步长

\[t_k=-\frac{\langle\mathbf{r}_k,\mathbf{p}_k\rangle}{\langle\mathbf{p}_k,\mathbf{A}\mathbf{p}_k\rangle}.\]

初始方向为 \(-\mathbf{r}_0\), 所以非零梯度下 \(t_0>0\).

按 \(\mathbf{x}_{k+1}=\mathbf{x}_k+t_k\mathbf{p}_k\)、\(\mathbf{r}_{k+1}=\mathbf{r}_k+t_k\mathbf{A}\mathbf{p}_k\) 的顺序更新位置和梯度残差, 检查残差容限后再计算方向系数, 并令 \(\mathbf{p}_{k+1}=-\mathbf{r}_{k+1}+\beta_{k+1}\mathbf{p}_k\). 这里 \(\mathbf{r}_k=\mathbf{A}\mathbf{x}_k-\mathbf{b}\); 若改用方程残差 \(\mathbf{b}-\mathbf{A}\mathbf{x}_k\), 残差及其更新式整体反号, 但实际搜索方向不变.

由共轭条件 \(\langle\mathbf{p}_{k+1},\mathbf{A}\mathbf{p}_k\rangle=0\) 得到相应方向系数. 在精确算术下, 残差互相正交、方向关于 \(\mathbf{A}\) 共轭, 并有

\[t_k=\frac{\lVert\mathbf{r}_k\rVert^2}{\langle\mathbf{p}_k,\mathbf{A}\mathbf{p}_k\rangle},\qquad \beta_{k+1}=\frac{\lVert\mathbf{r}_{k+1}\rVert^2}{\lVert\mathbf{r}_k\rVert^2}.\]

这些步长与方向系数正是 FR 法在正定二次目标和精确线搜索下的形式. \(n\) 维问题至多 \(n\) 步到达解; 浮点计算可能因失去正交性而需要更多步. 达到残差容限后停止, 不应继续计算零方向或零分母.

6.3 Quasi-Newton Optimization *

一维情形可对驻点方程 \(f'(x)=0\) 使用 Newton method. 这不自动判定极小值: 若极限驻点处 \(f''(x)>0\) 则可判定它是严格局部极小值; 迭代还要求分母非零并需局部收敛条件. 公式为:

\[x_{k+1}=x_k-\frac{f'(x_k)}{f''(x_k)}\]

多维情形可以使用 Newton 或 quasi-Newton method. 若 \(f\in C^2\), 在 \(\mathbf{x}_k\) 附近的 Taylor 展开为

\[f(\mathbf{x}_k+\Delta\mathbf{x}_k) =f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)^{\mathsf{T}}\Delta\mathbf{x}_k +\frac{1}{2}\Delta\mathbf{x}_k^{\mathsf{T}}\mathbf{H}(\mathbf{x}_k)\Delta\mathbf{x}_k +o(\lVert\Delta\mathbf{x}_k\rVert^2).\]

其中 \(\mathbf{H}\) 为 Hessian matrix:

\[\mathbf{H}(\mathbf{x})= \begin{pmatrix} \frac{\partial^2 f}{\partial x_1^2}&\frac{\partial^2 f}{\partial x_1\partial x_2}&\cdots&\frac{\partial^2 f}{\partial x_1\partial x_n}\\ \frac{\partial^2 f}{\partial x_2\partial x_1}&\frac{\partial^2 f}{\partial x_2^2}&\cdots&\frac{\partial^2 f}{\partial x_2\partial x_n}\\ \vdots&\vdots&\cdots&\vdots\\ \frac{\partial^2 f}{\partial x_n\partial x_1}&\frac{\partial^2 f}{\partial x_n\partial x_2}&\cdots&\frac{\partial^2 f}{\partial x_n^2} \end{pmatrix}\]

一般我们选用近似矩阵 \(\mathbf{B}\) 来代替 Hessian matrix.

\[\nabla f(\mathbf{x}_k+\Delta\mathbf{x}_k)\approx\nabla f(\mathbf{x}_k)+\mathbf{B}_k\Delta\mathbf{x}_k\]

求解线性方程组 \(\mathbf{B}_k\mathbf{p}_k=-\nabla f(\mathbf{x}_k)\), 并通过线搜索选取 \(t_k>0\), 令 \(\Delta\mathbf{x}_k=t_k\mathbf{p}_k\). 迭代为

\[\mathbf{x}_{k+1}=\mathbf{x}_k+\Delta\mathbf{x}_k.\]

若 \(\mathbf{B}_k\) 对称正定且梯度非零, 则 \(\nabla f(\mathbf{x}_k)^{\mathsf{T}}\mathbf{p}_k\lt0\), 因而 \(\mathbf{p}_k\) 是下降方向. 非凸问题中的不定近似矩阵需要修正或其他步长保护.

但现在近似矩阵 \(\mathbf{B}\) 还是未知的, 需要其他条件去计算得到 \(\mathbf{B}\) 的迭代关系, 不同的算法给出了不同的计算过程.

6.4 Broyden Algorithm *

这里将求非线性方程组的 Broyden 秩-1 更新用于 \(\nabla f=0\), 近似 Jacobian 为 \(\mathbf{B}_k\). 更新满足割线条件 \(\mathbf{B}_{k+1}\mathbf{s}_k=\mathbf{y}_k\), 但通常不保持 \(\mathbf{B}_k\) 的对称性或正定性, 因而不能仅凭它保证极小化的下降方向. 它也不等同于优化中保持对称性的 Broyden 更新族. 更新及逆更新为

\[\begin{aligned} \mathbf{B}_{k+1}&=\mathbf{B}_k+\frac{(\mathbf{y}_k-\mathbf{B}_k\mathbf{s}_k)\mathbf{s}_k^{\mathsf{T}}}{\mathbf{s}_k^{\mathsf{T}}\mathbf{s}_k},\\ \mathbf{B}_{k+1}^{-1}&=\mathbf{B}_k^{-1}+\frac{(\mathbf{s}_k-\mathbf{B}_k^{-1}\mathbf{y}_k)\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k^{-1}}{\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k^{-1}\mathbf{y}_k} \end{aligned}\]

其中

\[\begin{aligned} \mathbf{s}_k&=\Delta(\mathbf{x}_k)=\mathbf{x}_{k+1}-\mathbf{x}_k,\\ \mathbf{y}_k&=\nabla f(\mathbf{x}_{k+1})-\nabla f(\mathbf{x}_k). \end{aligned}\]

上述秩-1 更新要求 \(\mathbf{s}_k\ne0\). 逆更新还要求 \(\mathbf{B}_k\) 可逆且 \(\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k^{-1}\mathbf{y}_k\ne0\); 分母过小时不应直接使用.

6.5 BFGS(Broyden-Fletcher-Goldfarb-Shanno) Algorithm *

BFGS 算法对 Hessian 近似作秩-2 更新. 以下假设 \(\mathbf{B}_k\) 对称正定且 \(\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k>0\), 其中 \(\mathbf{s}_k=\mathbf{x}_{k+1}-\mathbf{x}_k\)、\(\mathbf{y}_k=\nabla f(\mathbf{x}_{k+1})-\nabla f(\mathbf{x}_k)\). 这两个条件保证下式分母非零, 并使更新后的矩阵仍对称正定:

\[\mathbf{B}_{k+1}=\mathbf{B}_k+\alpha\mathbf{u}\mathbf{u}^{\mathsf{T}}+\beta\mathbf{v}\mathbf{v}^{\mathsf{T}}\]

选择

\[\mathbf{u}=\mathbf{y}_k,\qquad\mathbf{v}=\mathbf{B}_k\mathbf{s}_k\]

于是割线条件可以写作

\[\mathbf{B}_{k+1}\mathbf{s}_k=\mathbf{y}_k\]

代入 \(\mathbf{B}_{k+1}\) 得

\[\mathbf{B}_k\mathbf{s}_k +\alpha\mathbf{y}_k(\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k) +\beta(\mathbf{B}_k\mathbf{s}_k)(\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k\mathbf{s}_k) =\mathbf{y}_k\]

可以看到

\[\begin{aligned} \alpha&=\frac{1}{\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k},\\ \beta&=-\frac{1}{\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k\mathbf{s}_k}. \end{aligned}\]

得到 \(\mathbf{B}_{k}\) 的迭代公式:

\[\mathbf{B}_{k+1} =\mathbf{B}_k +\frac{\mathbf{y}_k\mathbf{y}_k^{\mathsf{T}}}{\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k} -\frac{\mathbf{B}_k\mathbf{s}_k\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k^{\mathsf{T}}}{\mathbf{s}_k^{\mathsf{T}}\mathbf{B}_k\mathbf{s}_k}\]

逆矩阵:

\[\mathbf{B}_{k+1}^{-1} =\mathbf{B}_k^{-1} +\frac{(\mathbf{s}_k^{\mathsf{T}}\mathbf{y}_k+\mathbf{y}_k^{\mathsf{T}}\mathbf{B}_k^{-1}\mathbf{y}_k)\mathbf{s}_k\mathbf{s}_k^{\mathsf{T}}}{(\mathbf{s}_k^{\mathsf{T}}\mathbf{y}_k)^2} -\frac{\mathbf{B}_k^{-1}\mathbf{y}_k\mathbf{s}_k^{\mathsf{T}}+\mathbf{s}_k\mathbf{y}_k^{\mathsf{T}}\mathbf{B}_k^{-1}}{\mathbf{s}_k^{\mathsf{T}}\mathbf{y}_k}\]

下降方向配合 Wolfe 线搜索可保证 \(\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k>0\): 由 \(\mathbf{s}_k=t_k\mathbf{p}_k\) 和曲率条件得

\[\mathbf{y}_k^{\mathsf{T}}\mathbf{s}_k =t_k[\nabla f(\mathbf{x}_{k+1})-\nabla f(\mathbf{x}_k)]^{\mathsf{T}}\mathbf{p}_k \geq t_k(c_2-1)\nabla f(\mathbf{x}_k)^{\mathsf{T}}\mathbf{p}_k>0,\qquad0\lt c_2\lt1.\]

若数值计算未满足曲率条件, 应跳过或阻尼该次更新, 不能直接使用零分母或负曲率更新.

6.6 Linear Programming and the Simplex Method *

线性规划问题是在一组线性约束下优化线性目标函数. 标准形式可写为

\[\max z=\mathbf{c}^{\mathsf{T}}\mathbf{x}, \qquad \mathbf{A}\mathbf{x}\leq \mathbf{b}, \qquad \mathbf{x}\geq 0.\]

可行域是凸多面体. 对上式含非负约束的形式, 若可行且目标有有限最优值, 至少有一个最优解位于顶点. 单纯形法从基本可行解出发作换基, 通常移到目标值更好的相邻顶点; 退化换基也可能保持同一顶点及目标值, 需要防循环规则.

下面的单纯形表采用目标行 \(z-\mathbf{c}^{\mathsf{T}}\mathbf{x}=0\) 的符号约定. 若 \(\mathbf{b}\geq0\), 松弛变量构成初始基本可行解; 否则应先用第一阶段法等方法寻找基本可行解, 无可行解时终止. 在已有基本可行解的前提下, 步骤为:

  1. 把不等式约束加松弛变量, 化为等式约束.
  2. 建立初始单纯形表.
  3. 进入变量: 目标行中最负的系数所在列.
  4. 离开变量: 对进入列正系数作最小非负比值 \(b_i/a_{ij}\). 若该列没有正系数且其目标系数为负, 则目标沿该方向无界. 退化时可采用 Bland 等防循环选取规则.
  5. 主元行归一化, 并消去进入列其他元素.
  6. 若目标行所有系数非负, 当前解最优; 否则继续迭代.

Example 求

\[\max f=4x_1+6x_2\]

在约束

\[-x_1+x_2\leq 11, \qquad x_1+x_2\leq 27, \qquad 2x_1+5x_2\leq 90, \qquad x_1,x_2\geq0\]

下的最大值.

Solution 加松弛变量:

\[-x_1+x_2+s_1=11, \qquad x_1+x_2+s_2=27, \qquad 2x_1+5x_2+s_3=90.\]

初始单纯形表为

\(x_1\)\(x_2\)\(s_1\)\(s_2\)\(s_3\)\(b\)
\(-1\)110011
1101027
2500190
\(-4\)\(-6\)0000

进入列为 \(x_2\), 因为目标行中 \(-6\) 最负. 比值为

\[\frac{11}{1}=11, \qquad \frac{27}{1}=27, \qquad \frac{90}{5}=18.\]

最小非负比值在第一行, 主元为 1. 对主元列作行变换后继续迭代. 最终单纯形表可化为

\(x_1\)\(x_2\)\(s_1\)\(s_2\)\(s_3\)\(b\)
010\(-2/3\)\(1/3\)12
001\(7/3\)\(-2/3\)14
100\(5/3\)\(-1/3\)15
000\(8/3\)\(2/3\)132

因此最优解为

\[(x_1,x_2)=(15,12), \qquad f_{\max}=132.\]

6.7 Nelder-Mead Downhill Simplex Method *

下山单纯形法的核心思路十分直观, 可以理解为在 \(N\) 维空间中构造了一个可以伸缩变形的探针来寻找最低点. 算法计算单纯形每个顶点的函数值, 并通过反射, 扩张, 收缩和缩小四种操作, 让单纯形在解空间中不断移动变形, 逐渐向函数值最小的区域移动.

设目标函数 \(f:\mathbb R^n\to\mathbb R\), 取 \(n+1\) 个仿射无关顶点 \(x_1,\cdots,x_{n+1}\) 构成初始单纯形. 本节 \(x_i\in\mathbb R^n\) 表示整个顶点向量, 而非函数的第 \(i\) 个自变量.

首先对顶点函数值排序:

\[f(x_1)\leq f(x_2)\leq\cdots\leq f(x_{n+1})\]

\(f(x_{n+1})\) 是最大的, 因此是最差点. 计算除最差点以外的其余 \(n\) 个点的中心:

\[x_0=\frac{1}{n}\sum_{i=1}^{n}x_i\]

  1. 反射

    \[x_r=x_0+\alpha(x_0-x_{n+1})\]

    若 \(f(x_1) \leq f(x_r) \lt f(x_n)\), 用 \(x_{n+1} \rightarrow x_r\), 重复过程;

    若 \(f(x_r) \lt f(x_1)\), 进行扩张操作; 若 \(f(x_r) \geq f(x_n)\), 进行收缩操作.

  2. 扩张

    \[x_e=x_0+\gamma(x_r-x_0)\]

    若 \(f(x_{e})\lt f(x_{r})\), 则用 \(x_{e}\) 替换 \(x_{n+1}\), 说明扩张成功; 否则扩张失败, 保留反射点, 用 \(x_{r}\) 替换 \(x_{n+1}\), 重复过程.

  3. 收缩

    有两种情况:

    (a) 如果 \(f(x_{r})\lt f(x_{n+1})\), 说明反射点比最差点更好, 进行外部收缩:

    \[x_c=x_0+\rho(x_r-x_0)\]

    若 \(f(x_{c}) \leq f(x_{r})\), 则用 \(x_{c}\) 替换 \(x_{n+1}\), 重复过程; 否则收缩失败, 进行缩小操作.

    (b) 如果 \(f(x_{r}) \geq f(x_{n+1})\), 说明反射点比最差点更差, 进行内部收缩:

    \[x_c=x_0+\rho(x_{n+1}-x_0)\]

    若 \(f(x_{c})\lt f(x_{n+1})\), 则用 \(x_{c}\) 替换 \(x_{n+1}\), 重复过程; 否则收缩失败, 进行缩小操作.

  4. 缩小

    如果收缩失败, 说明当前单纯形的形状不理想, 需要进行缩小操作:

    \[x_i=x_1+\sigma(x_i-x_1),\qquad i=2,3,\cdots,n+1\]

    然后重复上述过程.

参数的选择有一个经典值:

\[\alpha=1,\qquad\gamma=2,\qquad\rho=0.5,\qquad\sigma=0.5\]

这些默认值在很多实际问题中表现良好, 但在特定问题中可能需要调整. 实现时同时检查单纯形直径和顶点函数值差, 并设置最大迭代或函数评估次数. 即使单纯形缩到一点, 也不能一般地据此证明已到达驻点或全局极小值; Nelder-Mead 有收敛到非极小点的反例.5

References

  1. SciPy, Akima1DInterpolator, C1 光滑性说明及 Notes 中的节点导数权重公式. ↩︎

  2. NIST, DLMF §3.4 Differentiation, 插值微分与三点公式. ↩︎

  3. NIST, DLMF §3.5 Quadrature, Simpson、Romberg 与 Gaussian 求积. ↩︎ ↩︎

  4. Netlib, zeroin.f, Brent 算法的实现, 包含割线退化、逆二次插值及步长保护. ↩︎

  5. J. C. Lagarias, J. A. Reeds, M. H. Wright and P. E. Wright, Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions, SIAM Journal on Optimization 9(1), 112–147 (1998), 引言及算法描述. ↩︎