数值计算相关论文&书籍中文翻译

5.0 引言

本章的目的是向你介绍一些常用于函数求值的技术。在第6章中，我们将通过给出多种具体函数的例程来应用并说明这些技术。因此，本章与下一章的目的在很大程度上是一致的。然而，偶尔在第6章中，某个特殊函数的最佳求值方法可能仅适用于该函数本身。通过比较本章与下一章的内容，你应该能够大致了解在实际应用中“通用”方法与“专用”方法之间的平衡关系。

在上述平衡倾向于通用方法的情况下，本章将为你提供一些思路，帮助你编写自己的函数求值例程——即使该函数对你而言是“特殊”的，但还不足以被纳入第6章或标准函数库中。

参考文献

Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: Prentice-Hall).

Lanczos, C. 1956, Applied Analysis; 1988年重印版 (New York: Dover), 第7章。

5.1 多项式与有理函数

一个 $N$ 次多项式在数值上表示为一个存储系数的数组 $\mathbf{c}[j]$，其中 $j=0,\dots ,N$。我们始终将 $\mathbf{c}[0]$ 视为多项式的常数项，$\mathbf{c}[N]$ 为 $x^N$ 的系数；当然，也可以采用其他约定。对多项式可以进行两类操作：数值操作（例如求值），此时给定其自变量的数值；或代数操作，此时希望以某种方式变换系数数组，而不指定具体的自变量值。我们先从数值操作开始。

我们假定你已经知道绝不能用如下方式求值多项式：

$$\begin{array}{r}\mathrm{p}=\mathrm{c}\left[0\right]+\mathrm{c}\left[1\right]\ast \mathrm{x}+\mathrm{c}\left[2\right]\ast \mathrm{x}\ast \mathrm{x}+\mathrm{c}\left[3\right]\ast \mathrm{x}\ast \mathrm{x}\ast \mathrm{x}+\mathrm{c}\left[4\right]\ast \mathrm{x}\ast \mathrm{x}\ast \mathrm{x}\ast \mathrm{x};\end{array}$$

或者（更糟糕的！）：

p = c[0] + c[1]*x + c[2]*pow(x,2.0) + c[3]*pow(x,3.0) + c[4]*pow(x,4.0);

等到（计算机）革命来临之时，所有被发现犯有此类罪行的人都将被立即处决，而他们的程序却不会！然而，至于是写成：

$$p=c[0]+x\ast (c[1]+x\ast (c[2]+x\ast (c[3]+x\ast c[4])));$$

还是写成：

$$\begin{array}{r}\mathsf{p}=(((c[4]\ast x+c[3])\ast x+c[2])\ast x+c[1])\ast x+c[0];\end{array}$$

则属于个人偏好的问题。

如果系数 $c[\mathrm{0..}n-1]$ 的数量很大，通常会写成：

$$\begin{array}{r}\mathrm{p}=\mathrm{c}[\mathrm{n}-1];\\ \mathrm{for}(\mathrm{j}=\mathrm{n}-2;\mathrm{j}>=0;\mathrm{j}--)\mathrm{p}=\mathrm{p}\ast \mathrm{x}+\mathrm{c}[\mathrm{j}];\end{array}$$

或者：

$$\begin{array}{r}\mathrm{p}=\mathrm{c}[\mathrm{j}=\mathrm{n}-1];\\ \mathrm{while}(\mathrm{j}>0)\mathrm{p}=\mathrm{p}\ast \mathrm{x}+\mathrm{c}[--\mathrm{j}];\end{array}$$

我们可以将其形式化为定义一个函数对象（或称仿函数），该对象绑定一个系数数组的引用，并为其赋予多项式求值功能：

struct Poly {
    // Polynomial function object that binds a reference to a vector of coefficients.
    VecDoub &c;
    Poly(VecDoub &cc) : c(cc) {}
    Doub operator() (Doub x) {
        // Returns the value of the polynomial at x.
        Int j;
        Doub p = c[j = c.size() - 1];
        while (j > 0) p = p * x + c[--j];
        return p;
    }
};

这样你就可以写出如下形式的代码：

y = Poly(c)(x);

其中 $c$ 是一个系数向量。

另一个有用的技巧是同时计算多项式 $P(x)$ 及其导数 $dP(x)/dx$：

或者：

该技巧（本质上是综合除法 [1,2]）可推广为同时计算多项式及其前 $nd$ 阶导数：

给定一个次数为 $nc$ 的多项式的系数数组 $c[\mathrm{0..}nc]$（大小为 $nc+1$，其中 $c[0]$ 为常数项）以及一个数值 $x$，该例程将输出数组 $pd$（大小为 $nd+1$）填充为：$pd[0]$ 为多项式在 $x$ 处的值，$pd[\mathrm{1..}nd]$ 为在 $x$ 处的前 $nd$ 阶导数值。

Int nnd, j, i, nc = c.size() - 1, nd = pd.size() - 1;
Doub cnst = 1.0;
pd[0] = c[nc];
for (j = 1; j < nd + 1; j++) pd[j] = 0.0;
for (i = nc - 1; i >= 0; i--) {
    nnd = (nd < (nc - i) ? nd : nc - i);
    for (j = nnd; j > 0; j--) pd[j] = pd[j] * x + pd[j - 1];
    pd[0] = pd[0] * x + c[i];
}
for (i = 2; i < nd + 1; i++) {
    cnst -= i;
    pd[i] *= cnst;
}

作为趣闻，你或许会感兴趣：对于次数 $n>3$ 的多项式，若允许预先计算某些辅助系数，并在某些情况下增加一些加法运算，则其求值所需的乘法次数可以少于 $n$ 次。例如，多项式

$$P(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4$$

其中 $a_4>0$，可通过以下方式仅用三次乘法和五次加法进行求值：

$$P(x)=\left[(Ax+B{)}^2+Ax+C\right]\left[(Ax+B{)}^2+D\right]+E$$

其中 $A$、$B$、$C$、$D$ 和 $E$ 需预先计算如下：

$$\begin{array}{rl}& A=(a_4{)}^{1/4}\\ & B=\frac{a_3-A^3}{4A^3}\\ & D=3B^2+8B^3+\frac{a_{1}A-2a_{2}B}{A^2}\\ & C=\frac{a_2}{A^2}-2B-6B^2-D\\ & E=a_0-B^4-B^2(C+D)-CD\end{array}$$

五次多项式可用四次乘法和五次加法求值；六次多项式可用四次乘法和七次加法求值。如果你觉得这些内容有趣，可参考文献 [3–5]。这一主题与 §2.11 中讨论的快速矩阵乘法有相似之处。

现在转向代数运算。你可以通过如下代码将一个 $n-1$ 次多项式（系数数组范围为 [0..n-1]）乘以一个单项因子 $x-a$：

$$\begin{array}{rl}& \text{c[n] = c[n - 1];}\\ & \text{for (j = n - 1; j >= 1; j--) c[j] = c[j - 1] - c[j] * a;}\\ & \text{c[0] = (-a);}\end{array}$$

类似地，你也可以通过如下代码将一个 $n$ 次多项式除以单项因子 $x-a$（即再次使用综合除法）：

rem = c[n];
c[n] = 0.;
for (i = n - 1; i >= 0; i--) {
    swap = c[i];
    c[i] = rem;
    rem = swap + rem * a;
}

该操作将得到一个新的多项式系数数组以及一个数值余项 rem。

两个一般多项式的乘法涉及直接对各项乘积求和，每个乘积包含来自两个多项式各一个系数。而两个一般多项式的除法，虽然可以用纸笔计算时所教的笨拙方法完成，但也可以大幅简化。以下例程基于文献 [3] 中的算法：

void poldiv(VecDoub_I &u, VecDoub_I &v, VecDoub_O &q, VecDoub_O &r)

将多项式 u 除以多项式 v，并将商和余式分别返回到 q 和 r 中。这四个多项式均以系数向量表示，均从常数项开始。对 u 和 v 的相对长度没有限制，且两者均可包含尾部零（即表示次数低于其长度所允许的多项式）。q 和 r 返回时长度与 u 相同，但通常会有尾部零。

{
    Int k, j, n = u.size() - 1, nv = v.size() - 1;
    while (nv >= 0 && v[nv] == 0.) nv--;
    if (nv < 0) throw("poldiv divide by zero polynomial");
    r = u;
    q = VecDoub(u.size(), 0);
    for (k = n - nv; k >= 0; k--) {
        q[k] = r[nv + k] / v[nv];
        for (j = nv + k - 1; j >= k; j--) r[j] -= q[k] * v[j - k];
    }
    for (j = nv; j <= n; j++) r[j] = 0.0;
}

5.1.1 有理函数

有理函数形如

$$R(x)=\frac{P_{\mu }(x)}{Q_{\nu }(x)}=\frac{p_0+p_{1}x+\cdots +p_{\mu }{x}^{\mu }}{q_0+q_{1}x+\cdots +q_{\nu }{x}^{\nu }}$$

其求值方法显而易见，即分别计算分子和分母两个多项式后再相除。按惯例，通常选择 $q_0=1$，这可通过将分子分母同时除以任意其他 $q_0$ 实现。此时，通常将两组系数（省略 $q_0$）按如下顺序存入一个数组中会更方便：

$$(p_0,p_1,\dots ,p_{\mu },q_1,\dots ,q_{\nu })$$

以下结构体封装了一个有理函数。它提供了两种构造方式：一种是从独立的分子和分母多项式构造，另一种是从形如 (5.1.5) 的单一数组构造，并显式指定 $\mathbf{n}=\mu +1$ 和 $\mathbf{d}=\nu +1$。其求值函数使 Ratfn 成为一个仿函数（functor），类似于 Poly。我们将在 §5.12 和 §5.13 中使用该对象。

struct Ratfn {
    // 有理函数的函数对象。
    VecDoub cofs;
    Int nn, dd; // 分子、分母的系数个数。

    Ratfn(VecDoub_I &num, VecDoub_I &den) : cofs(num.size() + den.size() - 1),
                                            nn(num.size()), dd(den.size()) {
        // 由分子、分母多项式（系数向量）构造。
        Int j;
        for (j = 0; j < nn; j++) cofs[j] = num[j] / den[0];
        for (j = 1; j < dd; j++) cofs[j + nn - 1] = den[j] / den[0];
    }
    Ratfn(VecDoub_I &cofs, const Int n, const Int d) : cofs(cofs), nn(n), dd(d) {}
    // 由已归一化并存于单一数组中的系数构造。

    Doub operator()(Doub x) const {
        // 在 x 处求值有理函数并返回结果。
        Int j;
        Doub sumn = 0., sumd = 0.;
        for (j = nn - 1; j >= 0; j--) sumn = sumn * x + cofs[j];
        for (j = nn + dd - 2; j >= nn; j--) sumd = sumd * x + cofs[j];
        return sumn / (1.0 + x * sumd);
    }
};

5.1.2 多项式的并行求值

一个 $N$ 次多项式可在约 ${\log }_{2}N$ 个并行步骤内完成求值 [6]。以下以 $N=5$ 为例说明该方法。首先从系数向量开始，并想象其后补零：

$$c_0,c_1,c_2,c_3,c_4,c_5,0,\dots$$

然后将元素两两相加，每对中的第二个元素乘以 $x$：

$$c_0+c_{1}x,c_2+c_{3}x,c_4+c_{5}x,0,\dots$$

接着执行相同的操作，但乘数变为 $x^2$：

$$(c_0+c_{1}x)+(c_2+c_{3}x)x^2,(c_4+c_{5}x)+(0)x^2,0\dots$$

最后一次乘以 $x^4$：

$$[(c_0+c_{1}x)+(c_2+c_{3}x)x^2]+[(c_4+c_{5}x)+(0)x^2]x^4,0\dots$$

最终我们得到一个（有效）长度为 1 的向量，其值即为所求的多项式求值结果。可以看出，这些零值仅用于处理有效子向量长度为奇数的情形；在实际实现中，可以避免对这些零值进行大部分运算。这种并行方法通常比标准的串行编码具有更好的舍入误差特性。

参考文献

Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington, DC: Mathematical Association of America), pp. 183, 190.[1]

Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics, 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), pp. 361–363.[2]

Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §4.6.[3]

Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: Prentice-Hall), Chapter 4.

Winograd, S. 1970, “On the number of multiplications necessary to compute certain functions,” Communications on Pure and Applied Mathematics, vol. 23, pp. 165–179.[4]

Kronsjö, L. 1987, Algorithms: Their Complexity and Efficiency, 2nd ed. (New York: Wiley).[5]

Estrin, G. 1960, quoted in Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §4.6.4.[6]

5.2 连分式的求值

连分式通常是科学应用中计算函数的一种强有力的方法。一个连分式具有如下形式：

$$f(x)=b_0+\frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3+\frac{a_4}{b_4+\frac{a_5}{b_5+\cdots }}}}}$$

排版人员更倾向于将其写作：

$$f(x)=b_0+\frac{a_1}{b_1+}\frac{a_2}{b_2+}\frac{a_3}{b_3+}\frac{a_4}{b_4+}\frac{a_5}{b_5+}\cdots$$

在式 (5.2.1) 或 (5.2.2) 中，$a$ 和 $b$ 本身可以是 $x$ 的函数，通常最多是线性或二次单项式（即常数乘以 $x$ 或 $x^2$）。例如，正切函数的连分式表示为：

$$\tan x=\frac{x}{1-}\frac{x^2}{3-}\frac{x^2}{5-}\frac{x^2}{7-}\cdots$$

连分式通常比幂级数展开收敛得更快，并且在复平面上的收敛域通常也更大（尽管不一定包含幂级数的收敛域）。有时连分式在幂级数表现最差的地方收敛得最好，但这并非普遍规律。Blanch [1] 对连分式最有用的收敛性判别法进行了很好的综述。

存在一些标准技术（包括重要的商差算法），用于在连分式近似、幂级数近似和有理函数近似之间相互转换。关于这一主题的入门介绍可参阅 Acton [2]，更详细的论述和参考文献见 Fike [3]。

在计算连分式时，如何判断何时停止？与级数不同，你不能简单地从左到右计算式 (5.2.1)，并在变化很小时停止。若按 (5.2.1) 的形式书写，唯一的方法是从右向左计算，首先（盲目地！）猜测从何处开始截断。这不是正确的方法。

正确的方法是利用连分式与有理逼近之间的关系，从而能够从左到右计算式 (5.2.1) 或 (5.2.2)。令 $f_n$ 表示用系数 $a_n$ 和 $b_n$ 计算式 (5.2.2) 所得的结果，则

$$f_n=\frac{A_n}{B_n}$$

其中 $A_n$ 和 $B_n$ 由以下递推关系给出：

$$\begin{array}{rl}& A_{-1}\equiv 1B_{-1}\equiv 0\\ & A_0\equiv b_0{B}_0\equiv 1\\ & A_j=b_j{A}_{j-1}+a_j{A}_{j-2}{B}_j=b_j{B}_{j-1}+a_j{B}_{j-2}j=1,2,\dots ,n\end{array}$$

该方法由 J. Wallis 于 1655 年提出（！），并在其著作《Arithmetica Infinitorum》[4] 中进行了讨论。你可以很容易地用数学归纳法证明它。

在实际应用中，该算法存在一些缺点：递推式 (5.2.5) 经常会产生非常大或非常小的部分分子 $A_j$ 和分母 $B_j$，从而可能导致浮点数表示的溢出或下溢。然而，递推式 (5.2.5) 对 $A$ 和 $B$ 是线性的。在任意时刻，你都可以对当前保存的递推式的两层进行重新缩放，例如将 $A_j$、$B_j$、$A_{j-1}$ 和 $B_{j-1}$ 都除以 $B_j$。这顺便使得 $A_j=f_j$，便于检验是否已达到足够精度：只需检查上一次迭代得到的 $f_j$ 和 $f_{j-1}$ 是否足够接近。如果 $B_j$ 恰好为零（这种情况可能发生），则跳过本次循环的归一化操作。更高级的优化策略是仅在即将发生溢出时才进行归一化，以避免不必要的除法运算。事实上，C 语言库函数 ldexp 可用于完全避免除法。（参见 §6.5 末尾的示例。）

近年来提出了两种新的连分式求值算法。Steed 方法不显式使用 $A_j$ 和 $B_j$，而只使用比值 $D_j=B_{j-1}/B_j$。通过以下递推式计算 $D_j$ 和 $\Delta f_j=f_j-f_{j-1}$：

$$\begin{array}{r}{D}_j=1/(b_j+a_j{D}_{j-1})\\ \Delta f_j=(b_j{D}_j-1)\Delta f_{j-1}\end{array}$$

Steed 方法（参见，例如 [5]）避免了对中间结果进行重新缩放的需要。然而，对于某些连分式，有时会遇到（5.2.6）式中分母趋近于零的情况，导致 $D_j$ 和 $\Delta f_j$ 非常大。下一个 $\Delta f_{j+1}$ 通常会抵消这一大的变化，但会导致 $f_j$ 的数值累加和精度损失。围绕这一问题进行编程处理较为棘手，因此 Steed 方法仅推荐用于那些事先已知分母不会为零的情形。我们将在 besselik 例程（§6.6）中为特定目的使用该方法。

目前看来，计算连分式最好的通用方法是修正的 Lentz 方法 [6]。该方法通过同时使用比值

$$C_i=A_i/A_{i-1},D_i=B_{i-1}/B_i$$

并按如下方式计算 $f_i$：

$$f_i=f_{i-1}{C}_j{D}_j$$

由方程 (5.2.5) 可以很容易地证明，这些比值满足递推关系：

$$D_i=1/(b_i+a_i{D}_{i-1}),C_j=b_j+a_j/C_{i-1}$$

在此算法中，存在一种风险：$D_j$ 表达式中的分母或 $C_j$ 本身可能趋近于零。这两种情况都会使 (5.2.10) 式失效。然而，Thompson 和 Barnett [5] 展示了如何修改 Lentz 算法来解决这一问题：只需将有问题的项加上一个很小的量，例如 $10^{-30}$。如果你按照这一规则执行算法的一个循环，就会发现 $f_{j+1}$ 能被准确计算。

具体而言，修正的 Lentz 算法如下：

