数值分析——最佳平方逼近及其计算（C++版本）

摘要本文基于李庆扬《数值分析》第五版系统梳理了最小二乘逼近中法方程的推导过程从平方误差积分最小化出发通过偏导数为零、内积符号引入和矩阵化最终得到法方程。随后给出法方程解存在唯一性的严谨证明并深入推导了误差公式与希尔伯特矩阵的由来。文章还提供了完整的 C 代码直观演示希尔伯特矩阵在 n10 时的病态效应微小扰动导致系数剧烈震荡、残差极小但解不可靠。最后指出工程上的终极解决方案——采用勒让德等正交多项式使格拉姆矩阵对角化、条件数趋近于 1从而彻底规避病态问题。1.李庆扬《数值分析》第五版资料2.法方程的推导第一步把问题写成“数学式子”目标找一个多项式让它逼近原函数 f(x)。误差标准平方误差的加权积分最小注意这里的I是关于未知系数 a0,a1,…,an 的多元函数。我们的任务就是找一组 a让这个积分最小。第二步求极值就是“偏导数为0”第三步把求和号拆开引入“内积”符号第四步从“方程组”变成“矩阵方程”3.证明1. 目标我们要证明什么2. 关键的代数展开3. 致命一击为什么第一项积分为 04. 最终的结论4.证明δ(x希尔伯特矩阵1.误差公式的严谨推导2.希尔伯特矩阵的推导5.C代码希尔伯特矩阵的病态效应1.代码#include iostream #include cmath #include iomanip #include functional #include vector #include stdexcept using namespace std; // 1. 通用数值积分工具 double simpsonIntegrate(functiondouble(double) f, double a, double b, int n 100000) { if (n % 2 ! 0) n; double h (b - a) / n; double sum f(a) f(b); for (int i 1; i n; i) { double x a i * h; sum (i % 2 0) ? 2 * f(x) : 4 * f(x); } return sum * h / 3.0; } double innerProduct(functiondouble(double) f, functiondouble(double) g, double a, double b) { return simpsonIntegrate([](double x) { return f(x) * g(x); }, a, b); } // 2. 通用线性方程组求解器高斯消元 列主元 // 替换之前的 solve2x2用于求解任意 n*n 的矩阵 vectordouble solveLinearSystem(vectorvectordouble A, vectordouble b) { int n A.size(); for (int i 0; i n; i) { // 列主元选取找出当前列绝对值最大的元素避免除数过小导致误差爆炸 int pivot i; for (int j i 1; j n; j) { if (abs(A[j][i]) abs(A[pivot][i])) pivot j; } swap(A[i], A[pivot]); swap(b[i], b[pivot]); if (abs(A[i][i]) lt; 1e-15) throw runtime_error(矩阵奇异); // 消元过程 for (int j i 1; j lt; n; j) { double factor A[j][i] / A[i][i]; for (int k i; k lt; n; k) { A[j][k] - factor * A[i][k]; } b[j] - factor * b[i]; } } // 回代过程 vectorlt;doublegt; x(n); for (int i n - 1; i gt; 0; --i) { double sum 0.0; for (int j i 1; j lt; n; j) { sum A[i][j] * x[j]; } x[i] (b[i] - sum) / A[i][i]; } return x; } // 3. 主函数 int main() { double a 0.0, b 1.0; auto f [](double x) { return sqrt(1.0 x * x); }; cout lt;lt; lt;lt; endl; cout lt;lt; 演示希尔伯特矩阵的病态效应 (n2 vs n10) lt;lt; endl; cout lt;lt; lt;lt; endl lt;lt; endl; // 实验一2阶逼近正常情况 { int n 2; cout lt;lt; 【实验一】n 2 阶逼近 (良态计算正常) lt;lt; endl; vectorlt;vectorlt;doublegt;gt; H(n, vectorlt;doublegt;(n, 0.0)); vectorlt;doublegt; d(n, 0.0); // 【修复处 1】给 lambda 加上按值捕获 [i] 和 [j] for (int i 0; i lt; n; i) { for (int j 0; j lt; n; j) { H[i][j] innerProduct([i](double x) { return pow(x, i); }, [j](double x) { return pow(x, j); }, a, b); } d[i] innerProduct(f, [i](double x) { return pow(x, i); }, a, b); } vectorlt;doublegt; coef solveLinearSystem(H, d); cout lt;lt; 解得的系数 a0, a1: ; for (double c : coef) cout lt;lt; fixed lt;lt; setprecision(4) lt;lt; c lt;lt; ; cout lt;lt; \n与教材例6吻合\n\n; } // 实验二10阶逼近病态灾难 { int n 10; // 10次多项式有11个系数 cout lt;lt; 【实验二】n 10 阶逼近 (希尔伯特矩阵病态测试) lt;lt; endl; vectorlt;vectorlt;doublegt;gt; H(n 1, vectorlt;doublegt;(n 1, 0.0)); vectorlt;doublegt; d(n 1, 0.0); // 【修复处 2】同样的用按值捕获 [i] 和 [j] for (int i 0; i lt; n; i) { for (int j 0; j lt; n; j) { H[i][j] innerProduct([i](double x) { return pow(x, i); }, [j](double x) { return pow(x, j); }, a, b); } d[i] innerProduct(f, [i](double x) { return pow(x, i); }, a, b); } // 求解 vectorlt;doublegt; coef solveLinearSystem(H, d); cout lt;lt; 解得的系数 a0 ~ a10 lt;lt; endl; cout lt;lt; scientific lt;lt; setprecision(4); for (int i 0; i lt; n; i) { cout lt;lt; a lt;lt; setw(2) lt;lt; i lt;lt; lt;lt; coef[i] lt;lt; endl; } // 计算真实的逼近误差 double error_sq innerProduct(f, f, a, b); for (int i 0; i lt; n; i) { error_sq - coef[i] * d[i]; } cout lt;lt; \n真实逼近误差 ||delta||^2 lt;lt; error_sq lt;lt; endl; cout lt;lt; 注意观察虽然误差看起来还行但系数已经彻底失控正负剧烈震荡\n\n; // 实验三微小扰动测试照妖镜 cout lt;lt; 【实验三】微小扰动测试在右端项 d 加入 1e-8 的误差 lt;lt; endl; vectorlt;doublegt; d_perturbed d; d_perturbed[0] 1e-8; // 仅在第0项加一个极小的误差 vectorlt;doublegt; coef_perturbed solveLinearSystem(H, d_perturbed); cout lt;lt; 扰动后的系数 a0 ~ a10观察变化幅度 lt;lt; endl; for (int i 0; i lt; n; i) { double change abs(coef_perturbed[i] - coef[i]); cout lt;lt; a lt;lt; setw(2) lt;lt; i lt;lt; lt;lt; coef_perturbed[i] lt;lt; (变化量: lt;lt; change lt;lt; ) lt;lt; endl; } cout lt;lt; \n【结论】 lt;lt; endl; cout lt;lt; 右端项仅仅改变了 1e-8但系数变化了极其巨大的倍数。 lt;lt; endl; cout lt;lt; 这就是希尔伯特矩阵在 n10 时的恐怖病态效应 lt;lt; endl; } return 0; }2.运行结果3.病态矩阵解释1.为什么说这绝对是病态铁证1微小的输入扰动导致输出的巨大震荡看你的【实验三】。你仅仅在右端项 d0​ 上加了10E−8的误差相当于小数点后8位的一点点抖动。再看系数变化a5​ 变化了 0.229放大了一千万倍a6​ 变化了 0.465a8​ 变化了 0.475输入 10e−8 的扰动输出 0.46 的变化放大倍率达到了惊人的 10e7 级别。这说明这个 11×11 的矩阵极其脆弱微小的舍入误差在求解过程中被无限放大了。铁证2系数正负剧烈震荡看你的【实验二】结果a01.0a1≈0a2≈0.5a3≈0a4−0.125a5≈0a60.069a7−0.012a8​−0.04a90.03a10−0.007。你会发现系数在正负之间疯狂横跳。虽然理论误差极小但高次项如 a8,a9,a10​根本没有衰减到接近零反而出现了不该有的非零值。这就是病态导致高阶系数“失控”的表现。2.为什么误差极小-4.2e-15但系数却彻底失控这是一个极其经典、容易让人困惑的问题。你可能会想“误差这么小难道不是好事吗”这里的“误差”是残差残差 Ha−d它衡量的是“方程是否被满足”而不是“解是否可靠”。在病态矩阵中存在一条“峡谷”在峡谷底部有很多个不同的 a它们都能让方程 Ha≈d 的残差几乎为零所以你看到了 10E−15 级别的误差。但当你施加一点点扰动10E−8解就会在峡谷壁上发生巨大的位移系数变化 0.46这就是解的不稳定性。在工程上“解不稳定”比“残差大”更可怕。这意味着你算出来的系数是不可信的哪怕方程看起来被满足得很好。3.为什么希尔伯特矩阵会这么恐怖4.工程上的终极解决方案正交多项式现在明白为什么教材 3.3.2 节要立刻引入正交多项式了吧这是绝境中的唯一生路。如果我们将基函数换成勒让德多项式Legendre Polynomials它们在区间 [−1,1] 上两两正交内积为0。格拉姆矩阵 G 会直接变成一个对角矩阵对角线上只有 2/(2k1)其余全是0。条件数 κ1完美良态。解方程退化成简单的除法。无论 n 取多大10、100、1000计算都会稳如泰山绝不失控。