数值分析：从数学理论到 C++ 工程实现

“数值分析是数学与计算机科学的握手之处。” —— 本专题聚焦插值、数值积分、非线性方程求解及其在 C++ 中的高性能工程实现，旨在培养严谨的误差意识与算法落地能力。

---

🪜 练习阶梯与评价标准

等级	难度目标	核心考察点	期望达成
● Level A	公式直接应用	插值余项、代数精度计算	掌握误差上界评估逻辑
● Level B	算法工程化	复化求积、Newton 迭代	完成鲁棒的 C++ 函数封装
● Level C	综合稳定性	样条插值、病态方程组	解决数值不稳定性与收敛问题

---

🎯 跨学科考点矩阵 (Knowledge Matrix)

知识模块	数学核心	计算机落地 (C++)	关联习题
插值论	Lagrange/Newton 插值, 龙格现象	动态多项式构建, 复杂度 $O(n^2)$	练习 1, 2
数值积分	Newton-Cotes, Romberg 外推	复化 Simpson, 递归精度控制	练习 3, 4
方程求根	二分法, Newton-Raphson 迭代	容差控制 (eps), 收敛阶监测	练习 5, 6
工程建模	最小二乘拟合, 样条曲线	矩阵库应用 (Eigen 思想), $O(n)$ 追赶法	练习 7, 8

---

📂 核心习题库

Level A：理论与误差分析

练习 1：Lagrange 插值余项估计

题目描述：利用 $n=2$ 阶 Lagrange 插值多项式 $L_2(x)$ 逼近 $f(x) = \sin(x)$，节点分布在 $[0, 1]$ 上。试求最大可能误差上界。

Check Solution (Formal Proof)

思维链 (Thought Chain)：

1. 建模：设节点为 $x_0, x_1, x_2 \in [0, 1]$。
2. 推导：
   • 余项公式：$R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)$。
   • 对于 $f(x) = \sin(x)$，$f^{(3)}(x) = -\cos(x)$，故 $\max |f^{(3)}(\xi)| = 1$。
   • 令 $\omega_3(x) = |x(x-x_1)(x-1)|$。若取等距节点 $x_1=0.5$，则 $\max \omega_3(x) = \frac{\sqrt{3}}{36} \approx 0.048$。
3. 计算： $|R_2(x)| \le \frac{1}{3!} \cdot 1 \cdot 0.048 \approx 0.008$。

结论：误差上界为 $0.008$。

---

Level B：C++ 算法实战

练习 2：Newton 插值法的工程实现

题目描述：实现 Newton 插值算法。给定 $n$ 个点 $(x_i, y_i)$，计算均差表并构造插值多项式。

Check Solution (C++ Implementation)

解题思路： Newton 插值相比 Lagrange 插值，在增加新节点时无需重新计算所有项，具有更好的增量性。

C++ 代码实现：

#include <iostream>
#include <vector>

/**
 * @brief Newton Interpolation Class
 * 实现均差表的构建与插值计算
 */
class NewtonInterpolator {
private:
    std::vector<double> x_coords;
    std::vector<double> coefficients; // 均差表对角线元素

public:
    void addPoint(double x, double y) {
        x_coords.push_back(x);
        int n = x_coords.size();
        std::vector<double> temp_y(n);
        
        // 此处为简化逻辑：实际工程中建议维护二维表或一维压缩表
        // 演示目的采用动态构建
        if (n == 1) {
            coefficients.push_back(y);
        } else {
            // 计算新的均差
            double current_f = y;
            std::vector<double> diffs = {y};
            // ... (省略复杂的 O(n^2) 均差表逻辑，直接展示核心)
        }
    }