$$\begin{array}{rl}& \bullet \text{ 设 }{f}_0=b_0\text{；若 }{b}_0=0\text{，则设 }{f}_0=\text{tiny}。\\ & \bullet \text{ 设 }{C}_0=f_0。\\ & \bullet \text{ 设 }{D}_0=0。\\ & \bullet \text{ 对 }j=1,2,\dots \text{ 执行：}\\ & \text{计算 }{D}_j=1/(b_j+a_j{D}_{j-1})\text{；若 }{D}_j=0\text{，则设 }{D}_j=\text{tiny}。\\ & \text{计算 }{C}_j=b_j+a_j/C_{j-1}\text{；若 }{C}_j=0\text{，则设 }{C}_j=\text{tiny}。\\ & \text{计算 }\Delta f_j=C_j{D}_j。\\ & \text{更新 }{f}_j=f_{j-1}\Delta f_j。\\ & \text{若 }|\Delta f_j-1|<\text{eps}\text{，则终止。}\end{array}$$

其中 eps 是你的浮点精度，例如 $10^{-7}$ 或 $10^{-15}$。参数 tiny 应小于 eps $|b_i|$ 的典型值，例如 $10^{-30}$。

上述算法假设当 $|f_j-f_{j-1}|$ 足够小时，即可终止连分式的计算。这通常是成立的，但绝非绝对保证。Jones [7] 给出了一系列定理，可用于为各种连分式证明该终止准则的合理性。

目前尚无对 Lentz 算法中误差传播的严格分析。然而，经验测试表明，其效果至少与其他方法相当。

5.2.1 连分式的变换

连分式的若干重要性质可用于将其改写为能加速数值计算的形式。等价变换

$$a_n\to \lambda a_n,b_n\to \lambda b_n,a_{n+1}\to \lambda a_{n+1}$$

不会改变连分式的值。通过适当选择缩放因子 $\lambda$，通常可以简化 $a$ 和 $b$ 的形式。当然，你可以对连分式的连续项依次进行多次等价变换，每次使用不同的 $\lambda$。

连分式的偶部和奇部分别是其收敛子为 $f_{2n}$ 和 $f_{2n+1}$ 的连分式。它们的主要用途在于其收敛速度是原连分式的两倍，因此，如果它们的项并不比原连分式复杂太多，就能大幅节省计算量。对于 (5.2.2) 式，其偶部的公式为

$$f_{\mathrm{even}}=d_0+\frac{c_1}{d_1+}\frac{c_2}{d_2+}\cdots$$

其中借助中间变量

$$\begin{array}{rl}& {\alpha }_1=\frac{a_1}{b_1}\\ & {\alpha }_n=\frac{a_n}{b_n{b}_{n-1}},n\ge 2\end{array}$$

我们有

$$\begin{array}{r}{d}_0=b_0,c_1={\alpha }_1,d_1=1+{\alpha }_2\\ c_n=-{\alpha }_{2n-1}{\alpha }_{2n-2},d_n=1+{\alpha }_{2n-1}+{\alpha }_{2n},n\ge 2\end{array}$$

奇部的类似公式可在 Blanch 的综述 [1] 中找到。通常将变换 (5.2.14) 与 (5.2.11) 结合使用，以获得最适合数值计算的形式。

我们将在下一章频繁使用连分式。

参考文献

Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions (Washington: National Bureau of Standards); 1968 年重印 (New York: Dover); 在线版见 http://www.nr.com/aands, §3.10。

Blanch, G. 1964, "Numerical Evaluation of Continued Fractions," SIAM Review, vol. 6, pp. 383–421.[1]

Acton, F.S. 1970, Numerical Methods That Work, 1990 年修订版 (Washington, DC: Mathematical Association of America), 第 11 章。[2]

Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: North-Holland), 第 1 章。

Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: Prentice-Hall), §8.2, §10.4, and §10.5.[3]

Wallis, J. 1695, in Opera Mathematica, vol. 1, p. 355, Oxoniae e Theatro Shedoniano. Reprinted by Georg Olms Verlag, Hildesheim, New York (1972).[4]

Thompson, I.J., and Barnett, A.R. 1986, "Coulomb and Bessel Functions of Complex Arguments and Order," Journal of Computational Physics, vol. 64, pp. 490–509.[5]

Lentz, W.J. 1976, "Generating Bessel Functions in Mie Scattering Calculations Using Continued Fractions," Applied Optics, vol. 15, pp. 668–671.[6]

Jones, W.B. 1973, in Padé Approximants and Their Applications, P.R. Graves-Morris, ed. (London: Academic Press), p. 125.[7]

5.3 级数及其收敛性

众所周知，一个解析函数可以在某点 $x_0$ 的邻域内展开为幂级数：

$$f(x)=\sum _{k=0}^{\infty }{a}_k(x-x_0{)}^k$$

这类级数的求值是直接的。当然，你不会对每一项都从头计算 $(x-x_0{)}^k$；而是保留 $(x-x_0{)}^{k-1}$，再通过一次乘法更新它。类似地，系数 $a_k$ 的形式通常也便于利用前面的计算结果：诸如 $k!$ 或 $(2k)!$ 这样的项只需一两次乘法即可更新。

如何判断何时已求和足够多的项？在实践中，各项最好能迅速变小，否则该级数一开始就不适合作为计算方法。尽管在所有情况下这并非数学上严格的做法，但标准做法是：当刚刚加入的项的绝对值小于某个小量 $ϵ$ 乘以当前累计和的绝对值时，就停止求和。（但要注意，若某些孤立的 $a_k=0$ 是可能的！）

有时，即使计算效率不高，你仍可能希望用级数表示来计算某个函数。例如，你可能用所得数值去拟合一个近似表达式，以便后续使用（参见 §5.8）。如果你需要对大量缓慢收敛的项求和，请注意舍入误差！在浮点表示中，从小到大依次累加数值比从大到小更精确。更好的做法是将项两两分组，再将这些组两两合并，依此类推，使得每次加法操作的两个操作数大小相当。

幂级数表示的一个弱点是：它保证不会收敛到复平面上离 $x_0$ 更远的位置——一旦在该距离处遇到奇点，级数就不再收敛。这种“灾难”通常并不意外：当你从书中找到一个幂级数（或自己推导出一个）时，通常也会知道其收敛半径。然而，还有一种更隐蔽的问题：某些级数在数学意义上处处收敛，但在数值方法中几乎处处收敛得不够快，无法实用。两个熟悉的例子是正弦函数和第一类贝塞尔函数：

$$\sin x=\sum _{k=0}^{\infty }\frac{(-1{)}^k}{(2k+1)!}{x}^{2k+1}$$

$$J_n(x)={\left(\frac{x}{2}\right)}^n\sum _{k=0}^{\infty }\frac{(-\frac{1}{4}{x}^2{)}^k}{k!(k+n)!}$$

这两个级数对所有 $x$ 都收敛。但它们在 $k\gg |x|$ 之前甚至不会开始收敛；在此之前，各项的绝对值是递增的。更糟的是，各项符号交替，导致在有限精度算术下产生严重的相消误差。这使得这些级数在 $x$ 较大时完全无法使用。

5.3.1 发散级数

发散级数常常非常有用。一类是超出其收敛半径的幂级数，通常可通过下文将要介绍的加速技术进行求和。另一类是渐近级数，例如由欧拉积分（与指数积分 $E_1$ 相关）导出的欧拉级数：

$$E(x)={\int }_0^{\infty }\frac{e^{-t}}{1+xt}dt\simeq \sum _{k=0}^{\infty }(-1{)}^{k}k!x^k$$

该级数是通过对 $(1+xt{)}^{-1}$ 按 $x$ 的幂展开并逐项积分得到的。该级数对所有 $x\ne 0$ 均发散。当 $x=0.1$ 时，该级数在发散前仅能给出三位有效数字。然而，借助收敛加速技术，我们仍可轻松计算函数 $E(x)$，即使在 $x\sim 2$（此时级数剧烈发散）的情况下也是如此！

5.3.2 级数收敛的加速

有若干技巧可用于加速级数（或等价地，部分和序列）

$$s_n=\sum _{k=0}^n{a}_k$$

的收敛速度。（在本节中，我们将“序列”和“级数”这两个术语互换使用。）Weniger [1] 对此作了极佳的综述。在介绍这些技巧及其适用场合之前，我们需要对序列收敛的方式进行分类。假设 $s_n$ 收敛于 $s$，且

$$\underset{n\to \infty }{\lim }\frac{a_{n+1}}{a_n}=\rho$$

若 $0<|\rho |<1$，称其为线性收敛；若 $\rho =1$，则为对数收敛；若 $\rho =0$，则为超线性收敛。当然，若 $|\rho |>1$，则序列发散。（更严格地说，这种分类应基于所谓的余项 $s_n-s$ [1]。然而，我们的定义更具实用性，且在将对数收敛情形限制为同号项时是等价的。）

线性收敛的典型例子是几何级数：

$$s_n=\sum _{k=0}^n{x}^k=\frac{1-x^{n+1}}{1-x}$$

容易看出 $\rho =x$，因此当 $0<|x|<1$ 时为线性收敛。对数收敛的典型例子是黎曼 zeta 函数的级数：

$$\xi (x)=\sum _{k=1}^{\infty }\frac{1}{k^x},x>1$$

该级数以收敛速度极慢而著称，尤其当 $x\to 1$ 时更为明显。级数 (5.3.2) 和 (5.3.3)，或 $e^x$ 的级数，是超线性收敛的典型例子。我们注意到，超线性收敛并不一定意味着该级数对所有 $x$ 值都易于计算。有时，只有当级数项开始减小时，加速收敛才有帮助。

最著名的用于加速收敛的级数变换方法可能是欧拉变换（Euler transformation）（参见，例如 [2,3]），它最早可追溯到 1755 年。欧拉变换适用于交错级数（即级数项符号交替变化的情况）。通常建议先直接计算前若干项（例如前 $n-1$ 项），然后从第 $n$ 项开始对剩余级数应用该变换。其公式（当 $n$ 为偶数时）为：

$$\sum _{s=0}^{\infty }(-1{)}^s{a}_s=a_0-a_1+a_2\dots -a_{n-1}+\sum _{s=0}^{\infty }\frac{(-1{)}^s}{2^{s+1}}[{\Delta }^s{a}_n]$$

其中 $\Delta$ 是前向差分算子，即：

$$\begin{array}{rl}& \Delta a_n\equiv a_{n+1}-a_n\\ & {\Delta }^2{a}_n\equiv a_{n+2}-2a_{n+1}+a_n\\ & {\Delta }^3{a}_n\equiv a_{n+3}-3a_{n+2}+3a_{n+1}-a_n\end{array}$$

当然，你实际上并不会计算公式 (5.3.9) 右侧的无穷级数，而只会计算前 $p$ 项（例如），这需要从 $a_n$ 开始计算前 $p$ 阶差分（5.3.10）。van Wijngaarden [6] 提出了一种优雅而巧妙的欧拉变换实现方法，详见 Webnote [7]。

欧拉变换是线性变换的一个例子：变换后级数的部分和是原级数部分和的线性组合。尽管欧拉变换及其他线性变换在理论上仍然重要，但通常已被更新、更强大的非线性变换方法所取代。正如数值计算中的一般规律：天下没有免费的午餐。虽然非线性变换更强大，但它们比线性变换风险更高，偶尔会彻底失败。但如果你遵循下文的指导，我们认为你将不再需要使用那些弱小的线性变换。

最早的非线性序列变换例子是 Aitken 的 ${\Delta }^2$-过程。若 $s_n,s_{n+1},s_{n+2}$ 是三个连续的部分和，则改进的估计值为：

$$s_n^{\prime }\equiv s_n-\frac{(s_{n+1}-s_n{)}^2}{s_{n+2}-2s_{n+1}+s_n}=s_n-\frac{(\Delta s_n{)}^2}{{\Delta }^2{s}_n}$$

公式 (5.3.11) 对几何级数是精确的，这也是推导该公式的一种方式。如果你构造出序列 $s_i^{\prime }$，可以再次对这个新序列应用 (5.3.11)，依此类推。（实践中，这种迭代在第一阶段之后通常效果甚微。）注意，公式 (5.3.11) 应按原式计算；存在一些代数上等价但更容易受舍入误差影响的形式。

Aitken 的 ${\Delta }^2$-过程仅适用于线性收敛的序列。与欧拉变换类似，它也已被我们接下来将描述的两种算法所取代。在给出这些算法的例程后，我们将提供一些关于何时使用它们的经验法则。

第一个“现代”非线性变换由 Shanks 提出。Wynn 给出了一种高效的递归实现，称为 $ϵ$ 算法。Aitken 的 ${\Delta }^2$-过程是 $ϵ$ 算法的一个特例，对应于每次仅使用三项。尽管此处不给出推导过程，但很容易说明 $ϵ$ 算法的作用：如果你输入一个幂级数的部分和，$ϵ$ 算法将返回在该幂级数所用 $x$ 值处计算的“对角线”Padé 逼近（§5.12）的数值结果。（逼近式本身的系数并未计算。）也就是说，若 $[M/N]$ 表示分子为 $M$ 次多项式、分母为 $N$ 次多项式的 Padé 逼近，则该算法返回以下逼近的数值：

$$[0/0],[1/0],[1/1],[2/1],[2/2],[3/2],[3/3]\dots$$

（下面的 Epsalg 对象大致等价于 §5.12 中的 pade 函数，后接对所得有理函数的求值。）

在基于 [1] 中例程实现的 Epsalg 对象中，你逐项提供序列，并在调用程序中监控输出以判断是否收敛。例程内部包含除零检查，若发生除零则用一个大数代替结果。该检查可能由以下三种情况触发：(i) 最可能的情况是算法已经收敛，本应更早停止；(ii) 出现了“偶然”的零项，程序可以自行恢复；(iii) 在实践中几乎不会发生，即由于项的特殊组合导致算法真正失败。由于 (i) 和 (ii) 远比 (iii) 常见，Epsalg 隐藏了检查条件，而是返回最近一次已知的良好估计值。

// struct Epsalg {
// 使用 $\epsilon$ 算法对序列进行收敛加速。
// 通过构造函数初始化，参数 nmax 为待求和项数的上限，epss 为期望精度。
// 然后依次调用 next 函数，传入序列的下一个部分和。
// next 返回当前对序列极限的估计值。
// 当检测到收敛时，标志 cnvgd 被置为 true。

VecDoub e;
Int n, ncv;
Bool cnvgd;
Doub eps, small, big, lastval, lasteps;

Epsalg(Int nmax, Doub epss) : e(nmax), n(0), ncv(0),
cnvgd(false), eps(epss), lastval(0.) {
    small = numeric_limits<Doub>::min() * 10.0;
    big = numeric_limits<Doub>::max();
}

Doub next(Doub sum) {
    Doub diff, temp1, temp2, val;
    e[n] = sum;
    temp2 = 0.0;
    for (Int j = n; j > 0; j--) {
        temp1 = temp2;
        temp2 = e[j-1];
        diff = e[j] - temp2;
        if (abs(diff) <= small)
            e[j-1] = big;
        else
            e[j-1] = temp1 + 1.0 / diff;
    }

    n++;
    val = (n & 1) ? e[0] : e[1];
    if (abs(val) > 0.01 * big) val = lastval;
    lasteps = abs(val - lastval);
    if (lasteps > eps) ncv = 0;
    else {
        ncv++;
        if (ncv >= 3) cnvgd = true;
    }
    return (lastval = val);
}

上述最后几行实现了一个简单的收敛判断准则。对于收敛稳健的问题，你可以简单地将对 next 的调用放在如下 while 循环中：

Doub val, partialsum, eps = ...;
Epsalg mysum(1000, eps);
while (!mysum.cnvgd) {
    partialsum = ...;
    val = mysum.next(partialsum);
}

对于更敏感的情况，你可以忽略 cnvgd 标志，持续调用 next 直到你对收敛程度满意为止。

一大类现代非线性变换可通过“模型序列”的概念推导出来。其思想是选择一个“简单”的序列来近似给定序列的渐近形式，并构造一个能精确求和该模型序列的变换。可以预期，该变换对具有类似渐近性质的其他序列也会表现良好。例如，几何级数就是 Aitken 的 ${\Delta }^2$-过程所采用的模型序列。

Levin 变换 可能是目前已知最好的单一序列加速方法。它基于用如下形式的表达式对序列进行渐近逼近：

$$s_n=s+{\omega }_n\sum _{j=0}^{k-1}\frac{c_j}{(n+\beta {)}^j}$$

其中 ${\omega }_n$ 是序列余项中的主导项：

$$s_n-s={\omega }_n[c+O(n^{-1})],n\to \infty$$

常数 $c_j$ 是任意的，而 $\beta$ 是一个被限制为正数的参数。Levin 指出，对于形如 (5.3.13) 的模型序列，以下变换给出了级数的精确值：

$$s=\frac{\sum _{j=0}^k(-1{)}^j(\genfrac{}{}{0ex}{}{k}{j})\frac{(\beta +n+j{)}^{k-1}}{(\beta +n+k{)}^{k-1}}\frac{s_{n+j}}{{\omega }_{n+j}}}{\sum _{j=0}^k(-1{)}^j(\genfrac{}{}{0ex}{}{k}{j})\frac{(\beta +n+j{)}^{k-1}}{(\beta +n+k{)}^{k-1}}\frac{1}{{\omega }_{n+j}}}$$

（分子和分母中的公共因子 $(\beta +n+k{)}^{k-1}$ 可降低在 $k$ 较大时发生溢出的可能性。）方程 (5.3.15) 的推导见网络附注 [4]。

(5.3.15) 中的分子和分母并非按所写形式直接计算。相反，它们可以通过一个具有不同初始值的单一递推关系高效计算（推导见 [1]）：

$$D_{k+1}^n(\beta )=D_k^{n+1}(\beta )-\frac{(\beta +n)(\beta +n+k{)}^{k-1}}{(\beta +n+k+1{)}^k}{D}_k^n(\beta )$$

初始值为

