数值计算方法：收敛性、稳定性与算法选择简述

数值计算方法：收敛性、稳定性与算法选择简述

数值计算研究的是如何用有限精度的机器运算，近似解决通常难以直接得到解析解的问题。一个可靠的数值算法不仅要给出结果，还必须回答三个问题：

1. 结果是否趋向目标问题的真实解；
2. 输入扰动和舍入误差是否会被放大；
3. 结果是否能够通过残差、误差估计或独立方法验证。

本文从统一的误差概念出发，依次介绍非线性方程求根、线性方程组、插值、数值积分和最小二乘拟合，并给出适合工程实现的算法框架。

1. 数值计算的三个基本概念

1.1 误差、残差与条件数

设精确结果为 \(x\)，计算结果为 \(\hat{x}\)。绝对误差为：

\[e_{abs}=|\hat{x}-x| \]

相对误差为：

\[e_{rel}=\frac{|\hat{x}-x|}{|x|},\qquad x\neq0 \]

当精确解未知时，通常使用残差。例如方程 \(f(x)=0\) 的残差是：

\[r_f=|f(\hat{x})| \]

线性系统 \(Ax=b\) 的残差是：

\[r_A=A\hat{x}-b \]

残差小并不自动等价于解误差小。问题本身的条件数决定了输入扰动对输出的敏感程度。对线性系统，有近似关系：

\[\frac{\|\delta x\|}{\|x\|} \lesssim \kappa(A)\frac{\|\delta b\|}{\|b\|} \]

条件数很大时，即使算法稳定，输入数据的微小误差也可能导致解发生明显变化。

1.2 收敛性与稳定性

• 收敛性：随着迭代次数增加，算法是否趋向目标解。
• 稳定性：输入误差、舍入误差和中间计算误差是否被过度放大。
• 一致性：离散算法是否在步长趋于零时逼近原问题。

对于迭代算法，收敛速度通常分为：

• 线性收敛：误差大致按固定比例缩小；
• 超线性收敛：后期收敛速度超过线性；
• 二次收敛：误差近似满足 e_{k+1}=Ce_k^2。

1.3 浮点数与停止条件

浮点数不能使用“两个结果完全相等”作为收敛判断。统一的停止条件应采用绝对容差和相对容差：

\[\|x_{k+1}-x_k\| \leq \varepsilon_{abs}+\varepsilon_{rel}\|x_{k+1}\| \]

同时设置最大迭代次数，并在出现 NaN、无穷大、残差持续增长或中间量超出范围时停止计算。

2. 非线性方程求根

2.1 问题形式

非线性方程求根要求找到满足：

\[f(x)=0 \]

的 x。例如：

\[f(x)=x^3-\sin x-12x+1 \]

其导数为：