    // 工业级实现：一次性构建
    NewtonInterpolator(const std::vector<double>& x, const std::vector<double>& y) 
        : x_coords(x) {
        int n = x.size();
        std::vector<std::vector<double>> f(n, std::vector<double>(n));
        for (int i = 0; i < n; i++) f[i][0] = y[i];
        
        for (int j = 1; j < n; j++) {
            for (int i = j; i < n; i++) {
                f[i][j] = (f[i][j-1] - f[i-1][j-1]) / (x[i] - x[i-j]);
            }
        }
        for (int i = 0; i < n; i++) coefficients.push_back(f[i][i]);
    }

    double evaluate(double x) {
        double res = coefficients[0];
        double term = 1.0;
        for (int i = 1; i < coefficients.size(); i++) {
            term *= (x - x_coords[i-1]);
            res += coefficients[i] * term;
        }
        return res;
    }
};

int main() {
    std::vector<double> x = {0, 1, 2};
    std::vector<double> y = {1, 2, 4}; // f(x) = 2^x 的离散点
    NewtonInterpolator ni(x, y);
    std::cout << "f(1.5) ≈ " << ni.evaluate(1.5) << std::endl;
    return 0;
}

练习 3：复化 Simpson 积分

题目描述：编写一个通用的数值积分函数，采用复化 Simpson 公式计算 $\int_a^b f(x)dx$。

Check Solution (C++ Implementation)

数学原理： $\int_a^b f(x)dx \approx \frac{h}{3} [f(a) + 4\sum f(x_{2i-1}) + 2\sum f(x_{2i}) + f(b)]$

C++ 代码实现：

#include <iostream>
#include <functional>
#include <cmath>

double simpson(std::function<double(double)> f, double a, double b, int n) {
    if (n % 2 != 0) n++; // Simpson 要求区间数为偶数
    double h = (b - a) / n;
    double s1 = 0, s2 = 0;
    
    for (int i = 1; i < n; i += 2) s1 += f(a + i * h);
    for (int i = 2; i < n; i += 2) s2 += f(a + i * h);
    
    return (h / 3.0) * (f(a) + 4.0 * s1 + 2.0 * s2 + f(b));
}

int main() {
    auto target = [](double x) { return std::sin(x); };
    double res = simpson(target, 0, M_PI, 100);
    std::cout << "Integral of sin(x) from 0 to PI: " << res << std::endl;
    return 0;
}

---

Level C：挑战与高级应用

练习 4：Newton-Raphson 迭代的鲁棒性

题目描述：实现 Newton 迭代法求解 $f(x)=0$。要求：

1. 具备最大迭代次数限制。
2. 具备容差控制。
3. 具备导数极小时的异常处理（防除零）。

Check Solution (C++ Robust Implementation)

C++ 代码实现：

#include <iostream>
#include <functional>
#include <stdexcept>
#include <cmath>

class NewtonSolver {
public:
    static double solve(std::function<double(double)> f, 
                        std::function<double(double)> df, 
                        double x0, 
                        double tol = 1e-7, 
                        int max_iter = 100) {
        double x = x0;
        for (int i = 0; i < max_iter; i++) {
            double fx = f(x);
            double dfx = df(x);
            
            if (std::abs(dfx) < 1e-12) {
                throw std::runtime_error("Derivative too small, possible singularity.");
            }
            
            double dx = fx / dfx;
            x -= dx;
            
            if (std::abs(dx) < tol) return x;
        }
        throw std::runtime_error("Newton method failed to converge.");
    }
};

int main() {
    auto f = [](double x) { return x * x - 2.0; };
    auto df = [](double x) { return 2.0 * x; };
    try {
        double root = NewtonSolver::solve(f, df, 1.0);
        std::cout << "sqrt(2) ≈ " << root << std::endl;
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << std::endl;
    }
    return 0;
}

---

🏆 训练建议

1. 精度敏感性：在 C++ 中务必使用 double 或 long double，并了解 std::numeric_limits<double>::epsilon() 的含义。
2. 稳定性高于速度：在数值算法中，防止分母过小引起的溢出往往比追求 $O(n)$ 复杂度更重要。
3. 跨学科思考：尝试将数值积分应用于物理模拟（如求质心），将插值应用于图形学曲线平滑。