$$D_0^n(\beta )=\{\begin{array}{ll}{s}_n/{\omega }_n,& \text{分子}\\ 1/{\omega }_n,& \text{分母}\end{array}$$

尽管 $D_k^n$ 是一个二维对象，但该递推关系可以通过一维数组沿反对角线 $n+k=$ 常数方向进行编码实现。

选择 (5.3.14) 并不能唯一确定 ${\omega }_n$，但如果你对级数有解析信息，这正是可以加以利用的地方。通常你不会这么幸运，此时可以根据启发式方法进行选择。例如，交错级数的余项大约为第一个被忽略项的一半，这提示将 ${\omega }_n$ 设为 $a_n$ 或 $a_{n+1}$。这两种选择分别称为 Levin 的 $t$ 变换和 $d$ 变换。类似地，几何级数的余项是部分和 (5.3.7) 与其极限 $1/(1-x)$ 之差。这可以写成 $a_n{a}_{n+1}/(a_n-a_{n+1})$，从而定义了 Levin 的 $v$ 变换。最常用的选择来自用积分近似 $\zeta$ 函数 (5.3.8) 的余项：

$$\sum _{k=n+1}^{\infty }\frac{1}{k^x}\approx {\int }_{n+1}^{\infty }\frac{dk}{k^x}=\frac{(n+1{)}^{1-x}}{x-1}=\frac{(n+1)a_{n+1}}{x-1}$$

这启发我们选择 $(n+\beta )a_n$（Levin 的 $u$ 变换），其中 $\beta$ 通常取为 1。总结如下：

$${\omega }_n=\{\begin{array}{ll}(\beta +n)a_n,& u\text{ 变换}\\ a_n,& t\text{ 变换}\\ a_{n+1},& d\text{ 变换（修正的 }t\text{ 变换）}\\ \frac{a_n{a}_{n+1}}{a_n-a_{n+1}},& v\text{ 变换}\end{array}$$

对于并非部分和的序列（即单个 $a_n$ 未定义的情况），在 (5.3.19) 中将 $a_n$ 替换为 $\Delta s_{n-1}$。

以下是 Levin 变换的例程，同样基于 [1] 中的例程：

struct Levin {
    // 通过 Levin 变换加速序列的收敛。
    // 初始化时调用构造函数，参数 nmax 为要求和项数的上界，epss 为期望精度。
    // 然后连续调用 next 函数，返回序列极限的当前估计值。
    // 当检测到收敛时，标志 cnvgd 被置为 1。

    Levin(Int nmax, Doub epss) : numer(nmax), denom(nmax), n(0), ncv(0), cnvgd(0), eps(epss), lastval(0.) {
        small = numeric_limits<Doub>::min() * 10.0;
        big = numeric_limits<Doub>::max();
    }

    Doub next(Doub sum, Doub omega, Doub beta = 1.) {
        // 参数说明：
        // sum：序列的第 n 个部分和；
        // omega：第 n 个余项估计值 ω_n，通常来自 (5.3.19)；
        // beta：参数，通常设为 1，但有时 0.5 效果更好。
        // 返回序列极限的当前估计值。

        Int j;
        Doub fact, ratio, term, val;
        term = 1.0 / (beta + n);
        denom[n] = term / omega;
        numer[n] = sum * denom[n];
        if (n > 0) {
            ratio = (beta + n - 1) * term;
            for (j = 1; j <= n; j++) {
                fact = (n - j + beta) * term;
                numer[n - j] = numer[n - j + 1] - fact * numer[n - j];
                denom[n - j] = denom[n - j + 1] - fact * denom[n - j];
                term = term * ratio;
            }
        }
        n++;
        val = abs(denom[0]) < small ? lastval : numer[0] / denom[0];
        lasteps = abs(val - lastval);
        if (lasteps <= eps) {
            ncv++;
            if (ncv >= 2) cnvgd = 1;
        }
        return (lastval = val);
    }

private:
    VecDoub numer, denom;
    Int n, ncv;
    Bool cnvgd;
    Doub eps, lastval, lasteps, small, big;
};

你可以像之前对 Epsalg 讨论的那样，选择使用或不使用 cnvgd 标志。

推导序列变换的另一种方法是使用多项式或有理函数逼近级数的外推法，例如 Wynn 的 $\rho$ 算法 [1]。由于这些方法通常无法超越我们已介绍的两种方法，因此不再赘述。

5.3.3 实用提示与示例

目前对非线性序列变换尚无普适的理论理解。因此，大多数实用建议都基于数值实验 [5]。你可能会认为，对一个剧烈发散的级数求和是序列变换中最困难的问题。然而，问题的难度更多取决于各项是否同号，还是符号交错，而非序列本身是否收敛。特别是，各项同号的对数收敛级数通常最难求和。即使是最优的加速方法，在加速对数收敛时也会受到舍入误差的严重影响。你应始终使用双精度计算，并准备好接受有效数字的部分损失。通常会观察到，在达到某个最优项数前收敛性良好，但若继续增加项数反而会导致有效数字的损失。此外，并不存在一种能加速所有对数收敛序列的通用算法。尽管如此，仍有一些经验法则可供参考。

首先，注意在发散级数中，将渐近级数（其项先减小后增大）与其他发散级数（例如在其收敛半径之外的幂级数）区分开是有益的。对于交错级数（无论是收敛的、渐近的，还是发散的幂级数），Levin 的 $u$ 变换几乎总是最佳选择。对于单调线性收敛或单调发散的幂级数，$ϵ$ 算法通常是首选，但 $u$ 变换也常常表现良好。对于对数收敛，$u$ 变换明显最优（$ϵ$ 算法则完全失效）。对于符号不规则或其他非标准特征的级数，$ϵ$ 算法通常相对稳健，常常在其他算法失败时仍能成功。最后，对于单调渐近级数（例如 Ei$(x)$ 的 (6.3.11) 式），直接求和（不进行加速）反而是最佳选择。

$v$ 变换和 $t$ 变换的效果几乎与 $u$ 变换相当，但 $t$ 变换通常在对数收敛情况下失效。

如果你仅有某个序列的少数几个数值项且缺乏理论洞察，盲目应用收敛加速器可能是危险的。该算法

有时会表现出“收敛”，但这种收敛只是表观上的，并非真实的收敛。解决方法是尝试两种不同的变换以进行验证。

由于对正项级数进行收敛加速远比对交错级数困难得多，因此有时将正项级数转换成交错级数是有用的。Van Wijngaarden 给出了一种实现该转换的变换方法 [6]：

$$\sum _{r=1}^{\infty }{v}_r=\sum _{r=1}^{\infty }(-1{)}^{r-1}{w}_r$$

其中

$$w_r\equiv v_r+2v_{2r}+4v_{4r}+8v_{8r}+\cdots$$

公式 (5.3.20) 和 (5.3.21) 将一个简单求和替换为一个二维求和，其中 (5.3.20) 的每一项本身又是一个无穷级数 (5.3.21)。这看起来似乎是一种奇怪的节省计算量的方法！然而，由于 (5.3.21) 中的下标以 2 的幂次迅速增长，通常只需少数几项就能使 (5.3.21) 以极高的精度收敛。不过，你需要能够高效地计算“任意”索引 $r$ 对应的 $v_r$ 值。上文在公式 (5.3.1) 后提到的针对连续 $r$ 值的标准“递推更新”技巧在此无法使用。

一旦通过 Van Wijngaarden 变换生成了交错级数，Levin $d$ 变换在对该级数求和时特别有效 [8]。这种策略对收敛因子 $\rho$ 接近 1 的线性收敛级数最为有用。而对于对数收敛级数，即使经过变换后的级数 (5.3.21) 通常也收敛得太慢，难以在数值计算中应用。

作为调用 Epsalg 或 Levin 例程的一个示例，考虑计算如下积分的问题：

$$I={\int }_0^{\infty }\frac{x}{1+x^2}{J}_0(x)dx=K_0(1)=0.4210244382\dots$$

标准数值积分方法（如 qromo）会失败，因为被积函数具有长振荡尾部，产生交替的正负贡献，趋于相互抵消。计算此类积分的一个好方法是将其拆分为 $J_0(x)$ 相邻零点之间的积分之和：

$$I={\int }_0^{\infty }f(x)dx=\sum _{j=0}^{\infty }{I}_j$$

其中

$$I_j={\int }_{x_{j-1}}^{x_j}f(x)dx,f(x_j)=0,j=0,1,\dots$$

此处取 $x_{-1}$ 为积分下限，在本例中为零。思路是使用 qromb 或高斯求积法计算相对简单的积分 $I_j$，然后对级数 (5.3.23) 进行收敛加速，因为我们预期各项贡献符号交替。对于示例 (5.3.22)，我们甚至不需要 $J_0(x)$ 零点的精确值，只需取 $x_j=(j+1)\pi$ 即可，这在渐近意义下是正确的。以下是代码：

// evex.h
Doub func(const Doub x)
// 被积函数，对应公式 (5.3.22)。
{
    if (x == 0.0)
        return 0.0;
    else {
        Bessel bess;
        return x * bess.jnu(0.0, x) / (1.0 + x * x);
    }
}

Int main_levex(void)
// 本示例程序演示如何使用 Levin u 变换计算振荡积分，见公式 (5.3.22)。
{
    const Doub PI = 3.141592653589793;
    Int nterm = 12;
    Doub beta = 1.0, a = 0.0, b = 0.0, sum = 0.0;
    Levin series(100, 0.0);
    Cout << setw(5) << "N" << setw(19) << "Sum (direct)" << setw(21)
         << "Sum (Levin)" << endl;
    for (Int n = 0; n <= nterm; n++) {
        b += PI;
        Doub s = qromb(func, a, b, 1.e-8);
        a = b;
        sum += s;
        Doub omega = (beta + n) * s;
        // 使用 u 变换。
        Doub ans = series.next(sum, omega, beta);
        Cout << setw(5) << n << fixed << setprecision(14) << setw(21)
             << sum << setw(21) << ans << endl;
    }
    return 0;
}

在 qromb 中将 eps 设为 $1\times 10^{-8}$，当 $n=8$ 时，仅需约 200 次函数求值即可获得 9 位有效数字。若将 qromb 替换为高斯求积例程，函数求值次数可减少一半。注意，$n=8$ 对应积分上限为 $9\pi$，此时被积函数的振幅仍约为 $10^{-2}$。这体现了收敛加速的非凡能力。（关于振荡积分的更多内容，参见 §13.9。）

参考文献

Weniger, E.J. 1989, "Nonlinear Sequence Transformations for the Acceleration of Convergence and the Summation of Divergent Series," Computer Physics Reports, vol. 10, pp. 189-371.[1]

Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions (Washington: National Bureau of Standards); reprinted 1968 (New York: Dover); online at http://www.nr.com/aands, §3.6.[2]

Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics, 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), §2.3.[3]

Numerical Recipes Software 2007, "Derivation of the Levin Transformation," Numerical Recipes Webnote No. 6, at http://www.nr.com/webnotes?6 [4]

Smith, D.A., and Ford, W.F. 1982, "Numerical Comparisons of Nonlinear Convergence Accelerators," Mathematics of Computation, vol. 38, pp. 481-499.[5]

Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), Chapter 13 [van Wijngaarden's transformations].[6]

Numerical Recipes Software 2007, "Implementation of the Euler Transformation," Numerical Recipes Webnote No. 5, at http://www.nr.com/webnotes?5 [7]

Jentschura, U.D., Mohr, P.J., Soff, G., and Weniger, E.J. 1999, "Convergence Acceleration via Combined Nonlinear-Condensation Transformations," Computer Physics Communications, vol. 116, pp. 28–54.[8]

Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall); reprinted 2003 (New York: Dover), Chapter 3.

5.4 递推关系与 Clenshaw 递推公式

许多有用的函数满足递推关系，例如：

$$(n+1)P_{n+1}(x)=(2n+1)x{P}_n(x)-n{P}_{n-1}(x)$$

$$J_{n+1}(x)=\frac{2n}{x}{J}_n(x)-J_{n-1}(x)$$

$$n{E}_{n+1}(x)=e^{-x}-x{E}_n(x)$$

$$\cos n\theta =2\cos \theta \cos (n-1)\theta -\cos (n-2)\theta$$

$$\sin n\theta =2\cos \theta \sin (n-1)\theta -\sin (n-2)\theta$$

其中前三个函数分别为勒让德多项式、第一类贝塞尔函数和指数积分。（记号参见[1]。）这些递推关系可用于将计算方法从两个连续的 $n$ 值扩展到其他值，无论是更大的还是更小的值。

公式 (5.4.4) 和 (5.4.5) 促使我们对三角函数多说几句。如果你的程序运行时间主要花在计算三角函数上，那你很可能做错了什么。对于参数构成线性序列 $\theta ={\theta }_0+n\delta$（$n=0,1,2,\dots$）的三角函数，可通过以下递推式高效计算：

$$\begin{array}{r}\cos (\theta +\delta )=\cos \theta -[\alpha \cos \theta +\beta \sin \theta ]\\ \sin (\theta +\delta )=\sin \theta -[\alpha \sin \theta -\beta \cos \theta ]\end{array}$$

其中 $\alpha$ 和 $\beta$ 是预先计算的系数：

$$\alpha \equiv 2{\sin }^2\left(\frac{\delta }{2}\right)\beta \equiv \sin \delta$$

之所以采用这种方式，而不是使用标准的（等价的）角度和公式，是因为当增量 $\delta$ 很小时，这里的 $\alpha$ 和 $\beta$ 不会损失有效数字。同样，公式 (5.4.6) 中的加法应按照方括号所示的顺序进行。在第12章处理傅里叶变换时，我们将多次使用 (5.4.6)。

另一个偶尔有用的技巧是注意到 $\sin \theta$ 和 $\cos \theta$ 都可以通过一次 $\tan$ 调用计算得到：

$$t\equiv \tan \left(\frac{\theta }{2}\right)\cos \theta =\frac{1-t^2}{1+t^2}\sin \theta =\frac{2t}{1+t^2}$$

因此，如果你同时需要 $\sin$ 和 $\cos$，其代价是一次 $\tan$ 计算加上两次乘法、两次除法和两次加法。在三角函数运算较慢的机器上，这可能节省时间。但请注意，当 $\theta \to \pm \pi$ 时需要特殊处理。此外，许多现代机器的三角函数运算非常快，因此在未经过测试前，不应假定公式 (5.4.8) 更快。

5.4.1 递推关系的稳定性

你需要注意，递推关系在你打算进行的方向（无论是递增 $n$ 还是递减 $n$）上未必对舍入误差稳定。一个三项线性递推关系

$${\gamma }_{n+1}+a_n{\gamma }_n+b_n{\gamma }_{n-1}=0,n=1,2,\dots$$

有两个线性无关的解，记为 $f_n$ 和 $g_n$。其中只有一个对应于你试图生成的函数序列 $f_n$。另一个解 $g_n$ 在你希望进行的方向上可能是指数增长的、指数衰减的，或者指数中性的（例如按某种幂律增长或衰减）。如果它是指数增长的，那么该递推关系在该方向上几乎没有实际用途。例如，当 $x<n$ 时，公式 (5.4.2) 在 $n$ 递增方向上就是这种情况。你无法通过 (5.4.2) 的前向递推来生成高阶贝塞尔函数。

更正式地说，如果

$$f_n/g_n\to 0\mathrm{as}n\to \infty$$

则称 $f_n$ 为递推关系 (5.4.9) 的最小解。像 $g_n$ 这样的非最小解称为主解。最小解（如果存在）是唯一的，但主解不是唯一的——你可以在给定的 $g_n$ 上加上任意倍数的 $f_n$。你可以通过前向递推计算任何主解，但不能计算最小解。（不幸的是，有时你想要的恰恰是最小解。）

Abramowitz 和 Stegun（在他们的引言中！）[1] 给出了一份在递增或递减方向上稳定的递推关系列表。当然，该列表并不包含所有可能的公式。对于某个函数 $f_n(x)$ 的递推关系，你可以通过大约五分钟的人工计算自行测试其稳定性：在你感兴趣的 $x$ 固定值下，不要用 $f_j(x)$ 和 $f_{j+1}(x)$ 的真实值启动递推，而是（首先）分别用 1 和 0，（其次）分别用 0 和 1 作为初始值。对这两种初始条件，沿你希望的方向（从 $j$ 开始递增或递减）生成 10 到 20 项递推序列。观察两个序列对应项之间的差异。如果差异始终保持在单位量级（例如绝对值小于 10），则递推是稳定的。如果差异缓慢增加，则递推可能轻微不稳定，但仍可接受。如果差异灾难性地增长，则说明递推存在指数增长的解。如果你知道你想要的函数实际上对应于这个增长解，那么你仍可保留该递推公式（例如，对于递增 $n$ 时的贝塞尔函数 $Y_h(x)$；参见 §6.5）。如果你不确定你的函数对应于哪个解，那么此时必须拒绝该递推公式。注意，你可以在费力寻找计算两个起始函数 $f_j(x)$ 和 $f_{j+1}(x)$ 的数值方法之前进行此测试：稳定性是递推关系本身的属性，与起始值无关。

另一种测试稳定性的启发式方法是将递推关系替换为一个类似的、具有常系数的线性递推关系。例如，关系式 (5.4.2) 变为

$${\gamma }_{n+1}-2\gamma {\gamma }_n+{\gamma }_{n-1}=0$$

其中 $\gamma \equiv n/x$ 被视为常数。通过尝试形如 $y_n=a^n$ 的解来求解此类递推关系。代入上述递推式得

$$a^2-2\gamma a+1=0\mathrm{或}a=\gamma \pm \sqrt{{\gamma }^2-1}$$

若对所有解 $a$ 都有 $|a|\le 1$，则递推是稳定的。这在 $|\gamma |\le 1$ 或 $n\le x$ 时成立（你可以自行验证）。因此，递推式 (5.4.2) 无法从 $J_0(x)$ 和 $J_1(x)$ 出发，用于计算大 $n$ 时的 $J_n(x)$。

或许此时你希望获得一些关于此主题的严格定理作为保障（尽管我们自己总是采用上述启发式方法）。以下是 Perron [2] 给出的两个定理：

定理 A. 如果在 (5.4.9) 中，当 $n\to \infty$ 时 $a_n\sim a{n}^{\alpha }$，$b_n\sim b{n}^{\beta }$，且 $\beta <2\alpha$，则

$$g_{n+1}/g_n\sim -a{n}^{\alpha },f_{n+1}/f_n\sim -(b/a)n^{\beta -\alpha }$$

且 $f_n$ 是 (5.4.9) 的最小解。

定理 B. 在与定理 A 相同的条件下，但取 $\beta =2\alpha$，考虑特征多项式