\[f'(x)=3x^2-\cos x-12 \]

2.2 二分法：可靠的区间方法

如果连续函数在 [a,b] 上满足：

\[f(a)f(b)<0 \]

根据介值定理，区间内至少存在一个根。每轮取中点：

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

根据 f(a)f(m) 的符号决定保留左半区间还是右半区间。

第 k 轮后区间长度为：

\[\frac{b-a}{2^k} \]

因此二分法的误差具有明确上界。它的缺点是只具有线性收敛速度，但优点是对初值不敏感、过程稳定、容易验证。

2.3 固定点迭代

将方程改写成：

\[x=\varphi(x) \]

然后使用：

\[x_{k+1}=\varphi(x_k) \]

对于示例函数，可以构造：

\[\varphi(x)=\frac{x^3-\sin x+1}{12} \]

局部收敛通常要求根附近满足：

\[|\varphi'(x)|<1 \]

更强的充分条件是在闭区间 I 上满足：

\[\varphi(I)\subseteq I \]

且：

\[\max_{x\in I}|\varphi'(x)|\leq q<1 \]

此时 φ 是压缩映射，迭代从区间内任意初值出发都收敛到唯一不动点。

固定点迭代的关键不在于“能否改写”，而在于如何选择具有收敛性的 φ。

2.4 Newton 法

在 x_k 处用切线近似函数：

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

得到：

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

当根是单根、初值足够接近根且函数满足一定光滑性条件时，Newton 法具有局部二次收敛：

\[|e_{k+1}|\approx C|e_k|^2 \]

但 Newton 法可能遇到：

• 导数接近零；
• 迭代点跳出有效区间；
• 初值落入另一个根的吸引域；
• 函数存在多个根或高阶根。

2.5 混合求根策略

实际应用中可以采用两阶段策略：

1. 使用二分法或其他区间方法保证迭代点进入目标根附近；
2. 使用 Newton 法快速提高精度；
3. 如果 Newton 步长跳出区间，则退回二分步；
4. 同时检查步长、残差、最大次数和数值合法性。

方法	导数	区间	局部速度	稳健性
二分法	否	是	线性	高
固定点迭代	否	否	取决于 φ	中
Newton 法	是	否	二次	依赖初值
混合方法	是或否	是	兼顾速度和稳定性	高

3. 线性方程组

3.1 增广矩阵与消元

线性系统写作：

\[Ax=b \]

通过增广矩阵：

\[\left[\begin{array}{cccc|c} a_{11}&a_{12}&\cdots&a_{1n}&b_1\\ a_{21}&a_{22}&\cdots&a_{2n}&b_2\\ \vdots&\vdots&&\vdots&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{nn}&b_n \end{array}\right] \]

可以同时对系数和右端项进行初等行变换。

3.2 Gauss 消元

第 i 列对第 j 行的消元因子为：

\[q_{ji}=\frac{a_{ji}}{a_{ii}},\qquad j>i \]

行变换为：

\[R_j\leftarrow R_j-q_{ji}R_i \]

前向消元得到上三角矩阵后，通过回代计算：

\[x_i=\frac{b_i-\sum_{j=i+1}^{n}a_{ij}x_j}{a_{ii}} \]

稠密矩阵的时间复杂度为 O(n^3)，存储复杂度为 O(n^2)。

3.3 列主元与 LU 分解

为了避免小主元放大误差，在第 i 列选择：

\[p=\arg\max_{j\geq i}|a_{ji}| \]

并交换第 p 行和第 i 行。带列主元的 LU 分解可以表示为：

\[PA=LU \]

当矩阵 A 不变而右端项有多个时，只需分解一次，然后对每个 b 执行：

\[Ly=Pb,\qquad Ux=y \]

这比每次重新消元更高效。

实现约束如下：

• 矩阵元素和交换临时变量统一使用浮点类型；
• 每次除法前检查主元是否接近零；
• 奇异矩阵直接返回失败状态；
• 求解完成后计算 r=A\hat{x}-b；
• 多右端项复用同一份分解结果。

3.4 迭代法的适用场景

大型稀疏系统不适合直接保存完整稠密矩阵。此时使用 Jacobi、Gauss-Seidel、SOR 或更高级的 Krylov 子空间方法，并结合稀疏存储和预条件技术。

4. 迭代法求解线性系统

4.1 统一迭代形式

线性迭代法可统一写作：

\[x^{(k+1)}=Bx^{(k)}+c \]

误差满足：

\[e^{(k+1)}=Be^{(k)} \]

收敛的核心条件为：

\[\rho(B)<1 \]

其中 ρ(B) 是迭代矩阵的谱半径。严格对角占优和对称正定等性质，可以作为部分方法的充分条件。

4.2 Jacobi 迭代

将矩阵分解为：

\[A=D+L+U \]

Jacobi 迭代为：

\[x^{(k+1)}=D^{-1}[b-(L+U)x^{(k)}] \]

逐分量形式为：

\[x_i^{(k+1)}= \frac{b_i-\sum_{j\neq i}a_{ij}x_j^{(k)}}{a_{ii}} \]

每个新分量都使用上一轮的旧向量，因此适合并行计算。

4.3 Gauss-Seidel 迭代

Gauss-Seidel 立即使用已经更新的分量：

\[x_i^{(k+1)}= \frac{b_i- \sum_{j<i}a_{ij}x_j^{(k+1)}- \sum_{j>i}a_{ij}x_j^{(k)}}{a_{ii}} \]

它通常比 Jacobi 收敛更快，但计算顺序存在依赖。

4.4 SOR 迭代

SOR 引入松弛因子 ω：

\[x_i^{(k+1)}=(1-\omega)x_i^{(k)}+ \frac{\omega}{a_{ii}} \left(b_i- \sum_{j<i}a_{ij}x_j^{(k+1)}- \sum_{j>i}a_{ij}x_j^{(k)}\right) \]

其中：

• ω=1 为 Gauss-Seidel；
• 0<ω<1 为欠松弛；
• 1<ω<2 为超松弛。

迭代器统一使用以下终止逻辑：

1. 计算更新量 ||x_{k+1}-x_k||；
2. 计算残差 ||Ax_{k+1}-b||；
3. 两者都低于阈值时返回成功；
4. 残差持续增长、出现非法浮点值或达到最大次数时返回未收敛。

5. 插值方法

设节点为：

\[(x_0,y_0),(x_1,y_1),\ldots,(x_{n-1},y_{n-1}) \]

插值要求构造函数 p(x)，使其满足 p(x_i)=y_i。

5.1 Lagrange 插值

\[P_{n-1}(x)=\sum_{i=0}^{n-1}y_iL_i(x) \]

\[L_i(x)=\prod_{j=0,j\neq i}^{n-1}\frac{x-x_j}{x_i-x_j} \]

误差为：

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

直接计算一个查询点的复杂度为 O(n^2)。当查询点较多时，使用重心形式：

\[w_i=\frac{1}{\prod_{j\neq i}(x_i-x_j)} \]

\[P(x)= \frac{\sum_i\frac{w_iy_i}{x-x_i}} {\sum_i\frac{w_i}{x-x_i}} \]

预处理权重需要 O(n^2)，每个查询点的计算降为 O(n)。节点必须非空、互异；查询点若恰好等于某个节点，则直接返回该节点值。

5.2 三次样条

三次样条在每个区间 [x_i,x_{i+1}] 上使用：

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

相邻区间满足函数值、一阶导数和二阶导数连续。自然边界条件为：

\[S''(x_0)=0,\qquad S''(x_{n-1})=0 \]

以节点二阶导数 M_i 为未知量时，内部节点形成三对角系统：

\[\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i \]

其中：

\[h_i=x_{i+1}-x_i \]

\[\mu_i=\frac{h_{i-1}}{h_{i-1}+h_i},\qquad \lambda_i=\frac{h_i}{h_{i-1}+h_i} \]

三对角系统使用追赶法求解，建表和求解复杂度均为 O(n)。节点必须严格递增；查询点必须位于 [x_0,x_{n-1}] 内；最右端点单独返回 y_{n-1}；区间定位使用二分查找。

与高次 Lagrange 多项式相比，三次样条通过局部多项式降低了端点振荡，代价是需要明确边界条件。

6. 数值积分

6.1 复化梯形公式

对区间 [a,b] 等分为 n 段：

\[h=\frac{b-a}{n} \]

复化梯形公式为：

\[T_n=\frac{h}{2}\left[f(a)+f(b)+2\sum_{k=1}^{n-1}f(a+kh)\right] \]

当 f 二阶连续可导时：

\[I-T_n=O(h^2) \]

6.2 Romberg 外推

若：

\[T(h)=I+c_2h^2+c_4h^4+c_6h^6+\cdots \]

则可以递推消除误差项：

\[R_{k,0}=T_{2^k} \]

\[R_{k,j}=R_{k,j-1}+ \frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1} \]

前几级分别产生分母 3、15 和 63。

Romberg 实现使用滚动数组保存递推表，相邻层结果满足：

\[|R_{k,k}-R_{k-1,k-1}|<\varepsilon \]

时停止。分段数使用正整数，最大外推层数固定，超过上限则返回未收敛状态。

对于：

\[f(x)=\frac{\sin x}{x} \]

在 x=0 处使用极限值 f(0)=1，避免直接计算 0/0。

7. 最小二乘多项式拟合

给定数据 (x_i,y_i)，用 m 次多项式：

\[p_m(x)=c_0+c_1x+\cdots+c_mx^m \]

拟合系数通过最小化：

\[\min_c\|Vc-y\|_2^2 \]

得到，其中 V 是 Vandermonde 矩阵。

7.1 QR 分解

拟合流程采用 Householder QR：

\[V=QR \]

将最小二乘问题转化为：

\[Rc=Q^Ty \]

该过程不显式构造 V^TV，避免条件数平方放大。计算完成后输出：

\[\mathrm{SSE}=\|Vc-y\|_2^2 \]

7.2 秩亏与 SVD

当设计矩阵接近秩亏时，QR 分解会暴露出很小的对角元素。此时切换到 SVD：

\[V=U\Sigma V^T \]

根据奇异值阈值确定有效秩，并截断过小奇异值对应的方向，从而抑制噪声放大。

7.3 模型评价

三次多项式的参数更多，在同一批数据上的 SSE 通常不会大于二次多项式，但模型选择不只依赖 SSE。最终评价同时使用：

• 验证集误差；
• 残差是否存在系统性趋势；
• 区间边界处是否出现异常振荡；
• 模型参数数量与数据量是否匹配。

当观测数据具有不同方差时，使用加权最小二乘：

\[\min_c\sum_iw_i[p_m(x_i)-y_i]^2 \]

8. 统一的验证框架

不同算法使用不同的误差指标，但验证流程可以统一：

1. 验证输入合法性：维度、节点、区间、主元和容差；
2. 执行有最大次数或最大层数限制的算法；
3. 记录每次迭代的近似值、步长和残差；
4. 根据停止条件判定成功、失败或未收敛；
5. 对最终结果执行独立残差检查；
6. 在已知参考值时额外计算绝对误差和相对误差。

各类问题的核心指标如下：

问题	主要指标
方程求根	`
线性方程组	`
迭代线性系统	残差曲线、谱半径、迭代次数
插值	最大绝对误差、节点误差、边界振荡
数值积分	外推差值、参考值误差、层数
曲线拟合	SSE、验证集误差、有效秩

9. 一组可运行的参考实现

下面给出一组相互独立的 C++17 参考实现。它们不是某个特定项目的补丁，而是按照前文的数值约束重新组织的最小实现：

• 迭代过程不使用递归；
• 所有迭代都有最大次数；
• 结果携带收敛状态和诊断信息；
• 线性系统使用列主元；
• 插值节点经过合法性检查；
• 积分过程使用最大外推层数限制。

9.1 统一的迭代结果

将算法结果和诊断信息放在同一个结构中，可以避免通过特殊数值表示失败状态：

#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

struct IterationResult {
    double value = std::numeric_limits<double>::quiet_NaN();
    double residual = std::numeric_limits<double>::infinity();
    int iterations = 0;
    bool converged = false;
    std::string status;
};

9.2 二分法与 Newton 混合求根

下面的实现先验证变号区间，再优先执行 Newton 步。如果 Newton 步越出当前区间，或者导数过小，就回退到二分步。这样既保留了二分法的区间约束，也利用了 Newton 法的局部二次收敛。

IterationResult solve_root(
    const std::function<double(double)>& f,
    const std::function<double(double)>& df,
    double left,
    double right,
    double x0,
    double abs_tol = 1e-12,
    double residual_tol = 1e-12,
    int max_iter = 100)
{
    const double nan = std::numeric_limits<double>::quiet_NaN();
    if (!(left < right) || max_iter <= 0 || abs_tol <= 0 || residual_tol <= 0)
        return { nan, nan, 0, false, "invalid_input" };

    double fa = f(left);
    double fb = f(right);
    if (!std::isfinite(fa) || !std::isfinite(fb) || fa * fb > 0.0)
        return { nan, nan, 0, false, "invalid_bracket" };

    double x = std::clamp(x0, left, right);
    double fx = f(x);

    for (int k = 1; k <= max_iter; ++k)
    {
        double next = std::numeric_limits<double>::quiet_NaN();
        double dfx = df(x);

        if (std::isfinite(dfx) && std::abs(dfx) > 1e-14)
            next = x - fx / dfx;

        if (!std::isfinite(next) || next <= left || next >= right)
            next = 0.5 * (left + right);

        double fnext = f(next);
        if (!std::isfinite(fnext))
            return { next, fnext, k, false, "non_finite" };

        if (std::abs(fnext) < residual_tol &&
            std::abs(next - x) < abs_tol)
            return { next, std::abs(fnext), k, true, "converged" };

        if (fa * fnext <= 0.0)
        {
            right = next;
            fb = fnext;
        }
        else
        {
            left = next;
            fa = fnext;
        }

        x = next;
        fx = fnext;
    }

    return { x, std::abs(fx), max_iter, false, "max_iterations" };
}

对于：

\[f(x)=x^3-\sin x-12x+1 \]

调用时传入 f 和 f'，同时提供一个已经验证过的变号区间。函数返回值中，status 区分输入错误、非有限数值、达到最大次数和正常收敛。

9.3 带列主元的 Gauss-LU 求解

下面的实现原地保存 L 和 U，并通过行交换记录置换。它适合中小规模稠密矩阵；多组右端项时，可以复用分解后的矩阵。

using Matrix = std::vector<std::vector<double>>;

std::vector<double> solve_lu(Matrix a, std::vector<double> b)
{
    const int n = static_cast<int>(a.size());
    if (n == 0 || static_cast<int>(b.size()) != n)
        throw std::invalid_argument("invalid linear system");

    for (const auto& row : a)
        if (static_cast<int>(row.size()) != n)
            throw std::invalid_argument("matrix must be square");

    for (int k = 0; k < n; ++k)
    {
        int pivot = k;
        for (int i = k + 1; i < n; ++i)
            if (std::abs(a[i][k]) > std::abs(a[pivot][k]))
                pivot = i;

        if (std::abs(a[pivot][k]) < 1e-14)
            throw std::runtime_error("singular or ill-conditioned matrix");

        std::swap(a[k], a[pivot]);
        std::swap(b[k], b[pivot]);

        for (int i = k + 1; i < n; ++i)
        {
            a[i][k] /= a[k][k];
            for (int j = k + 1; j < n; ++j)
                a[i][j] -= a[i][k] * a[k][j];
            b[i] -= a[i][k] * b[k];
        }
    }

    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i)
    {
        double value = b[i];
        for (int j = i + 1; j < n; ++j)
            value -= a[i][j] * x[j];
        x[i] = value / a[i][i];
    }
    return x;
}

该实现的消元阶段复杂度为 O(n^3)，回代为 O(n^2)。实际应用中还应单独计算 ||Ax-b||，因为求解器返回成功并不等于问题本身条件良好。

9.4 重心 Lagrange 插值

重心形式将节点权重预计算，使多次查询更加高效：

class BarycentricInterpolator {
public:
    BarycentricInterpolator(std::vector<double> x,
                            std::vector<double> y)
        : x_(std::move(x)), y_(std::move(y)), w_(x_.size())
    {
        if (x_.empty() || x_.size() != y_.size())
            throw std::invalid_argument("invalid interpolation data");

        for (size_t i = 0; i < x_.size(); ++i)
        {
            w_[i] = 1.0;
            for (size_t j = 0; j < x_.size(); ++j)
            {
                if (i == j)
                    continue;
                double gap = x_[i] - x_[j];
                if (gap == 0.0)
                    throw std::invalid_argument("duplicate interpolation node");
                w_[i] /= gap;
            }
        }
    }

    double operator()(double query) const
    {
        double numerator = 0.0;
        double denominator = 0.0;

        for (size_t i = 0; i < x_.size(); ++i)
        {
            if (query == x_[i])
                return y_[i];

            double term = w_[i] / (query - x_[i]);
            numerator += term * y_[i];
            denominator += term;
        }
        return numerator / denominator;
    }

private:
    std::vector<double> x_, y_, w_;
};

预处理复杂度为 O(n^2)，单次查询复杂度为 O(n)。节点数量很大或节点高度聚集时，还应考虑节点缩放、Chebyshev 节点以及更稳定的基函数表示。

9.5 Romberg 积分

Romberg 表可以使用一维数组滚动更新：

double romberg(
    const std::function<double(double)>& f,
    double a,
    double b,
    double tol = 1e-10,
    int max_level = 12)
{
    if (!(a < b) || tol <= 0 || max_level < 1)
        throw std::invalid_argument("invalid integration parameters");

    std::vector<double> table(max_level + 1);
    double h = b - a;
    table[0] = 0.5 * h * (f(a) + f(b));

    for (int level = 1; level <= max_level; ++level)
    {
        const int new_points = 1 << (level - 1);
        const double step = h / (1 << level);
        double midpoint_sum = 0.0;

        for (int i = 1; i <= new_points; ++i)
            midpoint_sum += f(a + (2 * i - 1) * step);

        const std::vector<double> previous = table;
        table[0] = 0.5 * table[0] + step * midpoint_sum;

        for (int order = 1; order <= level; ++order)
        {
            const double factor = std::pow(4.0, order);
            table[order] = table[order - 1] +
                (table[order - 1] - previous[order - 1]) / (factor - 1.0);
        }

        if (std::abs(table[level] - table[level - 1]) < tol)
            return table[level];

        h *= 0.5;
    }

    throw std::runtime_error("Romberg integration did not converge");
}

每一层新增采样点，避免重复计算已有的梯形节点。滚动数组使额外空间为 O(level)；函数求值次数仍随外推层数指数增长，因此 max_level 必须受到限制。

10. 算法选择

问题特征	算法选择	关键条件
只有变号区间	二分法	函数连续且端点异号
有可靠初值和导数	Newton 法	导数不能过小
稠密线性系统	列主元 LU	主元和条件数可接受
多个右端项	LU 分解复用	系数矩阵不变
大型稀疏系统	迭代法	收敛条件和预条件有效
少量低阶节点	Lagrange	节点互异、次数较低
需要平滑局部曲线	三次样条	节点有序、边界条件明确
光滑函数积分	Romberg	梯形误差展开有效
小规模拟合	Householder QR	设计矩阵满秩或近似满秩
秩亏拟合	截断 SVD	奇异值阈值明确

11. 工程化实现

数值程序遵循以下结构：

• 输入校验层：检查维度、范围、节点、容差和数据类型；
• 算法层：只执行数学计算，不直接打印；
• 诊断层：统一返回近似值、残差、迭代次数、停止原因和收敛状态；
• 验证层：独立计算残差、参考误差和稳定性指标；
• 展示层：根据诊断结果生成表格、曲线或日志。

每个迭代算法都具有固定的停止状态：

• converged：满足步长和残差条件；
• max_iterations：达到最大次数；
• invalid_input：输入不满足算法要求；
• singular_or_ill_conditioned：主元或奇异值过小；
• non_finite：出现 NaN 或无穷大。

测试覆盖以下边界：

• 初值已经是精确解；
• 导数或主元接近零；
• 奇异矩阵和近奇异矩阵；
• 插值点位于左右端点和区间外；
• 节点重复或未排序；
• 误差容差过小；
• 迭代不收敛；
• 设计矩阵秩亏。

12. 总结

数值计算算法的可靠性来自完整的闭环：

\[\text{数学模型} \rightarrow \text{收敛分析} \rightarrow \text{稳定实现} \rightarrow \text{误差验证} \]

求根问题关注不动点、导数和初值；线性系统关注主元、条件数和残差；插值问题关注节点分布和边界条件；积分问题关注误差展开和外推层数；拟合问题关注矩阵秩、模型复杂度和泛化误差。

本文采用的实现体系是：二分法与 Newton 的混合求根、列主元 LU、带残差检测的迭代线性求解、重心 Lagrange、自然三次样条、滚动数组 Romberg、Householder QR 以及秩亏时的截断 SVD。它们共同构成了一套从理论、实现到验证都闭合的数值计算方法。