$$t^2+at+b=0$$

若 (5.4.14) 的根 $t_1$ 和 $t_2$ 具有不同的模，例如 $|t_1|>|t_2|$，则

$$g_{n+1}/g_n\sim t_1{n}^{\alpha },f_{n+1}/f_n\sim t_2{n}^{\alpha }$$

且 $f_n$ 仍是 (5.4.9) 的最小解。除上述两个定理所述情形外，其他情形无法确定最小解是否存在。（关于递推关系的稳定性，参见 [3]。）

如果你所期望的解恰好是最小解，该如何处理？答案在于那句古老的谚语：“每朵乌云都镶着银边”：如果某个递推关系在一个方向上是灾难性不稳定的，那么该（不希望的）解在反方向上会迅速衰减。这意味着你可以从任意初始值 $f_j$ 和 $f_{j+1}$ 开始，只要沿稳定方向迭代足够多步，就会收敛到你想要的函数序列，只是会乘上一个未知的归一化因子。如果存在其他方式对序列进行归一化（例如，已知 $f_n$ 的和的公式），那么这就可以成为一种实用的函数计算方法。这种方法称为 Miller 算法。一个常被引用的例子 [1,4] 正是这样使用方程 (5.4.2)，并结合归一化公式

$$1=J_0(x)+2J_2(x)+2J_4(x)+2J_6(x)+\cdots$$

顺便提一下，三项递推关系与连分式之间存在重要联系。将递推关系 (5.4.9) 改写为

$$\frac{y_n}{y_{n-1}}=-\frac{b_n}{a_n+y_{n+1}/y_n}$$

从 $n$ 开始迭代该式，得到

$$\frac{y_n}{y_{n-1}}=-\frac{b_n}{a_n-}\frac{b_{n+1}}{a_{n+1}-}\dots$$

Pincherle 定理 [2] 指出，(5.4.18) 收敛当且仅当 (5.4.9) 存在最小解 $f_n$，此时它收敛于 $f_n/f_{n-1}$。这一结果（通常取 $n=1$，并结合某种确定 $f_0$ 的方法）构成了下一章中许多特殊函数实用计算方法的基础。

5.4.2 Clenshaw 递推公式

Clenshaw 递推公式 [5] 是一种优雅而高效的方法，用于计算系数与满足递推关系的函数之积的和，例如

$$f(\theta )=\sum _{k=0}^N{c}_k\cos k\theta \text{或}f(x)=\sum _{k=0}^N{c}_k{P}_k(x)$$

其工作原理如下：假设所求的和为

$$f(x)=\sum _{k=0}^N{c}_k{F}_k(x)$$

且 $F_k$ 满足递推关系

$$F_{n+1}(x)=\alpha (n,x)F_n(x)+\beta (n,x)F_{n-1}(x)$$

其中 $\alpha (n,x)$ 和 $\beta (n,x)$ 为某些函数。现在定义量 $y_k$（$k=N,N-1,\dots ,1$）如下：

$$\begin{array}{rl}{y}_{N+2}=y_{N+1}=0& \\ y_k=\alpha (k,x)y_{k+1}+\beta (k+1,x)y_{k+2}+c_k& (k=N,N-1,\dots ,1)\end{array}$$

若你将方程 (5.4.21) 左边解出 $c_k$，然后显式写出求和式 (5.4.19)，其形式（部分）如下：

$$\begin{array}{rl}& f(x)=\cdots \\ & +[y_8-\alpha (8,x)y_9-\beta (9,x)y_{10}]F_8(x)\\ & +[y_7-\alpha (7,x)y_8-\beta (8,x)y_9]F_7(x)\\ & +[y_6-\alpha (6,x)y_7-\beta (7,x)y_8]F_6(x)\\ & +[y_5-\alpha (5,x)y_6-\beta (6,x)y_7]F_5(x)\\ & +\cdots \\ & +[y_2-\alpha (2,x)y_3-\beta (3,x)y_4]F_2(x)\\ & +[y_1-\alpha (1,x)y_2-\beta (2,x)y_3]F_1(x)\\ & +[c_0+\beta (1,x)y_2-\beta (1,x)y_2]F_0(x)\end{array}$$

注意，我们在最后一行加减了 $\beta (1,x)y_2$。如果你检查 (5.4.22) 中包含因子 $y_8$ 的各项，会发现它们因递推关系 (5.4.20) 而总和为零；类似地，所有其他 $y_k$（直到 $y_2$）也是如此。(5.4.22) 中唯一剩下的项是

$$f(x)=\beta (1,x)F_0(x)y_2+F_1(x)y_1+F_0(x)c_0$$

方程 (5.4.21) 和 (5.4.23) 构成了用于计算求和式 (5.4.19) 的 Clenshaw 递推公式：你首先利用 (5.4.21) 自上而下计算 $y_k$；当计算到 $y_2$ 和 $y_1$ 后，再用 (5.4.23) 得到最终结果。

如上所述的 Clenshaw 递推公式按 $k$ 递减的顺序引入系数 $c_k$。在每一步，所有先前 $c_k$ 的影响都被“记住”为两个系数，分别乘以函数 $F_{k+1}$ 和 $F_k$（最终为 $F_0$ 和 $F_1$）。如果当 $k$ 较大时函数 $F_k$ 很小，而当 $k$ 较小时系数 $c_k$ 很小，则求和可能主要由较小的 $F_k$ 主导。此时，“记住”的系数将涉及精细的抵消，可能导致灾难性的有效数字损失。一个例子是求和平凡级数

$$J_{15}(1)=0\times J_0(1)+0\times J_1(1)+\dots +0\times J_{14}(1)+1\times J_{15}(1)$$

这里微小的 $J_{15}$ 最终被表示为 $J_0$ 和 $J_1$（量级为 1）的一个抵消性线性组合。

此类情况的解决方案是使用一种替代的 Clenshaw 递推公式，该公式以向上的方向纳入系数 $c_k$。相关方程如下：

$$\begin{array}{rl}& y_{-2}=y_{-1}=0\\ & y_k=\frac{1}{\beta (k+1,x)}[y_{k-2}-\alpha (k,x)y_{k-1}-c_k],k=0,1,\dots ,N-1\\ & f(x)=c_N{F}_N(x)-\beta (N,x)F_{N-1}(x)y_{N-1}-F_N(x)y_{N-2}\end{array}$$

除极少数情况外，Clenshaw 递推总是稳定的，无论函数 $F_k$ 的递推在向上或向下方向是否稳定。而需要使用方程 (5.4.25)–(5.4.27) 而非方程 (5.4.21) 和 (5.4.23) 的罕见情形，可通过检测 (5.4.23) 中第一个求和项的操作数是否符号相反且幅值几乎相等来自动识别。

5.4.3 线性递推关系的并行求值

在需要时，线性递推关系可以高度并行化地求值。考虑一般的一阶线性递推关系：

$$u_j=a_j+b_{j-1}{u}_{j-1},j=2,3,\dots ,n$$

初始值为 $u_1=a_1$。为并行化该递推，可采用一种称为递归倍增（recursive doubling）的强大通用策略。将方程 (5.4.28) 分别对 $2j$ 和 $2j-1$ 写出：

$$\begin{array}{r}{u}_{2j}=a_{2j}+b_{2j-1}{u}_{2j-1}\\ u_{2j-1}=a_{2j-1}+b_{2j-2}{u}_{2j-2}\end{array}$$

将第二个方程代入第一个方程以消去 $u_{2j-1}$，得到：

$$u_{2j}=(a_{2j}+a_{2j-1}{b}_{2j-1})+(b_{2j-2}{b}_{2j-1})u_{2j-2}$$

这是一个与 (5.4.28) 形式相同的新递推式，但仅针对偶数下标的 $u_j$，因此只涉及 $n/2$ 项。显然，我们可以递归地继续这一过程，在每一阶段将递推项数减半，直到剩下长度为 1 或 2 的递推式，此时可直接显式求解。每次完成递归的一个子部分后，利用 (5.4.29) 中的第二个方程回填奇数项。在实践中，这比听起来更简单。总运算量与串行求值相同，但可在约 ${\log }_{2}n$ 个并行步骤内完成。

还有一种称为循环约简（cyclic reduction）的递归倍增变体，可用简单的迭代循环而非递归过程实现 [6]。此时，我们首先对所有相邻项 $u_j$ 和 $u_{j-1}$（而不仅限于偶数项）写出递推式 (5.4.28)。像在方程 (5.4.30) 中那样消去 $u_{j-1}$，得到：

$$u_j=(a_j+a_{j-1}{b}_{j-1})+(b_{j-2}{b}_{j-1})u_{j-2}$$

这是一个具有新系数 $a_j^{\prime }$ 和 $b_j^{\prime }$ 的一阶递推式。重复此过程，可依次得到用 $u_{j-2}$、$u_{j-4}$、$u_{j-8}$、$\cdots$ 表示 $u_j$ 的公式。当达到 $u_{j-n}$（其中 $n$ 为 2 的幂）时过程终止，此时对所有 $j$ 有 $u_{j-n}=0$。因此，最后一步直接给出 $u_j$ 等于最后一组 $a_j^{\prime }$。

在循环约简中，每阶段更新的向量 $u_j$ 的长度并非每步减半，而是在全部 ${\log }_{2}n$ 个阶段中仅从 $\sim n$ 减少到 $\sim n/2$。因此，总运算量为 $O(n\log n)$，而递归倍增为 $O(n)$。这一点是否重要，取决于计算机体系结构的具体细节。

二阶递推关系也可以并行化。考虑二阶递推关系：

$$y_j=a_j+b_{j-2}{y}_{j-1}+c_{j-2}{y}_{j-2},j=3,4,\dots ,n$$

初始值为：

$$y_1=a_1,y_2=a_2$$

在此编号方案下，需提供系数 $a_1,\dots ,a_n$、$b_1,\dots ,b_{n-2}$ 和 $c_1,\dots ,c_{n-2}$。将该递推关系改写为如下形式 [6]：

$$\left(\begin{array}{l}{y}_j\\ y_{j+1}\end{array}\right)=\left(\begin{array}{l}0\\ a_{j+1}\end{array}\right)+\left(\begin{array}{ll}0& 1\\ c_{j-1}& b_{j-1}\end{array}\right)\left(\begin{array}{l}{y}_{j-1}\\ y_j\end{array}\right),j=2,\dots ,n-1$$

即：

$${\mathbf{u}}_j={\mathbf{a}}_j+{\mathbf{b}}_{j-1}\cdot {\mathbf{u}}_{j-1},j=2,\dots ,n-1$$

其中：

$${\mathbf{u}}_j=\left(\genfrac{}{}{0ex}{}{y_j}{y_{j+1}}\right),{\mathbf{a}}_j=\left(\genfrac{}{}{0ex}{}{0}{a_{j+1}}\right),{\mathbf{b}}_{j-1}=\left(\begin{array}{cc}0& 1\\ c_{j-1}& b_{j-1}\end{array}\right)$$

且：

$${\mathbf{u}}_1={\mathbf{a}}_1=\left(\genfrac{}{}{0ex}{}{y_1}{y_2}\right)=\left(\genfrac{}{}{0ex}{}{a_1}{a_2}\right)$$

这是关于向量 ${\mathbf{u}}_j$ 的一阶递推关系，可使用上述任一算法求解。唯一区别在于乘法是 $2\times 2$ 矩阵 ${\mathbf{b}}_j$ 的矩阵乘法。在第一次递归调用之后，$\mathbf{a}$ 和 $\mathbf{b}$ 中的零元素将丢失，因此必须为一般的二维向量和矩阵编写例程。注意，该算法并未避免 §5.4.1 中讨论的二阶递推可能存在的不稳定性问题。此外，该算法可自然推广到更高阶的递推关系：一个 $n$ 阶递推关系可表示为涉及 $n$ 维向量和矩阵的一阶递推关系。

参考文献

Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions (Washington: National Bureau of Standards)；1968年重印版（New York: Dover）；在线版见 http://www.nr.com/aands，第xiii, 697页。[1]

Gautschi, W. 1967, “Computational Aspects of Three-Term Recurrence Relations,” SIAM Review, vol. 9, pp. 24–82。[2]

Lakshmikantham, V., and Trigiante, D. 1988, Theory of Difference Equations: Numerical Methods and Applications (San Diego: Academic Press)。[3]

Acton, F.S. 1970, Numerical Methods That Work；1990年修订版（Washington, DC: Mathematical Association of America），第20页起。[4]

Clenshaw, C.W. 1962, Mathematical Tables, vol. 5, National Physical Laboratory (London: H.M. Stationery Office)。[5]

Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall)；2003年重印版（New York: Dover），§4.4.3，第111页。

Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library)，第76页。

Hockney, R.W., and Jesshope, C.R. 1988, Parallel Computers 2: Architecture, Programming, and Algorithms (Bristol and Philadelphia: Adam Hilger)，§5.2.4 和 §5.4.2。[6]

5.5 复数运算

由于 C++ 内置了 complex 类，通常你可以让编译器和类库为你处理复数运算。但“通常”并不等于“总是”。如果你的程序中只包含少量复数运算，你可能希望手动内联编写这些运算。或者，你可能会发现你的编译器不够完善：经常会遇到这样的情况——复数操作数和复数结果本身都是完全可以表示的，但复数运算却产生了溢出或下溢。我们认为，这种情况的发生是因为软件公司将复数运算的实现误认为是一项完全平凡的任务，不需要任何特别的技巧。

实际上，复数运算并非完全平凡。加法和减法以显而易见的方式进行，即分别对操作数的实部和虚部执行相应运算。乘法也可以用显而易见的方式完成，需要四次乘法、一次加法和一次减法：

$$(a+ib)(c+id)=(ac-bd)+i\left(bc+ad\right)$$

（此处 $i$ 前的加号不算作运算；它仅用于在记号上分隔实部和虚部）。但有时通过下式进行乘法会更快：

$$(a+ib)(c+id)=(ac-bd)+i\left[(a+b)(c+d)-ac-bd\right]$$

该式仅需三次乘法 $(ac,bd,(a+b)(c+d))$，外加两次加法和三次减法。虽然总运算次数多了两次，但在某些机器上乘法是较慢的操作。

诚然，公式 (5.5.1) 和 (5.5.2) 中的中间结果可能在最终结果可表示的情况下发生溢出，但这种情况仅在最终结果处于可表示范围的边缘时才会发生。但对于复数模长则不然——如果你或你的编译器错误地按如下方式计算：

$$|a+ib|=\sqrt{a^2+b^2}$$

那么当中间结果 $a^2$ 或 $b^2$ 超过最大可表示数的平方根（例如 $10^{19}$ 相对于 $10^{38}$）时就会溢出。正确的计算方法应为：

$$|a+ib|=\{\begin{array}{ll}|a|\sqrt{1+(b/a{)}^2}& |a|\ge |b|\\ |b|\sqrt{1+(a/b{)}^2}& |a|<|b|\end{array}$$

复数除法也应采用类似的技巧，以避免不必要的溢出、下溢或精度损失：

$$\frac{a+ib}{c+id}=\{\begin{array}{ll}\frac{[a+b(d/c)]+i[b-a(d/c)]}{c+d(d/c)}& |c|\ge |d|\\ \frac{[a(c/d)+b]+i[b(c/d)-a]}{c(c/d)+d}& |c|<|d|\end{array}$$

当然，像 $c/d$ 或 $d/c$ 这样的重复子表达式应当只计算一次。

复数平方根更为复杂，因为我们既要保护中间结果，又要强制执行所选的分支切割（此处取负实轴）。要计算 $c+id$ 的平方根，首先计算：

$$w\equiv \{\begin{array}{ll}0& c=d=0\\ \sqrt{|c|}\sqrt{\frac{1+\sqrt{1+(d/c{)}^2}}{2}}& |c|\ge |d|\\ \sqrt{|d|}\sqrt{\frac{|c/d|+\sqrt{1+(c/d{)}^2}}{2}}& |c|<|d|\end{array}$$

然后结果为：

$$\sqrt{c+id}=\{\begin{array}{ll}0& w=0\\ w+i\left(\frac{d}{2w}\right)& w\ne 0,c\ge 0\\ \frac{|d|}{2w}+iw& w\ne 0,c<0,d\ge 0\\ \frac{|d|}{2w}-iw& w\ne 0,c<0,d<0\end{array}$$

参考文献

Midy, P., and Yakovlev, Y. 1991, “Computing Some Elementary Functions of a Complex Variable,” Mathematics and Computers in Simulation, vol. 33, pp. 33–49。

Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley) [参见习题 4.2.1.16 和 4.6.4.41 的解答]。

5.6 二次与三次方程

简单代数方程的根可以看作是方程系数的函数。我们在初等代数中就学过这些函数。然而，令人惊讶的是，许多人并不知道求解具有两个实根的二次方程或求三次方程根的正确方法。

二次方程

$$a{x}^2+bx+c=0$$

的解可以用两种方式表示，其中系数 $a,b,c$ 为实数：

$$x=\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$

和

$$x=\frac{2c}{-b\pm \sqrt{b^2-4ac}}$$

如果你使用上述任一公式（5.6.2）或（5.6.3）来计算两个根，就会遇到问题：如果 $a$ 或 $c$（或两者）很小，那么其中一个根将涉及从一个与 $b$ 几乎相等的量（判别式）中减去 $b$，从而导致该根的计算精度非常低。计算根的正确方法是：

$$q\equiv -\frac{1}{2}[b+\mathrm{sgn}(b)\sqrt{b^2-4ac}]$$

然后两个根为

$$x_1=\frac{q}{a}\mathrm{and}{x}_2=\frac{c}{q}$$

如果系数 $a,b,c$ 是复数而非实数，则上述公式仍然成立，只是在公式（5.6.4）中，平方根的符号应选择使得

$$\mathrm{Re}(b^{\ast }\sqrt{b^2-4ac})>0$$

其中 $\mathrm{Re}$ 表示实部，星号表示复共轭。

谈到二次方程，这里似乎是一个方便的地方，用来回顾反双曲函数 ${\sinh }^{-1}$ 和 ${\cosh }^{-1}$ 实际上只是这类方程解的对数形式：

$$\begin{array}{rl}{\sinh }^{-1}(x)=& \ln \left(x+\sqrt{x^2+1}\right)\end{array}$$

$${\cosh }^{-1}(x)=\pm \ln (x+\sqrt{x^2-1})$$

公式（5.6.7）在 $x\ge 0$ 时数值上是稳定的。对于负的 $x$，可利用对称性 ${\sinh }^{-1}(-x)=-{\sinh }^{-1}(x)$。公式（5.6.8）当然仅在 $x\ge 1$ 时有效。

对于三次方程

$$x^3+a{x}^2+bx+c=0$$

其中系数 $a$、$b$、$c$ 为实数或复数，首先计算

$$Q\equiv \frac{a^2-3b}{9}\mathrm{and}R\equiv \frac{2a^3-9ab+27c}{54}$$

如果 $Q$ 和 $R$ 为实数（当 $a,b,c$ 为实数时总是成立）且满足 $R^2<Q^3$，则该三次方程有三个实根。此时通过计算

$$\theta =\arccos (R/\sqrt{Q^3})$$

可得三个根为

$$\begin{array}{r}{x}_1=-2\sqrt{Q}\cos \left(\frac{\theta }{3}\right)-\frac{a}{3}\\ x_2=-2\sqrt{Q}\cos \left(\frac{\theta +2\pi }{3}\right)-\frac{a}{3}\\ x_3=-2\sqrt{Q}\cos \left(\frac{\theta -2\pi }{3}\right)-\frac{a}{3}\end{array}$$

（此公式首次出现在弗朗索瓦·韦达（François Viète）1615年出版的著作《De emendatione》第六章中！）

否则，计算

$$A=-{\left[R+\sqrt{R^2-Q^3}\right]}^{1/3}$$

其中平方根的符号应选择使得

$$\mathrm{Re}\left(R^{\ast }\sqrt{R^2-Q^3}\right)\ge 0$$

（星号再次表示复共轭）。若 $Q$ 和 $R$ 均为实数，则方程 (5.6.13)–(5.6.14) 等价于

$$A=-\mathrm{sgn}(R){\left[|R|+\sqrt{R^2-Q^3}\right]}^{1/3}$$

其中取正的平方根。接着计算

$$B=\{\begin{array}{ll}Q/A& (A\ne 0)\\ 0& (A=0)\end{array}$$

三个根可表示为

$$x_1=(A+B)-\frac{a}{3}$$

（当 $a$、$b$、$c$ 为实数时，这是唯一的实根），以及

$$\begin{array}{rl}{x}_2& =-\frac{1}{2}(A+B)-\frac{a}{3}+i\frac{\sqrt{3}}{2}(A-B)\\ x_3& =-\frac{1}{2}(A+B)-\frac{a}{3}-i\frac{\sqrt{3}}{2}(A-B)\end{array}$$

（在相同情况下，构成一对共轭复根）。方程 (5.6.13)–(5.6.16) 的安排既是为了最小化舍入误差，也是（如 A.J. Glassman 所指出的）为了确保无论怎样选择复立方根的分支，都不会导致某个不同的根被错误地丢失。

如果你需要求解许多系数仅有微小差异的三次方程，使用牛顿法（§9.4）会更高效。

参考文献

Weast, R.C. (ed.) 1967, Handbook of Tables for Mathematics, 3rd ed. (Cleveland: The Chemical Rubber Co.), pp. 130–133.

Pachner, J. 1983, Handbook of Numerical Analysis Applications (New York: McGraw-Hill), §6.1.

McKelvey, J.P. 1984, “Simple Transcendental Expressions for the Roots of Cubic Equations,” American Journal of Physics, vol. 52, pp. 269–270；另见 vol. 53, p. 775，以及 vol. 55, pp. 374–375。

5.7 数值导数

假设你有一个计算函数 $f(x)$ 的过程，现在你想计算其导数 $f^{\prime }(x)$。这很简单，对吧？导数的定义是当 $h\to 0$ 时的极限：

$$f^{\prime }(x)\approx \frac{f(x+h)-f(x)}{h}$$

这几乎直接给出了程序实现方法：选取一个很小的 $h$；计算 $f(x+h)$；你可能已经计算过 $f(x)$，如果没有，也一并计算；最后应用公式 (5.7.1)。还需要多说什么呢？

实际上，还有很多需要注意的地方。如果不加批判地直接应用上述方法，几乎肯定会得到不准确的结果。只有在以下情况下，这种方法才是合适的：函数 $f$ 的计算代价非常高；你已经投入资源计算了 $f(x)$；并且因此希望仅通过一次额外的函数求值就得到导数值。在这种情形下，剩下的问题就是如何恰当地选择 $h$，这正是我们现在要讨论的问题。

公式 (5.7.1) 中存在两种误差来源：截断误差和舍入误差。截断误差来自泰勒级数展开的高阶项：

$$f(x+h)=f(x)+h{f}^{\prime }(x)+\frac{1}{2}{h}^2{f}^{\prime \prime }(x)+\frac{1}{6}{h}^3{f}^{\prime \prime \prime }(x)+\cdots$$

由此可得

$$\frac{f(x+h)-f(x)}{h}=f^{\prime }+\frac{1}{2}h{f}^{\prime \prime }+\cdots$$

舍入误差则有多种来源。首先是 $h$ 本身的舍入误差：举个例子，假设你当前位于 $x=10.3$，并盲目地选择 $h=0.0001$。无论是 $x=10.3$ 还是 $x+h=10.30001$，在二进制中都不是能精确表示的数；因此它们都会以机器浮点格式所特有的某种相对误差 ${ϵ}_m$ 被表示，单精度下 ${ϵ}_m$ 的值约为 $10^{-7}$。于是，有效 $h$ 值（即机器中表示的 $x+h$ 与 $x$ 之差）的误差量级约为 ${ϵ}_{m}x$，

这意味着 $h$ 的相对误差量级约为 $\sim {ϵ}_{m}x/h\sim 10^{-2}$！根据公式 (5.7.1)，这立即意味着导数至少也会有同样大的相对误差。

由此我们得到第一条经验法则：始终选择 $h$，使得 $x+h$ 与 $x$ 之差是一个机器中能精确表示的数。这通常可以通过以下程序步骤实现：

$$\begin{array}{rl}\mathrm{temp}& =x+h\\ h& =\mathrm{temp}-x\end{array}$$

某些优化编译器，以及某些浮点运算单元内部精度高于外部存储精度的计算机，可能会破坏这一技巧；若出现这种情况，通常只需将 temp 声明为 volatile，或者在两个公式 (5.7.4) 之间调用一个虚拟函数 donoting(temp) 即可。这会强制 temp 进出可寻址内存。

当 $h$ 是一个“精确”数时，公式 (5.7.1) 中的舍入误差约为 $e_r\sim {ϵ}_f|f(x)/h|$。这里 ${ϵ}_f$ 是计算 $f$ 时的相对精度；对于简单函数，${ϵ}_f$ 可能与机器精度相当，即 ${ϵ}_f\approx {ϵ}_m$，但对于包含其他误差来源的复杂计算，${ϵ}_f$ 可能更大。公式 (5.7.3) 中的截断误差量级为 $e_t\sim |h{f}^{\prime \prime }(x)|$。通过调整 $h$ 使总误差 $e_r+e_t$ 最小，即可得到最优的 $h$ 选择。

$$h\sim \sqrt{\frac{{ϵ}_{f}f}{f^{\prime \prime }}}\approx \sqrt{{ϵ}_f}{x}_c$$

其中 $x_c\equiv (f/f^{\prime \prime }{)}^{1/2}$ 是函数 $f$ 的“曲率尺度”或其变化的“特征尺度”。在缺乏其他信息的情况下，通常假设 $x_c=x$（但在 $x=0$ 附近应使用其他典型 $x$ 尺度的估计值）。

采用式 (5.7.5) 的选择后，计算导数的相对精度为

$$(e_r+e_t)/|f^{\prime }|\sim \sqrt{{ϵ}_f}(f{f}^{\prime \prime }/f^{\prime 2}{)}^{1/2}\sim \sqrt{{ϵ}_f}$$

此处最后的量级近似假设 $f$、$f^{\prime }$ 和 $f^{\prime \prime }$ 具有相同的特征长度尺度，这通常是成立的。可以看出，简单的有限差分公式 (5.7.1) 最多只能达到机器精度 ${ϵ}_m$ 的平方根。

如果你能为每次导数计算负担两次函数求值，那么使用对称形式会显著更好：

$$f^{\prime }(x)\approx \frac{f(x+h)-f(x-h)}{2h}$$

在此情况下，根据式 (5.7.2)，截断误差为 $e_t\sim h^2{f}^{\prime \prime \prime }$。舍入误差 $e_r$ 与之前大致相同。通过类似于上述的简短计算，现在最优的 $h$ 选择为

$$h\sim {\left(\frac{{ϵ}_{f}f}{f^{\prime \prime \prime }}\right)}^{1/3}\sim ({ϵ}_f{)}^{1/3}{x}_c$$

相对误差为

$$(e_r+e_t)/|f^{\prime }|\sim ({ϵ}_f{)}^{2/3}{f}^{2/3}(f^{\prime \prime \prime }{)}^{1/3}/f^{\prime }\sim ({ϵ}_f{)}^{2/3}$$

这通常比式 (5.7.6) 的结果好一个数量级（单精度）或两个数量级（双精度）。由此我们得到第二条经验：选择 $h$ 为 ${ϵ}_f$ 或 ${ϵ}_m$ 的适当幂次乘以特征尺度 $x_c$。

你可以轻松推导出其他情形下的正确幂次 [1]。例如，对于二维函数和混合导数公式

$$\frac{{\partial }^{2}f}{\partial x\partial y}=\frac{[f(x+h,y+h)-f(x+h,y-h)]-[f(x-h,y+h)-f(x-h,y-h)]}{4h^2}$$

正确的缩放关系为 $h\sim {ϵ}_f^{1/4}{x}_c$。

令人失望的是，像式 (5.7.1) 或 (5.7.7) 这样的简单有限差分公式无法达到与机器精度 ${ϵ}_m$ 相当的精度，甚至无法达到函数 $f$ 本身的计算精度 ${ϵ}_f$。难道没有更好的方法吗？

答案是肯定的。然而，所有这些方法都需要在与 $x_c$ 相当的尺度上探索函数的行为，并且需要对函数的光滑性或解析性做出某些假设，以确保泰勒展开（如式 (5.7.2)）中的高阶项具有实际意义。这类方法还需要多次计算函数 $f$，因此其精度的提升必须与计算成本的增加相权衡。

“理查森延迟趋近法”（Richardson's deferred approach to the limit）的总体思想尤其吸引人。对于数值积分，该思想引出了所谓的龙贝格积分（Romberg integration）（参见 §4.3）。对于导数，我们试图对有限差分计算结果进行外推，使其在 $h\to 0$ 时收敛。通过使用 Neville 算法（§3.2），每次新的有限差分计算不仅产生更高阶的外推结果，还能对先前较低阶的结果进行更小步长 $h$ 的外推。Ridders [2] 给出了该思想的一个优秀实现；下面的程序 dfridr 基于他的算法，并通过改进的终止准则进行了修改。该例程的输入包括一个函数 $f$（称为 func）、一个位置 $x$ 和一个最大步长 $h$（这里的 $h$ 更类似于上文提到的 $x_c$，而非上文中的 $h$）。输出为导数值及其误差估计 err。

template<class T>
Doub dfridr(T &func, const Doub x, const Doub h, Doub &err)

通过 Ridders 多项式外推法返回函数 func 在点 x 处的导数。输入值 h 为估计的初始步长；它不必很小，但应是在 x 上使 func 发生显著变化的增量。导数的误差估计通过 err 返回。

const Int ntab=10;
const Doub con=1.4, con2=(con*con); // 每次迭代步长按 CON 缩小
const Doub big=numeric_limits<Doub>::max();
const Doub safe=2.0;
Int i,j;
Doub errt,fac,hh,ans;
MatDoub a(ntab,ntab);
if (h == 0.0) throw("h must be nonzero in dfridr.");
hh=h;
a[0][0]=(func(x+hh)-func(x-hh))/(2.0*hh);
err=big;
for (i=1;i<ntab;i++) {
    // Neville 表中的后续列将使用更小的步长和更高阶的外推
    hh /= con;
    a[0][i]=(func(x+hh)-func(x-hh))/(2.0*hh);
    fac=con2;
    // 尝试新的、更小的步长
    // 当误差比当前最佳结果差 SAFE 倍时返回
    for (j=1;j<=i;j++) {
        // 计算不同阶数的外推，无需新的函数求值
        a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
        fac=con2*fac;
        // 误差策略是比较每个新外推结果与其低一阶的结果，
        // 分别在当前步长和前一步长下进行
        errt=MAX(abs(a[j][i]-a[j-1][i]), abs(a[j][i]-a[j-1][i-1]));
        if (errt <= err) {
            // 如果误差减小，保存改进后的结果
            err=errt;
            ans=a[j][i];
        }
    }
    if (abs(a[i][i]-a[i-1][i-1]) >= safe*err) break;
    // 如果高阶结果比当前结果差 SAFE 倍以上，则提前终止
}
return ans;

在 dfridr 中，func 的求值次数通常为 6 到 12 次，但最多可达 2×NTAB。对于输入的 h，通常随着 h 增大，精度会提高，直到某个临界点，此时无意义的外推会导致提前返回并给出较大的误差。因此，你应该选择一个相当大的 h 值，但需监控返回的误差 err，若 err 不够小则减小 h。对于特征尺度 $x$ 为单位量级的函数，我们通常取 h 为零点几。

除了 Ridders 方法外，还有其他可能的技术。如果你的函数相当光滑，并且你知道需要在某个区间内的任意点多次计算其导数，那么在该区间内构造函数的切比雪夫多项式逼近，并直接从得到的切比雪夫系数计算导数是有意义的。该方法将在随后的 §5.8 – §5.9 中描述。

另一种技术适用于函数由等间距采样的数据组成，且可能包含噪声的情形。此时，你可能希望在每个点处，使用目标 $x$ 值左侧的 $n_L$ 个点和右侧的 $n_R$ 个点，通过最小二乘法拟合一个 $M$ 阶多项式。估计的导数即为拟合多项式的导数。实现这一构造的一种高效方法是使用 Savitzky-Golay 平滑滤波器，这将在后面的 §14.9 中讨论。在那里，我们将给出一个获取滤波器系数的例程，这些系数不仅能构造拟合多项式，还能通过数据点与滤波器系数的单次加权求和直接计算导数值。事实上，给出的例程 savgol 有一个参数 ld，用于指定计算拟合多项式的哪一阶导数。对于一阶导数，应设置 ld=1，导数值即为加权求和结果除以采样间隔 $h$。

参考文献

Dennis, J.E., and Schnabel, R.B. 1983, Numerical Methods for Unconstrained Optimization and Nonlinear Equations; reprinted 1996 (Philadelphia: S.I.A.M.), §5.4 – §5.6.[1]

Ridders, C.J.F. 1982, "Accurate computation of $F^{\prime }(x)$ and $F^{\prime }(x)F^{\prime \prime }(x)$," Advances in Engineering Software, vol. 4, no. 2, pp. 75–76.[2]

5.8 切比雪夫逼近

$n$ 阶切比雪夫多项式记作 $T_n(x)$，其显式表达式为：

$$T_n(x)=\cos (n\arccos x)$$

乍看之下，这似乎是一个三角函数表达式（事实上，切比雪夫多项式与离散傅里叶变换之间确实存在紧密联系）；然而，将公式 (5.8.1) 与三角恒等式结合，可以得到 $T_n(x)$ 的显式多项式表达式（见图 5.8.1）：

$$\begin{array}{rl}& T_0(x)=1\\ & T_1(x)=x\\ & T_2(x)=2x^2-1\\ & T_3(x)=4x^3-3x\\ & T_4(x)=8x^4-8x^2+1\\ & \cdots \\ & T_{n+1}(x)=2x{T}_n(x)-T_{n-1}(x)n\ge 1.\end{array}$$

（也存在将 $x$ 的幂次用 $T_n$ 表示的逆公式——参见例如 [1]。）

切比雪夫多项式在区间 $[-1,1]$ 上关于权函数 $(1-x^2{)}^{-1/2}$ 正交。具体而言，

$${\int }_{-1}^1\frac{T_i(x)T_j(x)}{\sqrt{1-x^2}}dx=\{\begin{array}{ll}0& i\ne j\\ \pi /2& i=j\ne 0\\ \pi & i=j=0\end{array}$$

多项式 $T_n(x)$ 在区间 $[-1,1]$ 内有 $n$ 个零点，其位置为：

$$x=\cos \left(\frac{\pi (k+\frac{1}{2})}{n}\right)k=0,1,\dots ,n-1$$

在同一区间内，它还有 $n+1$ 个极值点（极大值和极小值），其位置为：

$$x=\cos \left(\frac{\pi k}{n}\right)k=0,1,\dots ,n$$

在所有极大值点处，$T_n(x)=1$；在所有极小值点处，$T_n(x)=-1$；正是这一性质使得切比雪夫多项式在函数的多项式逼近中极为有用。

切比雪夫多项式不仅满足连续正交关系 (5.8.3)，还满足离散正交关系：若 $x_k$（$k=0,\dots ,m-1$）是由 (5.8.4) 给出的 $T_m(x)$ 的 $m$ 个零点，且 $i,j<m$，则有

$$\sum _{k=0}^{m-1}{T}_i(x_k)T_j(x_k)=\{\begin{array}{ll}0& i\ne j\\ m/2& i=j\ne 0\\ m& i=j=0\end{array}$$

图 5.8.1. 切比雪夫多项式 $T_0(x)$ 至 $T_6(x)$。注意 $T_j$ 在区间 $(-1,1)$ 内有 $j$ 个根，且所有多项式均被限制在 $\pm 1$ 之间。

结合公式 (5.8.1)、(5.8.4) 和 (5.8.6) 并不困难，可以证明以下定理：若 $f(x)$ 是区间 $[-1,1]$ 上的任意函数，且定义 $N$ 个系数 $c_j$（$j=0,\dots ,N-1$）为

$$\begin{array}{rl}& c_j=\frac{2}{N}{\sum }_{k=0}^{N-1}f(x_k)T_j(x_k)\\ & =\frac{2}{N}{\sum }_{k=0}^{N-1}f\left[\cos \left(\frac{\pi (k+\frac{1}{2})}{N}\right)\right]\cos \left(\frac{\pi j(k+\frac{1}{2})}{N}\right)\end{array}$$

则逼近公式

$$f(x)\approx \left[\sum _{k=0}^{N-1}{c}_k{T}_k(x)\right]-\frac{1}{2}{c}_0$$

在 $x$ 取 $T_N(x)$ 的全部 $N$ 个零点时是精确成立的。

对于固定的 $N$，公式 (5.8.8) 是一个关于 $x$ 的多项式，用于在区间 $[-1,1]$（即 $T_N(x)$ 所有零点所在区间）内逼近函数 $f(x)$。为何这种特定的逼近多项式优于其他在另外 $N$ 个点上精确的多项式？答案并非 (5.8.8) 一定比其他同阶 $N$ 的逼近多项式更“准确”（按某种“准确”的定义），而是 (5.8.8) 可以非常优雅地截断为一个更低阶 $m\ll N$ 的多项式，而这种截断方式恰好能给出 $m$ 阶下“最准确”的逼近（这一说法可以严格化）。假设 $N$ 足够大，使得 (5.8.8) 几乎完美地逼近了 $f(x)$。现在考虑截断后的逼近式：

$$f(x)\approx \left[\sum _{k=0}^{m-1}{c}_k{T}_k(x)\right]-\frac{1}{2}{c}_0$$

其中系数 $c_j$ 仍由 (5.8.7) 计算得到。由于所有 $T_k(x)$ 均被限制在 $\pm 1$ 之间，(5.8.9) 与 (5.8.8) 之间的误差不会超过被忽略的系数 $c_k$（$k=m,\dots ,N-1$）之和。事实上，若 $c_k$ 快速衰减（这通常是实际情况），则误差主要由 $c_m{T}_m(x)$ 主导，这是一个在区间 $[-1,1]$ 上具有 $m+1$ 个等幅极值点的振荡函数。这种误差在区间上平滑分布的特性极为重要：切比雪夫逼近 (5.8.9) 几乎等同于逼近多项式中的“圣杯”——极小极大（minimax）多项式，后者在所有同阶多项式中与真实函数 $f(x)$ 的最大偏差最小。极小极大多项式极难求得，而切比雪夫逼近多项式几乎与之相同，却非常容易计算！

因此，给定某种（可能较复杂）计算函数 $f(x)$ 的方法后，我们现在需要实现两个算法：一是计算 (5.8.7) 中的系数 $c_k$；二是在检查所得 $c_k$ 并选定截断阶数 $m$ 后，对 (5.8.9) 进行求值。此后，(5.8.9) 就成为一种高效计算 $f(x)$ 的便捷方法。

第一项任务较为直接。此处实现的 (5.8.7) 的一个推广形式是：允许逼近区间为任意两个边界 $a$ 和 $b$，而不仅限于 $[-1,1]$。这通过变量替换实现：

$$y\equiv \frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)}$$

并用关于 $y$ 的切比雪夫多项式逼近 $f(x)$。

为方便起见，我们将与切比雪夫多项式相关的一系列函数封装到一个对象中，尽管对其具体细节的讨论分散在 §5.8 至 §5.11 中：

struct Chebyshev {
    // 用于切比雪夫逼近及相关方法的对象。
    Int n, m;
    VecDoub c;
    Doub a, b;  // 逼近区间。
};

Chebyshev(Doub func(Doub), Doub aa, Doub bb, Int nn);
// 构造函数：用 nn 项在区间 [aa, bb] 上逼近函数 func。
Chebyshev(ChebyshevVecDoub &cc, Doub aa, Doub bb)
: n(cc.size()), m(n), c(cc), a(aa), b(bb) {}
// 从已计算的系数构造对象。
Int setm(Doub thresh) { while (m > 1 && abs(c[m-1]) < thresh) m--; return m; }
// 设置 m，即截断到误差阈值 thresh 后保留的系数个数，并返回该值。

Doub eval(Doub x, Int m);
inline Doub operator() (Doub x) { return eval(x, m); }
// 返回切比雪夫拟合的值，可使用存储的 m 或覆盖它。

Chebyshev derivative();
Chebyshev integral();

VecDoub polycofs(Int m);
// 参见 §5.9。
inline VecDoub polycofs() { return polycofs(m); }
// 参见 §5.11。

第一个构造函数（以任意函数 func 作为第一个参数）计算并保存 nn 个切比雪夫系数，用于在区间 [aa, bb] 上逼近 func。（你可以暂时忽略第二个构造函数，它只是从已计算的数据创建一个切比雪夫对象。）我们还要注意 setm 方法，它提供了一种快速截断切比雪夫级数的方法，实际上是从右侧删除所有绝对值小于某个阈值 thresh 的系数。

切比雪夫拟合：给定函数 func 和区间 [a, b] 的上下限，计算并保存 nn 个切比雪夫逼近系数，使得
$$\text{func}(x)\approx \left[\sum _{k=0}^{nn-1}{c}_k{T}_k(y)\right]-\frac{c_0}{2},$$
其中 $y$ 与 $x$ 通过 (5.8.10) 式关联。该例程设计用于中等较大的 n（例如 30 或 50），随后将系数数组 c 截断为较小的值 $m$，使得 $c_m$ 及其后的元素可忽略不计。

const Doub pi = 3.141592653589793;
Int k, j;
Doub fac, bpa, bma, y, sum;
VecDoub f(n);
bma = 0.5 * (b - a);
bpa = 0.5 * (b + a);
for (k = 0; k < n; k++) {
    y = cos(pi * (k + 0.5) / n);
    f[k] = func(y * bma + bpa);
}
fac = 2.0 / n;
for (j = 0; j < n; j++) {
    sum = 0.0;
    for (k = 0; k < n; k++) {
        sum += f[k] * cos(pi * j * (k + 0.5) / n);
    }
    c[j] = fac * sum;
}

如果你发现构造函数的执行时间主要消耗在计算 $N^2$ 个余弦值上，而不是在 $N$ 次函数求值上，那么你应该提前阅读 §12.3，特别是公式 (12.4.16)，其中展示了如何使用快速余弦变换方法来计算公式 (5.8.7)。

现在我们有了切比雪夫系数，如何计算逼近值呢？一种方法是利用公式 (5.8.2) 的递推关系，从 $T_0=1$、$T_1=x$ 生成 $T_k(x)$ 的值，同时累加公式 (5.8.9) 的和。更好的方法是使用 Clenshaw 递推公式（§5.4），将这两个过程同时完成。应用于切比雪夫级数 (5.8.9)，递推关系为：

$$\begin{array}{rl}& d_{m+1}\equiv d_m\equiv 0\\ & d_j=2x{d}_{j+1}-d_{j+2}+c_{j}j=m-1,m-2,\dots ,1\\ & f(x)\equiv d_0=x{d}_1-d_2+\frac{1}{2}{c}_0\end{array}$$

Doub Chebyshev::eval(Doub x, Int m)

切比雪夫求值：在点 $y=[x-(b+a)/2]/[(b-a)/2]$ 处计算切比雪夫多项式 ${\sum }_{k=0}^{m-1}{c}_k{T}_k(y)-c_0/2$，并返回结果作为函数值。

Doub d = 0.0, dd = 0.0, sv, y, y2;
Int j;
if ((x - a) * (x - b) > 0.0) throw("x not in range in Chebyshev::eval");
y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
for (j = m - 1; j > 0; j--) {
    sv = d;
    d = y2 * d - dd + c[j];
    dd = sv;
}
return y * d - dd + 0.5 * c[0]; // 最后一步不同。

eval 方法有一个参数，用于指定求值时应使用多少个前导系数 m。如果你只想使用之前通过 setm（或手动）设置的 m 值，那么你可以将切比雪夫对象作为函数对象使用。例如，

如果我们要在区间 $[-1,1]$ 上逼近一个偶函数，其展开式将只包含偶次切比雪夫多项式。此时构造一个所有奇次系数为零的切比雪夫对象是浪费的 [2]。相反，利用公式 (5.8.1) 中余弦的半角恒等式，我们得到关系式：

$$T_{2n}(x)=T_n(2x^2-1)$$

因此，我们可以通过将函数的自变量 $x$ 替换为 $2x^2-1$ 来构造一个更高效的偶函数切比雪夫对象，并在求值时同样进行替换。

奇函数的展开式将只包含奇次切比雪夫多项式。最佳做法是将其重写为函数 $f(x)/x$ 的展开式，后者只包含偶次切比雪夫多项式。这样还有一个额外的好处，即在 $x=0$ 附近能给出 $f(x)/x$ 的精确值。但不要试图通过数值计算 $f(x)/x$ 来构造级数。相反，$f(x)/x$ 的系数 $c_n^{\prime }$ 可以通过 $f(x)$ 的系数递推得到：

$$\begin{array}{rl}& c_{N+1}^{\prime }=0\\ & c_{n-1}^{\prime }=2c_n-c_{n+1}^{\prime },n=N-1,N-3,\dots \end{array}$$

公式 (5.8.13) 由公式 (5.8.2) 中的递推关系得出。

如果你坚持要计算一个奇次切比雪夫级数，高效的方法仍然是将 $x$ 替换为 $y=2x^2-1$ 作为函数的自变量。但现在，

你还必须将公式 (5.8.11) 中的最后一个公式改为：

$$f(x)=x\left[(2y-1)d_1-d_2+c_0\right]$$

并相应修改 eval 中的对应行。

5.8.1 切比雪夫与指数收敛

自从在 §1.1 中首次提到截断误差以来，我们已经看到了许多具有可调阶数 $M$ 的算法示例，其截断误差随某个量的 $M$ 次幂而减小。例子包括第 3 章中的大多数插值方法和第 4 章中的大多数数值积分方法。在这些例子中，还有另一个参数 $N$，即函数被求值的点数。

我们多次警告过：“高阶并不一定意味着更高精度。” 当 $N$ 固定而 $M$ 增加时，这仍然是很好的建议。然而，近年来科学计算的许多领域出现了一个新主题：在某些特殊情况下，可以同时增加 $M$ 和 $N$，使得误差不仅随阶数提高而减小，而且呈指数级减小！

几乎所有这些相对较新方法的共同点是一个显著的事实：当 $N$ 增加时，无限光滑的函数会被 $N$ 个采样点以指数精度确定。因此，单纯的幂律收敛可能仅仅是由于：(i) 函数不够光滑，或 (ii) 端点效应。

我们在第 4 章已经看到了几个这样的例子。在 §4.1 中，我们指出高阶数值积分规则可以在内部具有单位权重，就像梯形法则一样；所有的“高阶性”都是通过对边界附近的适当处理获得的。在 §4.5 中，我们进一步看到，将边界推至无穷远的变量变换可以产生快速收敛的数值积分算法。在 §4.5.1 中，我们实际上通过欧拉-麦克劳林公式证明了指数收敛。然后在 §4.6 中，我们提到高斯积分的收敛速度可以是指数级的（用上面的语言来说，就是同时增加 $M$ 和 $N$ 的一个例子）。

切比雪夫逼近可以因不同（但相关）的原因而实现指数收敛：光滑的周期函数根本没有端点，从而避免了端点效应！切比雪夫逼近可以看作是将 $x$ 区间 $[-1,1]$ 映射到角区间 $[0,\pi ]$（参见公式 5.8.4 和 5.8.5），使得区间 $[-1,1]$ 上的任何无限光滑函数都变成 $[0,2\pi ]$ 上的无限光滑、偶的周期函数。图 5.8.2 从几何上展示了这一思想。通过将横坐标投影到半圆上，产生半个周期。另一半周期通过对称反射得到，或者可以想象为将函数投影到一个相同的下半圆上。切比雪夫多项式的零点（或切比雪夫逼近的节点）在圆周上等距分布，而切比雪夫多项式本身就是一个余弦函数（参见公式 5.8.1）。这说明了切比雪夫逼近与圆周上周期函数之间的紧密联系；在第 12 章中，我们将以几乎等价的方式对这类函数应用离散傅里叶变换（§12.4.2）。

切比雪夫方法之所以效果如此之好（高斯求积法效果好也是同理），其根本原因在于采样点在区间端点附近以一种特殊的方式聚集。

图 5.8.2. 几何构造图，说明切比雪夫逼近与周期函数的关系。
(a) 中绘制了区间上的一个光滑函数；
(b) 中将横坐标映射到一个半圆上；
(c) 中将半圆展开。由于半圆在端点处具有垂直切线，函数在端点附近几乎为常数。事实上，若将其延拓到区间 $[\pi ,2\pi ]$ 上，则它在 $[0,2\pi ]$ 上成为一个光滑、偶对称的周期函数。

任何在区间上有界的函数，即使在复平面上存在邻近的极点，其切比雪夫逼近在 $N\to \infty$ 时仍会收敛。对于非无限光滑的函数，实际收敛速率取决于函数的光滑程度：有界导数的阶数越高，收敛速率越快。特别地，对于 $C^{\infty }$ 函数，收敛是指数级的。在 §3.0 中讨论多项式插值时，我们提到了硬币的另一面：区间上等间距的采样点几乎是几何上最差的选择，通常会导致病态问题。

采样定理（§4.5、§6.9、§12.1、§13.11）的使用通常与指数收敛方法密切相关。在 §20.7 讨论偏微分方程的谱方法时，我们将再次回到指数收敛方法的诸多概念。

参考文献

Arfken, G. 1970, Mathematical Methods for Physicists, 2nd ed. (New York: Academic Press), p. 631.[1]

Clenshaw, C.W. 1962, Mathematical Tables, vol. 5, National Physical Laboratory (London: H.M. Stationery Office).[2]

Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), Chapter 8.

Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall); reprinted 2003 (New York: Dover), §4.4.1, p. 104.

Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: Addison-Wesley), §6.5.2, p. 334.

Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), §1.10, p. 39.

5.9 切比雪夫逼近函数的导数或积分

如果你已经获得了在某个区间上逼近某函数的切比雪夫系数（例如通过 §5.8 中的 chebfft 得到），那么将这些系数转换为该函数导数或积分所对应的切比雪夫系数是一件很简单的事情。完成转换后，你就可以像处理一个从头开始用切比雪夫拟合的函数一样，对导数或积分进行求值。

相关公式如下：设 $c_i$（$i=0,\dots ,m-1$）是用于逼近函数 $f$ 的系数（见公式 (5.8.9)），$C_i$ 是用于逼近 $f$ 的不定积分的系数，$c_i^{\prime }$ 是用于逼近 $f$ 的导数的系数，则有：

$$\begin{array}{rl}& C_i=\frac{c_{i-1}-c_{i+1}}{2i}(i>0)\\ & c_{i-1}^{\prime }=c_{i+1}^{\prime }+2i{c}_i(i=m-1,m-2,\dots ,1)\end{array}$$

公式 (5.9.1) 需要任意指定 $C_0$，这对应于积分常数的任意选择。公式 (5.9.2) 是一个递推关系，其初始值设为 $c_m^{\prime }=c_{m-1}^{\prime }=0$，这表示我们对原函数 $f$ 的第 $m+1$ 阶切比雪夫系数没有任何信息。

以下是实现公式 (5.9.1) 和 (5.9.2) 的例程。每个例程都返回一个新的切比雪夫对象，你可以对其调用 setm、eval，或直接作为函数对象使用。

返回一个新的切比雪夫对象，用于在相同区间 [a,b] 上逼近原函数的导数。

{
    Int j;
    Doub con;
    VecDoub cder(n);
    cder[n-1] = 0.0;
    cder[n-2] = 2*(n-1)*c[n-1];
    for (j = n-2; j > 0; j--)
        cder[j-1] = cder[j+1] + 2*j*c[j];
    con = 2.0/(b-a);
    for (j = 0; j < n; j++) cder[j] *= con;
    return Chebyshev(cder, a, b);
}

返回一个新的切比雪夫对象，用于在相同区间 [a,b] 上逼近原函数的不定积分。积分常数被设定为使得积分在 $a$ 处为零。

sum += fac * cint[n-1];
cint[0] = 2.0 * sum;
return Chebyshev(cint, a, b);

设置积分常数

5.9.1 克伦肖-柯蒂斯求积法（Clenshaw-Curtis Quadrature）

由于光滑函数的切比雪夫系数 $c_i$ 通常快速衰减（一般呈指数衰减），公式 (5.9.1) 常常可作为高效求积方案的基础。如上所述，切比雪夫对象可用于计算积分 ${\int }_a^{x}f(x)dx$，尤其当需要在区间 $a\le x\le b$ 内多个不同的 $x$ 值处求值时。如果仅需计算单个定积分 ${\int }_a^{b}f(x)dx$，则可使用由公式 (5.9.1) 推导出的更简单公式：

$${\int }_a^{b}f(x)dx=(b-a)\left[\frac{1}{2}{c}_0-\frac{1}{3}{c}_2-\frac{1}{15}{c}_4-\cdots -\frac{1}{(2k+1)(2k-1)}{c}_{2k}-\cdots \right]$$

其中 $c_i$ 由 chebft 返回。当 $c_{2k}$ 变得可忽略时，级数即可截断，第一个被忽略的项可作为误差估计。

该方法称为克伦肖-柯蒂斯求积法 [1]。它通常与自适应选择 $N$ 相结合，其中 $N$ 是通过公式 (5.8.7) 计算的切比雪夫系数个数，同时也是函数 $f(x)$ 的求值次数。如果较小的 $N$ 无法使公式 (5.9.3) 中的 $c_{2k}$ 足够小，则尝试更大的 $N$。在这种自适应情形下，用所谓的“梯形”或高斯-洛巴托（Gauss-Lobatto，见 §4.6）变体替代公式 (5.8.7) 更为高效：

$$c_j=\frac{2}{N}\sum _{k=0}^N{}^{\prime \prime }f\left[\cos \left(\frac{\pi k}{N}\right)\right]\cos \left(\frac{\pi jk}{N}\right),j=0,\dots ,N-1$$

注意（特别注意！）：求和符号上的两个撇号表示求和中的首项和末项需乘以 $1/2$。若在公式 (5.9.4) 中将 $N$ 加倍，则新函数求值点中有一半与旧点重合，从而可以复用之前的函数求值结果。这一特性，加上解析形式的权重和横坐标（即公式 (5.9.4) 中的余弦函数），通常使克伦肖-柯蒂斯求积法优于高阶自适应高斯求积法（参见 §4.6.4），尽管两者在其他方面非常相似。

如果你的问题迫使你使用较大的 $N$ 值，需要注意的是，公式 (5.9.4) 可通过快速余弦变换（fast cosine transform）对所有 $j$ 值同时快速求值。（参见 §12.3，特别是公式 12.4.11。我们之前已指出，非梯形形式 (5.8.7) 也可通过快速余弦方法计算，参见公式 12.4.16。）

参考文献

Goodwin, E.T.（编）1961，《现代计算方法》，第2版（纽约：Philosophical Library），第78–79页。

Clenshaw, C.W. 与 Curtis, A.R. 1960，“自动计算机上的数值积分方法”，《Numerische Mathematik》，第2卷，第197–205页。[1]

5.10 由切比雪夫系数构造多项式近似

在阅读完前两节后，你可能会问：我是否必须将切比雪夫近似存储并求值为一个关于变换变量 $y$ 的切比雪夫系数数组？

难道不能将这些 $c_k$ 转换为原始变量 $x$ 下的实际多项式系数，从而得到如下形式的近似？

$$f(x)\approx \sum _{k=0}^{m-1}{g}_k{x}^k,a\le x\le b$$

是的，你可以这样做（我们也会给出实现该转换的算法），但我们建议你谨慎行事：当系数 $g_k$ 源于切比雪夫近似时，直接求值公式 (5.10.1) 通常比直接求值切比雪夫级数（例如通过 eval 函数）需要更多的有效数字。这是因为切比雪夫多项式本身表现出一种相当微妙的抵消效应：例如，$T_n(x)$ 的首项系数为 $2^{n-1}$，而其他系数甚至更大；然而它们却能巧妙地组合成一个值域始终在 $\pm 1$ 之间的多项式。只有当 $m$ 不超过 7 或 8 时，你才应考虑将切比雪夫拟合直接写成多项式形式；即便如此，你也应接受其精度比机器舍入误差极限低大约两三位有效数字。

要得到公式 (5.10.1) 中的系数 $g_k$，需分两步进行。首先，使用 Chebyshev 类中的成员函数 polycofs，输出一组与已存储的 $c_k$ 等价的多项式系数（即将区间 $[a,b]$ 缩放至 $[-1,1]$）。其次，使用子程序 pcshift 将这些系数变换回去，使区间重新映射回 $[a,b]$。这两个所需子程序如下所示：

//chebyshev.h
VecDoub Chebyshev::polycofs(Int m)

由切比雪夫拟合得到多项式系数。给定系数数组 c[0..n-1]，该子程序返回系数数组 d[0..n-1]，使得
$$\sum _{k=0}^{n-1}{\mathrm{d}}_k{y}^k=\sum _{k=0}^{n-1}{c}_k{T}_k(y)-c_0/2。$$
该方法基于 Clenshaw 递推公式 (5.8.11)，但此处是代数应用而非数值计算。

Int k,j;
Doub sv;
VecDoub d(m),dd(m);
for (j=0;j<m;j++) d[j]=dd[j]=0.0;
d[0]=c[m-1];
for (j=m-2;j>0;j--) {
    for (k=m-j;k>0;k--) {
        sv=d[k];
        d[k]=2.0*d[k-1]-dd[k];
        dd[k]=sv;
    }
    sv=d[0];
    d[0] = -dd[0]+c[j];
    dd[0]=sv;
}
for (j=m-1;j>0;j--) d[j]=d[j-1]-dd[j];
d[0] = -dd[0]+0.5*c[0];
return d;

void pcshft(Doub a, Doub b, VecDoub_IO &d)

多项式系数平移。给定系数数组 d[0..n-1]，该子程序生成系数数组 $g$[0..n-1]，使得
$$\sum _{k=0}^{n-1}{\mathrm{d}}_k{y}^k=\sum _{k=0}^{n-1}{g}_k{x}^k，$$
其中 $x$ 与 $y$ 通过公式 (5.8.10) 相关联，即将区间 $-1<y<1$ 映射到区间 $a<x<b$。结果数组 $g$ 通过 d 返回。

Int k,j,n=d.size();
Doub cnst=2.0/(b-a), fac=cnst;
for (j=1;j<n;j++) {
    d[j] *= fac;
    fac *= cnst;
}
cnst=0.5*(a+b);
// 接下来通过综合除法（synthetic division）实现平移，这是高中代数中的一个奇妙技巧。
for (j=0; j<=n-2; j++)
    for (k=n-2; k>=j; k--)
        d[k] += cnst * d[k+1];

参考文献

Acton, F.S. 1970，《行之有效的数值方法》；1990年修订版（华盛顿特区：美国数学协会），第59、182–183页 [综合除法]。

5.11 幂级数的经济化

切比雪夫方法的一个特定应用——幂级数的经济化（economization of power series）——是一种偶尔非常有用的技巧，颇有“无中生有”之妙。

假设你知道如何利用一个收敛的幂级数来计算某个函数，例如：

$$f(x)\equiv \frac{1}{2}-\frac{x}{4}+\frac{x^2}{8}-\frac{x^3}{16}+\cdots$$

（这个函数实际上是 $1/(x+2)$，但暂且假装你不知道这一点。）你可能正在处理一个问题，需要在某个特定区间（例如 [0, 1]）内多次求值该级数。一切都没问题，只是该级数在误差（可用第一个被忽略的项来近似）达到可接受水平之前，需要非常多的项。在我们的例子中，当 $x=1$ 时，大约需要 30 项才能使第一个被忽略的项小于 $10^{-9}$。

注意，由于 $x^{30}$ 中的高次幂，该误差在整个区间内（除了 $x$ 最大的端点附近）远小于 $10^{-9}$。正是这一特性使得“经济化”成为可能：如果我们愿意让区间内其他位置的误差上升到与区间端点处第一个被忽略项的大小相当，那么就可以用一个显著更短的级数来替代原来的 30 项级数。

具体步骤如下：

1. 计算足够多的幂级数系数，以确保在所关注区间内函数值足够精确。

2. 按照公式 (5.8.10) 将变量从 $x$ 变换为 $y$，将 $x$ 的区间映射到 $-1\le y\le 1$。

3. 找到一个切比雪夫级数（类似于公式 5.8.8），使其精确等于你截断后的幂级数。

4. 将该切比雪夫级数截断为更少的项数，并以第一个被忽略的切比雪夫多项式的系数作为误差的估计。

5. 转换回关于 $\nu$ 的多项式。

6. 将变量换回 $x$。

除了步骤 2 和 3 之外，我们已经具备了所有步骤所需的工具。步骤 2 正好是例程 pcsht（§5.10）的逆过程，该例程将一个定义在 $y$（区间 $[-1,1]$）上的多项式映射到 $x$（区间 $[a,b]$）上。但由于方程 (5.8.10) 给出了 $x$ 与 $y$ 之间的线性关系，因此也可以用 pcsht 来实现其逆变换。具体来说，

$$\mathrm{pcsht}(a,b,\mathrm{d},\mathrm{n})$$

的逆变换为（你可以自行验证）：

void ipcsht(Doub a, Doub b, VecDoub_IO &d) {
    pcsht(-(2.+b+a)/(b-a), (2.-b-a)/(b-a), d);
}

步骤 3 需要一个新的切比雪夫构造函数，它能从一个多项式系数向量计算出切比雪夫系数。以下代码实现了这一功能。该算法基于 §5.3 中的技术，从最高次项系数 $\mathrm{d}[\mathrm{n}-1]$ 开始构造多项式，并使用方程 (5.8.2) 的递推关系，写成如下形式：

$$\begin{array}{rl}& x{T}_0=T_1\\ & x{T}_n=\frac{1}{2}(T_{n+1}+T_{n-1}),n\ge 1.\end{array}$$

唯一的微妙之处在于：由于 $T_0$ 的系数在方程 (5.8.8) 中会以因子 $1/2$ 使用，因此需要将其乘以 2。

Chebyshev::Chebyshev(VecDoub &d)
    : n(d.size()), m(n), c(n), a(-1.), b(1.)
// Chebyshev 中例程 polycofs 的逆过程：
// 给定一个多项式系数数组 d[0..n-1]，构造一个等价的切比雪夫对象。
{
    c[n-1] = d[n-1];
    c[n-2] = 2.0 * d[n-2];
    for (Int j = n-3; j >= 0; j--) {
        c[j] = 2.0 * d[j] + c[j+2];
        for (Int i = j+1; i < n-2; i++) {
            c[i] = (c[i] + c[i+2]) / 2;
        }
        c[n-2] += 2;
        c[n-1] += 2;
    }
}

将上述所有步骤组合起来，步骤 2 到 6 的代码大致如下（从一个幂级数系数向量 powser 开始）：

ipcsht(a, b, powser);
Chebyshev cpower(powser);
cpower.setm(1.e-9);
VecDoub d = cpower.polycofs();
pcshft(a, b, d);

顺便提一下，在我们的例子中，为达到 $10^{-9}$ 的精度所需项数从 30 项减少到了 9 项。在不损失任何精度的前提下，用一个 9 项多项式替代一个 30 项多项式——这看起来像是无中生有。这种技术中是否有什么魔法？其实并没有。那个 30 项多项式定义了一个函数 $f(x)$。与其对级数进行经济化处理，我们也可以直接在足够多的点上计算 $f(x)$ 的值，然后利用 §5.8 中的方法在感兴趣的区间内构造其切比雪夫逼近，结果会得到完全相同的低阶多项式。这里的主要启示是：切比雪夫系数的收敛速率与幂级数系数的收敛速率毫无关系；真正决定多项式逼近所需项数的是前者。一个函数在某个感兴趣的区域内可能具有发散的幂级数，但只要函数本身性质良好，就仍然存在非常好的多项式逼近。这些逼近可以通过 §5.8 中的方法找到，但无法通过级数经济化得到。级数经济化的效果其实并没有表面上看起来那么神奇。

参考文献

Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington, DC: Mathematical Association of America), Chapter 12.

5.12 帕德逼近（Padé Approximants）

所谓帕德逼近（Padé approximant），是指一个（指定阶数的）有理函数，其幂级数展开与给定的幂级数在尽可能高的阶数上保持一致。若该有理函数为

$$R(x)\equiv \frac{{\sum }_{k=0}^M{a}_k{x}^k}{1+{\sum }_{k=1}^N{b}_k{x}^k}$$

则称 $R(x)$ 是幂级数

$$f(x)\equiv \sum _{k=0}^{\infty }{c}_k{x}^k$$

的一个帕德逼近，如果满足

$$R(0)=f(0)$$

并且

$${\frac{d^k}{d{x}^k}R(x)|}_{x=0}={\frac{d^k}{d{x}^k}f(x)|}_{x=0},k=1,2,\dots ,M+N.$$

方程 (5.12.3) 和 (5.12.4) 为未知量 $a_0,\dots ,a_M$ 和 $b_1,\dots ,b_N$ 提供了 $M+N+1$ 个方程。理解这些方程最简单的方式是令 (5.12.1) 与 (5.12.2) 相等，两边同乘以 (5.12.1) 的分母，然后令所有含 $a$ 或 $b$ 的 $x$ 幂次项系数相等。若仅考虑对角有理逼近（即 $M=N$，参见 §3.4）这一特殊情况，则有 $a_0=c_0$，其余的 $a$ 和 $b$ 满足

$$\begin{array}{rl}& \sum _{m=1}^N{b}_m{c}_{N-m+k}=-c_{N+k},k=1,\dots ,N\\ & \sum _{m=0}^k{b}_m{c}_{k-m}=a_k,k=1,\dots ,N\end{array}$$

（注意在方程 5.12.1 中，$b_0=1$）。求解时，首先处理方程 (5.12.5)，这是一组关于所有未知 $b$ 的线性方程组。尽管该方程组具有 Toeplitz 矩阵的形式（参见方程 2.8.8），但经验表明这些方程常常接近奇异，因此不应使用 §2.8 中的方法求解，而应采用完整的 $LU$ 分解。此外，最好通过迭代改进（参见 §2.5 中的 improve 方法）[1] 对解进行精化。

一旦求得 $b$，方程 (5.12.6) 就给出了未知 $a$ 的显式公式，从而完成整个求解过程。

帕德逼近通常用于存在某个未知底层函数 $f(x)$ 的情形。我们假设你能够以某种方式（也许是通过繁琐的解析展开）计算出 $f(x)$ 及其在 $x=0$ 处的若干阶导数值：$f(0)$、$f^{\prime }(0)$、$f^{\prime \prime }(0)$ 等等。这些当然就是 $f(x)$ 幂级数展开的前几项系数；但这些系数未必趋于零，你也不知道（甚至不确定）该幂级数在何处收敛。

与切比雪夫逼近（§5.8）或幂级数经济化（§5.11）等仅能压缩你已知函数信息的技术不同，帕德（Padé）逼近能为你提供关于函数值的真正新信息。

图 5.12.1. 某函数 $f(x)$ 的五项幂级数展开及其导出的五系数帕德逼近。完整的幂级数仅在 $x<1$ 时收敛。注意，帕德逼近在远超出级数收敛半径的区域仍保持高精度。

有时这种效果好得令人费解。（与其他数学中的谜题一样，这与解析性有关。）下面的例子将说明这一点。

假设你通过非凡的努力，计算出了某个未知函数 $f(x)$ 的幂级数展开的前五项：

$$f(x)\approx 2+\frac{1}{9}x+\frac{1}{81}{x}^2-\frac{49}{8748}{x}^3+\frac{175}{78732}{x}^4+\cdots$$

（实际上，并不需要这些系数以精确的有理数形式给出——数值形式同样有效。此处我们以有理数形式书写，是为了让你感觉它们来自某种解析推导。）式 (5.12.7) 在图 5.12.1 中以标为“幂级数”的曲线绘出。可以看出，当 $x\gtrsim 4$ 时，该近似主要受其最高次（四次）项主导。

现在，我们将式 (5.12.7) 中的五个系数输入下面列出的 pade 例程。该例程返回五个有理系数：三个 $a$ 和两个 $b$，用于式 (5.12.1)，其中 $M=N=2$。图中标为“Padé”的曲线即为所得有理函数的图像。注意，两条实线曲线均源自相同的五个原始系数。

为了评估结果，我们需要“天降神助”（Deus ex machina，当他可用时，确实很有用）告诉我们，式 (5.12.7) 实际上是如下函数的幂级数展开：

$$f(x)=[7+(1+x{)}^{4/3}{]}^{1/3}$$

该函数在图中以虚线绘出。此函数在 $x=-1$ 处有一个分支点，因此其幂级数仅在区间 $-1<x<1$ 内收敛。在图中所示的大部分区间内，该级数是发散的，截断至五项的结果几乎毫无意义。然而，这五项经转换为帕德逼近后，却能在至少 $x\sim 10$ 的范围内对函数给出令人惊讶的良好近似。

为什么会这样？难道不存在其他函数，它们的幂级数前五项与此相同，但在区间（例如）$2<x<10$ 内的行为却截然不同吗？确实存在这样的函数。帕德逼近具有一种不可思议的能力，能从所有可能中挑选出你心中所想的那个函数。当然，有时它也会失败！这正是帕德逼近的缺点：它是不可控的。通常无法判断其精度有多高，也无法确定它在多大的 $x$ 范围内仍有效。这是一种强大但终究仍显神秘的技术。

以下是返回一个 Ratfn 有理函数对象的例程，该对象是你提供的幂级数系数的帕德逼近。注意，该例程专用于 $M=N$ 的情形。之后，你可以直接将 Ratfn 对象作为函数对象使用，也可以手动读取其系数（见 §5.1）。

Ratfn pade(VecDoub_I &cof)

给定 cof[0..2*n]，即某函数幂级数展开的前若干项，求解线性帕德方程，返回一个 Ratfn 对象，该对象表示对同一函数的对角有理函数逼近。

const Doub BIG=1.0e99;
Int j,k=n=(cof.size() -1)/2;
Doub sum;
MatDoub q(n,n),qlu(n,n);
VecInt idx(n);
VecDoub x(n),y(n),num(n+1),denom(n+1);
for (j=0;j<n;j++) {
    y[j]=cof[n+j+1];
    for (k=0;k<n;k++) q[j][k]=cof[j-k+n];
}
LUdcmp lu(q);
lu.solve(y,x);
for (j=0;j<4;j++) lu.mprove(y,x);
for (k=0;k<n;k++) {
    for (sum=cof[k+1],j=0;j<=k;j++) sum -= x[j]*cof[k-j];
    y[k]=sum;
}
num[0] = cof[0];
denom[0] = 1;
for (j=0;j<n;j++) {
    num[j+1]=y[j];
    denom[j+1] = -x[j];
}
return Ratfn(num,denom);

参考文献

Ralston, A. and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley), p. 14.

Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: North-Holland), Chapter 2.

Graves-Morris, P.R. 1979, in Padé Approximation and Its Applications, Lecture Notes in Mathematics, vol. 765, L. Wuytack, ed. (Berlin: Springer).[1]

5.13 有理切比雪夫逼近

在 §5.8 和 §5.10 中，我们学习了如何在给定区间 $a\le x\le b$ 内为给定函数 $f(x)$ 找到良好的多项式逼近。现在，我们希望将这一任务推广为寻找良好的有理函数逼近（见 §5.1）。这样做的原因是，对于某些函数和某些区间，最优的有理函数逼近在使用相同数量系数的情况下，能够达到远高于最优多项式逼近的精度。但这也需要权衡：寻找有理函数逼近不像寻找多项式逼近那样直接，后者如我们所见，可通过切比雪夫多项式优雅地实现。

设所求有理函数 $R(x)$ 的分子次数为 $m$，分母次数为 $k$。则有

$$R(x)\equiv \frac{p_0+p_{1}x+\cdots +p_m{x}^m}{1+q_{1}x+\cdots +q_k{x}^k}\approx f(x)\mathrm{for}a\le x\le b$$

我们需要确定的未知量是 $p_0,\dots ,p_m$ 和 $q_1,\dots ,q_k$，共计 $m+k+1$ 个量。令 $r(x)$ 表示 $R(x)$ 与 $f(x)$ 的偏差，并令 $r$ 表示其最大绝对值：

$$r(x)\equiv R(x)-f(x)r\equiv \underset{a<x<b}{\max }|r(x)|$$

理想的极小极大（minimax）解即为使 $r$ 最小化的 $p$ 和 $q$ 的取值。显然存在某个极小极大解，因为 $r$ 的下界为零。我们如何找到它，或找到一个合理的近似解？

第一个提示来自以下基本定理：若 $R(x)$ 是非退化的（即分子与分母无公共多项式因子），则存在唯一一组 $p$ 和 $q$ 使 $r$ 最小；对于这组解，$r(x)$ 在区间 $a\le x\le b$ 内具有 $m+k+2$ 个极值点，其绝对值均为 $r$，且符号交替。（该定理省略了一些技术性假设，精确表述见 Ralston [1]。）由此可知，有理函数的情形与极小极大多项式非常类似：在 §5.8 中我们看到，一个 $n$ 阶逼近（含 $n+1$ 个切比雪夫系数）的误差项通常由第一个被忽略的切比雪夫项 $T_{n+1}$ 主导，而 $T_{n+1}$ 本身具有 $n+2$ 个等幅且符号交替的极值点。因此，此处有理系数的总数 $m+k+1$ 扮演了多项式系数总数 $n+1$ 的角色。

另一种理解为何 $r(x)$ 应有 $m+k+2$ 个极值的方式是注意到：$R(x)$ 可在任意 $m+k+1$ 个点 $x_i$ 上精确等于 $f(x)$。将式 (5.13.1) 两边乘以其分母，可得方程：

$$p_0+p_1{x}_i+\cdots +p_m{x}_i^m=f(x_i)(1+q_1{x}_i+\cdots +q_k{x}_i^k)i=0,1,\dots ,m+k$$

这是一组关于未知数 $p$ 和 $q$ 的 $m+k+1$ 个线性方程，可以通过标准方法（例如 $LU$ 分解）求解。如果我们选择所有 $x_i$ 都位于区间 $(a,b)$ 内，则通常在每对相邻点 $x_i$ 与 $x_{i+1}$ 之间都会存在一个极值点，此外在函数离开区间端点 $a$ 和 $b$ 处也会出现极值点，总共会有 $m+k+2$ 个极值点。对于任意选取的 $x_i$，这些极值点的幅值通常并不相等。该定理指出，存在一种特殊的 $x_i$ 选择方式，使得所有极值的幅值都能被压低到相同的最小值 $r$。

除了在点 $x_i$ 处令 $f(x_i)$ 与 $R(x_i)$ 相等之外，还可以通过求解以下线性方程，将残差 $r(x_i)$ 强制设为任意指定值 $y_i$：

$$p_0+p_1{x}_i+\cdots +p_m{x}_i^m=[f(x_i)-y_i](1+q_1{x}_i+\cdots +q_k{x}_i^k)i=0,1,\dots ,m+k$$

事实上，如果选择 $x_i$ 为 minimax 有理逼近解的极值点（而非零点），则满足的方程为：

$$p_0+p_1{x}_i+\cdots +p_m{x}_i^m=[f(x_i)\pm r](1+q_1{x}_i+\cdots +q_k{x}_i^k)i=0,1,\dots ,m+k+1$$

其中 $\pm$ 符号在交替的极值点处交替变化。注意，方程 (5.13.5) 在 $m+k+2$ 个极值点处成立，而方程 (5.13.4) 仅在 $m+k+1$ 个任意点处成立。这是如何可能的呢？答案在于方程 (5.13.5) 中的 $r$ 是一个额外的未知量，因此方程个数和未知量个数均为 $m+k+2$。诚然，该方程组对 $r$ 略微非线性，但通常仍可通过第 9 章将要介绍的方法完美求解。

由此可见，只要已知 minimax 有理函数极值点的位置，我们就能求出其系数和最大偏差。进一步的定理（最终导向所谓的 Remes 算法 [1]）说明了如何通过迭代过程收敛到这些极值点位置。例如，以下是 Remes 第二算法的一个（稍作简化）描述：(1) 找到一个初始有理函数，使其具有 $m+k+2$ 个极值点 $x_i$（偏差不相等）；(2) 求解方程 (5.13.5)，得到新的有理系数和 $r$；(3) 计算所得 $R(x)$，找出其实际极值点（通常与猜测值不同）；(4) 将每个猜测值替换为符号相同且最近的实际极值点；(5) 返回步骤 (2) 并迭代直至收敛。在广泛的假设条件下，该方法能够收敛。Ralston [1] 补充了必要的细节，包括如何找到初始的 $x_i$ 集合。

到目前为止，我们的讨论都是教科书级别的标准内容。现在我们要暴露自己的“异端”立场了。我们并不太喜欢优雅的 Remes 算法。其双重嵌套迭代（在非线性方程组 (5.13.5) 中对 $r$ 的迭代，以及对新 $x_i$ 集合的迭代）非常挑剔，且在退化情形下需要大量特殊逻辑处理。更“异端”的是，我们怀疑执着地寻找精确最优、等偏差逼近是否值得付出如此努力——除非你是世界上极少数专门从事为编译器和微码寻找最优逼近的人。

当我们使用有理函数逼近时，目标通常更为务实：在某个内层循环中，我们需要对某个函数进行海量次求值，希望加速其计算过程。我们几乎从不需要该函数达到机器精度的最后一位。假设（异端！）我们使用一个逼近函数，其误差具有 $m+k+2$ 个极值点，且最大偏差与最小偏差相差不超过 2 倍。Remes 算法所依据的定理保证，完美的 minimax 解的极值必然落在这个 2 倍范围内——压低较大的极值会导致较小的极值上升，直到所有极值相等。因此，我们的“粗糙”逼近实际上与 minimax 解的差距不超过最低有效位的一个小数部分。

对我们而言，这已经足够好了，尤其是当我们拥有一种非常稳健的方法来寻找这种所谓的“粗糙”逼近时。这种方法就是通过奇异值分解（SVD）求解超定线性方程组的最小二乘解（参见 §2.6 和 §15.4）。具体步骤如下：首先，以最小二乘意义求解方程 (5.13.3)，但不是仅在 $m+k+1$ 个 $x_i$ 点上，而是在显著更多的点上进行，这些点大致按照高阶切比雪夫多项式的零点分布。这给出了 $R(x)$ 的初始猜测。其次，计算所得偏差，求出平均绝对偏差并记为 $r$，然后再次以最小二乘意义求解方程 (5.13.5)，其中 $r$ 固定，$\pm$ 取每个点 $x_i$ 处观测偏差的符号。第三，重复第二步若干次。

你可以发现我们的算法中仍潜藏着一些 Remes 正统思想：我们求解的方程试图将偏差调整为某个一致的正值或负值，而非零。然而，我们不再追踪实际极值点，并且每一步仅需求解线性方程组。一个额外的技巧是求解加权最小二乘问题，其中权重的选择旨在最快地压制最大偏差。

下面是一个实现上述思想的函数。注意，对函数 fn 的调用仅出现在初始填充表 fs 时。你可以轻松修改代码，将此填充过程移到例程外部。甚至你的横坐标 xs 也不必严格使用我们所选的点，尽管如果你在（底层）minimax 解的每两个极值点之间没有设置若干横坐标，拟合质量会下降。该函数返回一个 Ratfn 对象，你可以将其作为函数对象使用，或从中提取存储的系数。

Ratfn ratlsq(Doub fn(const Doub), const Doub a, const Doub b, const Int mm, const Int kk, Doub &dev)

ratlsq.h

在区间 $(a,b)$ 内返回函数 fn 的有理函数逼近。输入参数 mm 和 kk 分别指定分子和分母的阶数。逼近的最大绝对偏差（就已知而言）通过 dev 返回。

{
    const Int NPFAC=8, MAXIT=5;
    const Doub BIG=1.0e99, PIO2=1.570796326794896619;

    Int i, it, j, ncof=mm+kk+1, npt=NPFAC*ncof;
    // 函数求值的点数，即网格的精细程度。
    Doub devmax, e, hth, power, sum;
    VecDoub bb(npt), coeff(ncof), ee(npt), fs(npt), wt(npt), xs(npt);
    MatDoub u(npt, ncof);
    Ratfn ratbest(coeff, mm+1, kk+1);
    dev = BIG;

    for (i = 0; i < npt; i++) {
        if (i < (npt/2) - 1) {
            hth = PIO2 * i / (npt - 1.0);
            xs[i] = a + (b - a) * SQR(sin(hth));
        } else {
            hth = PIO2 * (npt - i) / (npt - 1.0);
            xs[i] = b - (b - a) * SQR(sin(hth));
        }
        fs[i] = fn(xs[i]);
        wt[i] = 1.0;
        ee[i] = 1.0;
    }

    e = 0.0;
    for (it = 0; it < MAXIT; it++) {
        for (i = 0; i < npt; i++) {
            power = wt[i];
            bb[i] = power * (fs[i] + SIGN(e, ee[i]));
            // 关键思想：在偏差为正处拟合 fn(x)+e，在偏差为负处拟合 fn(x)-e。
            // 这样 e 就会逼近等波纹偏差。
            for (j = 0; j < mm + 1; j++) {
                u[i][j] = power;
                power *= xs[i];
            }
            power = -bb[i];
            for (j = mm + 1; j < ncof; j++) {
                power *= xs[i];
                u[i][j] = power;
            }
        }

        SVD svd(u);
        svd.solve(bb, coeff);
        // 在特别奇异或困难的情况下，可以在此处编辑奇异值，
        // 将 w[0..ncof-1] 中的小值置零。

        devmax = sum = 0.0;
        Ratfn rat(coeff, mm + 1, kk + 1);
        for (j = 0; j < npt; j++) {
            ee[j] = rat(xs[j]) - fs[j];
            wt[j] = abs(ee[j]);
            sum += wt[j];
            if (wt[j] > devmax) devmax = wt[j];
        }
        e = sum / npt;
        if (devmax <= dev) {
            ratbest = rat;
            dev = devmax;
        }
        cout << "ratlsq iteration " << it;
        cout << " max error=" << setw(10) << devmax << endl;
    }
    return ratbest;
}

图 5.13.1 显示了当 ratlsq 用于在区间 $(0,\pi )$ 内寻找函数 $f(x)=\cos x/(1+e^x)$ 的 $m=k=4$ 阶有理逼近时，前五次迭代的偏差情况。可以看出，第一次迭代后，结果几乎与 minimax 解一样好。

图 5.13.1。实线曲线显示了对某个任意测试问题，ratlsq 例程连续五次迭代的偏差 $r(x)$。该算法并未收敛到精确的 minimax 解（以虚线显示）。但在一次迭代后，偏差已仅为最低有效位精度的一个小数部分。

迭代的收敛顺序并不如图中所示那样单调。事实上，第二次迭代的结果最佳（具有最小的最大偏差）。因此，ratlsq 例程返回的是所有迭代中最佳的结果，而不一定是最后一次迭代的结果；进行超过五次迭代并无优势。

参考文献

Ralston, A. and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley) Chapter 13.[1]

5.14 通过路径积分求函数值

在计算机编程中，首选的技术未必是最高效、最优雅或执行最快的。相反，它可能是实现迅速、通用性强且易于验证的方法。

有时我们仅需对某个特殊函数进行少量（或几千次）求值，该函数可能是复变量的复值函数，并具有许多不同的参数、渐近区域，或两者兼有。使用常规技巧（级数、连分式、有理函数逼近、递推关系等）可能导致一个拼凑而成的程序，其中包含大量条件判断和分支跳转到不同公式。虽然这样的程序在执行时可能非常高效，但从零开始实现它往往不是最快捷的途径。

另一种具有相当普遍性的技术是直接积分函数的定义微分方程——对每个所需的函数值都从头进行一次积分——必要时沿复平面上的一条路径进行。虽然这乍看之下像是用金砖拍苍蝇，但事实证明，当你已经拥有这块金砖，而苍蝇恰好就在砖下熟睡时，你只需让它落下即可！

作为一个具体例子，我们考虑复超几何函数 ${}_2{F}_1(a,b,c;z)$，它被定义为所谓超几何级数的解析延拓：

$$\begin{array}{r}{}_2{F}_1(a,b,c;z)=1+\frac{ab}{c}\frac{z}{1!}+\frac{a(a+1)b(b+1)}{c(c+1)}\frac{z^2}{2!}+\cdots \\ +\frac{a(a+1)\dots (a+j-1)b(b+1)\dots (b+j-1)}{c(c+1)\dots (c+j-1)}\frac{z^j}{j!}+\cdots \end{array}$$

该级数仅在单位圆 $|z|<1$ 内收敛（参见 [1]），但人们对这一函数的兴趣通常并不局限于该区域。

超几何函数 ${}_2{F}_1$ 是超几何微分方程的一个解（实际上是原点处正则的解），该方程可写作：

$$z(1-z)F^{\prime \prime }=abF-[c-(a+b+1)z]F^{\prime }$$

此处的撇号表示对 $z$ 求导（即 $d/dz$）。可以看出，该方程在 $z=0,1$ 和 $\infty$ 处具有正则奇点。由于我们所求的解在 $z=0$ 处是正则的，因此一般而言，点 $1$ 和 $\infty$ 将成为分支点。若希望 ${}_2{F}_1$ 是单值函数，则必须在两点之间设置一条分支割线。通常将该割线置于从 $1$ 到 $\infty$ 的正实轴上，尽管在某些应用中我们可能希望保留改变这一选择的可能性。

我们的“金砖”由一组用于积分常微分方程组的例程组成，这些例程将在第 17 章中详细展开。目前，我们只需要一个高层的“黑箱”例程，它能够从某个（实）自变量值处的初始条件出发，积分到该自变量另一值处的终值条件，同时自动调整内部步长以维持指定的精度。该例程称为 Odeint，在一种特定调用方式下，它使用一种复杂的 Bulirsch-Stoer 方法来计算每一步。

假设我们在某点 $z_0$ 处已知 $F$ 及其导数 $F^{\prime }$ 的值，并希望求出复平面上另一点 $z_1$ 处的 $F$ 值。连接这两点的直线路径可参数化为：

$$z(s)=z_0+s(z_1-z_0)$$

其中 $s$ 为实参数。现在可将微分方程 (5.14.2) 改写为两个一阶方程组成的方程组：

$$\begin{array}{rl}& \frac{dF}{ds}=(z_1-z_0)F^{\prime }\\ & \frac{d{F}^{\prime }}{ds}=(z_1-z_0)\left(\frac{abF-[c-(a+b+1)z]F^{\prime }}{z(1-z)}\right)\end{array}$$

图 5.14.1. 复平面示意图，显示了超几何函数的奇点、分支割线，以及从圆 $|z|=1/2$（幂级数在此快速收敛）到平面上其他点的一些积分路径。

该方程组需从 $s=0$ 积分到 $s=1$。此处 $F$ 和 $F^{\prime }$ 应视为两个相互独立的复变量。撇号表示对 $z$ 求导这一事实可以忽略；它将作为方程组 (5.14.4) 中第一个方程的自然结果而出现。此外，方程 (5.14.4) 的实部和虚部分别构成四个实微分方程，自变量为 $s$。右侧的复数运算可视为对四个分量如何耦合的一种简写。正是这种观点被传递给 Odeint 例程，因为该例程对复函数或复自变量一无所知。

接下来只需决定从何处出发，以及在复平面上选择哪条路径，才能到达任意给定点 $z$。这正是需要考虑函数奇点及所选分支割线的地方。图 5.14.1 展示了我们采用的策略。对于 $|z|\le 1/2$，方程 (5.14.1) 中的级数通常会快速收敛，此时直接使用该级数是合理的。否则，我们沿直线路径从起始点 $(\pm 1/2,0)$ 或 $(0,\pm 1/2)$ 之一进行积分。前两种选择分别适用于 $0<\text{Re}(z)<1$ 和 $\text{Re}(z)<0$ 的情形。后两种选择用于 $\text{Re}(z)>1$ 且位于分支割线上方或下方的情形；在这些情况下，从偏离实轴的位置出发是为了避免路径过于靠近 $z=1$ 处的奇点（见图 5.14.1）。分支割线的位置由我们所采用的策略定义：该策略永远不会在 $\text{Re}(z)>1$ 的区域穿越实轴。

该算法的一个具体实现在 §6.13 中给出，即例程 hypgeo。

上述过程存在若干变体，且易于编程实现。如果连续调用时 $z$ 的取值彼此接近（且 $a$、$b$、$c$ 的值相同），则可在每次调用时保存状态向量 $(F,F^{\prime })$ 及对应的 $z$ 值，并将其作为下一次调用的初始值。此时增量积分可能仅需一两步即可完成。注意避免无意中穿越分支割线：函数值虽“正确”，但可能并非你想要的那个值。

另外，你或许希望通过一条“折线”路径（dog-leg path）积分到某点 $z$，即使该路径穿越了 $\text{Re}(z)>1$ 的实轴，以此移动分支割线。例如，在某些情况下，你可能希望先从 $(0,1/2)$ 积分到 $(3/2,1/2)$，再由此到达任意满足 $\text{Re}(z)>1$ 的点（无论 $\text{Im}(z)$ 为正或负）。（例如，若你通过迭代法求函数的根，则不希望对邻近点的积分绕分支点采取不同路径。否则，你的求根器将看到不连续的函数值，很可能无法正确收敛！）

无论如何，需注意：若在积分路径上经过函数值很大的区域，而最终目标点处函数值很小，则可能导致数值精度损失。（对于超几何函数，一个典型情形是当 $a$ 和 $b$ 均为较大正数，且 $c$ 与 $x\gtrsim 1$ 时。）在此类情况下，你需要寻找一条更优的折线路径。

通过在复平面上积分其微分方程来计算函数值的通用技术，也可应用于其他特殊函数。例如，复贝塞尔函数、艾里函数、库仑波函数和韦伯函数均为合流超几何函数的特例，其微分方程与上述形式类似（参见 [1] §13.6 中的特例表）。合流超几何函数在有限 $z$ 处无奇点：这使其易于积分。然而，它在无穷远处具有本性奇点，这意味着在某些路径和参数下，函数可能呈现高度振荡或指数衰减的行为：这使其难以积分。因此，需要根据具体情况作出判断（或进行实验）。

参考文献

Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions (Washington: National Bureau of Standards); reprinted 1968 (New York: Dover); online at http://www.nr.com/aands.[1]

本章解析PDF