简介：广东工业大学数值分析课程实验四的Python实现资源，围绕超定方程组最小二乘求解与最小二乘拟合两大核心任务，适合正在学习数值分析、科学计算或需要完成同类上机实验的本科生参考。实验报告完整记录了题目三两个超定方程组的最小二乘解推导过程，以及题目六基于线性组合模型的拟合函数求解，并拓展至石油产量预测应用场景。报告包含算法原理说明、正规方程组建立步骤、具体数值结果与图表分析，同时附有基于NumPy、scikit-learn与Matplotlib的Python源代码，方便读者对照复现。资源为单个docx格式文档，压缩包共1个文件，大小112KB，内容结构清晰，从问题描述、算法介绍到结果讨论层层递进。已有212人在CSDN学习或下载过这份实验报告，对理解最小二乘方法在实际数据拟合与预测中的应用具有直接参考价值。

1. 数值分析实验四：一份把最小二乘从公式推到代码仓库的 Python 报告

如果你正在抄数值分析作业，或者正被超定方程组、最小二乘拟合、多项式预测三座大山压着，这份某高校数值分析课程的实验四报告值得你下载下来逐行跑一遍。它不是一个教学 Demo，而是一份“拿到题目、写出算法、用 numpy 和 sklearn 落地、再对结果做残差分析”的完整闭环，覆盖了三类高频考点：超定方程组的最小二乘解、可线性化的非线性拟合、以及多模型对比预测。整份报告只有三个 Python 脚本，但把最小二乘的正规方程解法、线性回归模型、polyfit 与 leastsq 两种拟合路径全串起来了。适合正在修数值分析、需要交上机报告，或者想看看别人怎么组织“算法 + 代码 + 结果分析”的开发者参考。

2. 超定方程组求最小二乘解：正规方程与 numpy 的配合

2.1 先从理论说起：为什么 Ax=b 无解时反而有“最优解”

超定方程组就是方程个数大于未知数个数的线性方程组。比如题目给出的第一组：

• 2x1 + 4x2 = 11
• 3x1 - 5x2 = 3
• x1 + 2x2 = 6
• 2x1 + x2 = 7

4 个方程、2 个未知数，绝大多数情况下不存在一组 x 让四个等式同时成立。这在测量、拟合、反演问题里是常态——数据多了，系统就“过定”了。数值分析里处理它的标准思路不是硬解，而是找一个让所有方程残差平方和最小的解，也就是最小二乘解。

关键定理：x* 是 Ax=b 的最小二乘解的充分必要条件是 x* 满足正规方程 A^T A x = A^T b。推导不复杂：目标函数是 ||Ax - b||²，对 x 求导并令梯度为零，就得到这个结论。所以代码实现的路径非常固定：先算 A^T A 和 A^T b，再解一个 n 元线性方程组。注意这里有个“小坑”藏在细节里：A^T A 是方阵，多数情况下是可逆的，但条件数可能很差，后面避坑章节我会专门讲。

2.2 numpy 实现：四行代码求解两个方程组

报告里的源码用 numpy 把整个过程写得很干净。我把它拆开讲：

import numpy as np # 超定方程组 a A = np.array([[2, 4], [3, -5], [1, 2], [2, 1]]) # 系数矩阵 A，4x2 b = np.array([[11], [3], [6], [7]]) # 右端向量 b，4x1 # 构造正规方程 ATA = np.dot(A.T, A) # A^T A，得到 2x2 ATb = np.dot(A.T, b) # A^T b，得到 2x1 # 求解正规方程 print(np.linalg.solve(ATA, ATb))

逻辑说明：np.dot(A.T, A) 先做矩阵乘法得到 2x2 的系数矩阵，np.dot(A.T, b) 得到 2x1 的右端项，最后 np.linalg.solve 用高斯消元解这个 2 元方程组。这里必须注意 A.T 的维度——A 是 4x2，转置后是 2x4，dot 的结果才是 2x2；如果顺序写反，numpy 会直接抛维度不匹配的异常。

第二个方程组也完全同一个套路，只是 A 变成了 5x3 的矩阵：

A1 = np.array([[1, -2, 4], [1, -1, 1], [1, 0, 0], [1, 1, 1], [1, 2, 4]]) # 5x3 系数矩阵 b2 = np.array([[0], [1], [2], [1], [0]]) # 5x1 右端向量 x1 = np.dot(A1.T, A1) b3 = np.dot(A1.T, b2) print(np.linalg.solve(x1, b3))

参数说明：A1 的行对应 5 个方程，列对应 x1、x2、x3 三个未知数。这个结构很典型——列数等于未知数个数，行数等于数据点数。运行结果分别是：

• 方程组 a 的最小二乘解：x1=3.04029304，x2=1.24175824
• 方程组 b 的最小二乘解：x1=1.65714286，x2=0，x3=-0.42857143

2.3 结果验证：别只信 solve 的输出

拿到解之后，我习惯立刻做一个残差验证，这是很多人会跳过的步骤。把解代回原方程组，看每个方程的残差：

x_opt = np.linalg.solve(ATA, ATb) # 上一步解出的最优参数 residual = np.dot(A, x_opt) - b # 原始方程组的残差向量 print(residual.reshape(-1)) # 观察每个方程的误差

逻辑说明：residual 是 4x1 的向量，每个分量对应一个方程的残差。最小二乘解不保证某个残差为零，但会让所有残差的平方和最小。验证的意义在于：如果某个残差特别大（比如比其他大一个数量级），说明原始数据里可能存在异常点，或者你的正规方程条件数很差，解并不可靠。这是“数值解”和“解析解”之间微妙的地方——计算机给的答案，你得自己判断该不该信。

3. 最小二乘拟合实战：数据变换让非线性回归变成线性回归

3.1 题目本质：形如 y = c0 + c1/x 的拟合

题目六给了一组观察数据：x 取 1、2、4、5，对应 y 为 0.33、0.40、0.44、0.45。要求找一个形如 y = c0 + c1/x 的函数做最小二乘拟合。这里有个关键点：y 关于参数 c0、c1 是线性的，但关于自变量 x 不是线性的——它是个非线性函数。数值分析里管这类问题叫“可线性化的非线性拟合”，处理手法非常经典：做变量替换。

令 X' = 1/x，Y' = 1/y，原函数 y = c0 + c1/x 两边取倒数，就变成 1/y = 1/(c0 + c1/x)。这里报告的做法是直接把数据点 (x, y) 变换成 (1/x, 1/y)，然后对变换后的点做线性拟合。换句话说，模型变成 Y' = a·X' + b，其中 X' = 1/x，Y' = 1/y，待求的 a、b 就是线性模型的斜率和截距。这个变换的思路是：不改变参数的最小二乘意义，但把问题从非线性回归降维成线性回归，可以直接套用标准工具。

3.2 sklearn 实现：LinearRegression 的完整用法

import numpy as np from sklearn.linear_model import LinearRegression import matplotlib.pyplot as plt # 原始观察数据：[x, y] X = np.array([[1, 0.33], [2, 0.40], [4, 0.44], [5, 0.45]]) X_n = 1 / X # 对所有元素取倒数，变成 [1/x, 1/y] x_train = X_n[:, :1] # 取第一列作为训练特征（1/x） y_train = X_n[:, 1:2] # 取第二列作为目标值（1/y） a = np.array(x_train).reshape(-1) # 转成一维数组，方便画图 b = np.array(y_train).reshape(-1) # 转成一维数组，方便画图 lineModel = LinearRegression() lineModel.fit(x_train, y_train) # 拟合线性模型 y_predict = lineModel.predict(x_train) c1 = lineModel.coef_[0][0] # 斜率 c0 = lineModel.intercept_[0] # 截距 print("1/y=%.4f 1/x+%.4f" % (c1, c0)) # 画图显示 plt.scatter(a, b, c="blue") # 原始变换后的散点 plt.plot(x_train, y_predict, c="red") # 拟合直线 plt.xlabel('1/x') plt.ylabel('1/y') plt.show()

逻辑说明：X_n = 1 / X 利用 numpy 的广播机制对矩阵每个元素取倒数，一步完成两个变量的变换。x_train 和 y_train 保持二维列向量形状，这是 sklearn 的 LinearRegression 接口的要求——特征必须是二维数组，哪怕只有一列。fit 之后，coef_ 是斜率数组，intercept_ 是截距，两者就是变换后线性模型的参数。

运行结果：斜率 c1=1.00931，截距 c0=2.0143。还原到原函数就是 1/y = 1.00931/x + 2.0143，也就是说 y = 1 / (1.00931/x + 2.0143)。注意这里有个容易搞混的地方——sklearn 拟合的是变换后的空间，不是原空间。如果想把参数带回原模型，必须做逆变换，不能直接拿 c0、c1 当原函数的系数用。

3.3 为什么数据变换要小心：损失函数的意义变了

这里有一个值得反复琢磨的点：原始最小二乘目标是最小化 Σ(y_i - f(x_i))²，但变换后再拟合，最小化的是 Σ(1/y_i - 拟合值)²。两个目标函数不一样，得到的参数也不一样。严格来说，报告里的做法是一种工程近似——它假设变换后的残差分布更接近高斯，实际处理中这种“先变换再拟合”是常见手段，但如果数据里 y 接近零，1/y 会爆炸，变换就翻车了。更好的做法是用 scipy 的 curve_fit 直接拟合原函数，后面我会给出替代方案。

4. 石油产量预测：四种多项式拟合与残差对比

4.1 任务描述与数据

实验四的数据是 1994 到 2003 十年间全球石油日产量（百万桶/天），要求分别用直线、抛物线、立方曲线拟合这 10 个数据点，写出具体的函数表达式，做残差分析比较哪种拟合最好，再外推预测 2010 年的产量。

数据本身有轻微波动但整体上行，从 67.052 涨到 76.777。题目问“哪一种拟合最好”——这不是看谁曲线贴得近，而是看谁在“预测”这个任务上更可信。报告里用三个指标衡量：二范数的平方（即残差平方和）、残差的标准差、残差绝对值的均值。这三个指标从不同侧面描述拟合误差：二范数平方放大大误差，标准差体现离散程度，绝对均值给一个直观的平均偏差。

4.2 polyfit 与手动定义函数：两种写法的对比

报告代码里同时出现了两种拟合实现：np.polyfit 和 scipy.optimize.leastsq。前者是多项式拟合的快捷方式，后者是通用最小二乘参数估计。两种写法各有适用场景：

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import leastsq x = np.arange(0, 10, 1, dtype='float') # 年份偏移量：0 代表 1994，9 代表 2003 y = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777], dtype='float') # 写法一：polyfit 直接拟合多项式 an = np.polyfit(x, y, 1) # 1 次多项式（直线） print("2010年预测：%.4f" % (np.polyval(an, 2010 - 1994))) p1 = np.poly1d(an) print(p1) # 打印多项式表达式 # 写法二：leastsq 通用拟合 param0 = [0, 0] def linear_fun(s, x): k1, b = s return k1 * x + b def dist(a, fun, x, y): return fun(a, x) - y var = leastsq(dist, param0, args=(linear_fun, x, y)) print('[linear] 二范数：%.4f' % ((y - linear_fun(var[0], x))**2).sum())

逻辑说明：polyfit 返回系数数组，poly1d 把系数包装成可直接调用的多项式对象；leastsq 则需要你手动定义函数形式，传入初始参数猜测和残差函数。两者解出来的参数在数值上几乎一致，区别在于 polyfit 只适合多项式，leastsq 可以拟合任意自定义函数。param0 是初始值，这里给 [0, 0] 在简单线性问题上没问题，但如果拟合更复杂的模型，初始值选不好 leastsq 可能收敛到局部最优。

4.3 四种拟合的残差对比：谁更好不能只看图

报告的表 1、表 2 给出了明确的数值对比：

从数值上看，立方曲线的三个指标都是最小的，拟合精度最高。这符合直观预期——数据本身有波动，多项式次数越高，越能捕捉局部变化。但注意，拟合精度高不等于预测可靠。表 3 的 2010 年预测值暴露了问题：

抛物线预测值比 2003 年还低，出现了下降趋势，和实际经济发展趋势明显不符；立方曲线给出 99.04，增长幅度远大于前面十年的趋势。为什么立方曲线拟合误差最小、外推却最离谱？这就是高次多项式外推的经典陷阱——多项式在数据区间外会按照最高次项主导的方向飞速增长，次数越高外推越猛。报告结论选择了立方曲线，理由是“结合经济社会发展趋势”，但从纯数学角度看，外推预测本身就带有很强的主观性。

4.4 画图技巧：坐标偏移的细节

报告代码里有一个小技巧值得注意：x 用的是 0 到 9 的偏移量，但画图时用 x+1994 把横轴还原成真实年份。plt.axis([1994, 2003, min(y)-0.5, max(y)+0.5]) 设置坐标范围，plt.grid(True) 打开网格方便读值。预测 2010 年时用 2010-1994=16 作为输入，这是个容易出错的细节——如果你直接用年份 2010 代入 polyfit 得到的多项式，结果会完全不对，因为拟合时 x 的基准是偏移量而不是真实年份。

5. 常见问题与避坑：本实验最容易翻车的五个点

5.1 plt.axis 范围设置导致图形空白

现象：代码设置了 plt.axis([1994, 2003, min(y)-0.5, max(y)+0.5])，但画出来是空白图。

原因：leastsq 循环里 plt.plot(x+1994, ...) 用的是偏移量加年份，而如果你在 plt.show() 之前先调用了 plt.axis，坐标轴范围会被锁定，而后面的 plot 数据超出范围时就只能看到局部。另一个更隐蔽的原因是：axis 函数接收的参数顺序是 [xmin, xmax, ymin, ymax]，如果写反了，图形会直接错乱。

解决：先 plot 完所有曲线，最后再设置 axis；或者干脆不设置，让 matplotlib 自动缩放。

5.2 中文字体 SimHei 不生效

现象：plt.title 里的中文显示成方块，负号显示成方块。

原因：matplotlib 默认字体不支持中文，负号用的是 Unicode 的减号，默认字体渲染不出来。代码里已经给出了标准解法：mpl.rcParams['font.sans-serif'] = ['SimHei'] 指定黑体，mpl.rcParams['axes.unicode_minus'] = False 让负号用 ASCII 的 '-'。

解决：如果换成 Mac 或 Linux，SimHei 可能不存在，需要改成 'WenQuanYi Zen Hei' 或 'Noto Sans CJK SC'。建议用 plt.rcParams['font.family'] = 'sans-serif' 加一个系统中文字体列表兜底。

5.3 用年份直接拟合多项式，预测值错得离谱

现象：用 np.polyfit(x, y, 1) 拟合时，x 直接传真实年份，得到的多项式系数大得吓人，预测 2010 年输出一个完全不符合量级的值。

原因：真实年份是 1994 到 2003，数值很大，多项式拟合时高次项会产生严重的数值不稳定。np.polyfit 内部用的是以 x 的幂为基的范德蒙德矩阵，x 太大条件数爆炸。报告中用 x = np.arange(0, 10) 做偏移就是为了规避这个问题。

解决：所有拟合和预测都在偏移后的坐标下进行。如果非要显示真实年份，只在画图时加偏移量。

5.4 leastsq 初始参数给太差，拟合结果不稳定

现象：换一组初始参数，leastsq 拟合出来的结果差别很大，甚至报错。

原因：leastsq 用的是 Levenberg-Marquardt 算法，对初始值敏感。题目里的初始值给的是 [0, 0]，对直线和抛物线问题足够了，但如果拟合更复杂的函数，比如带指数项或对数项的模型，初始值离真实解太远会收敛到局部极小值。

解决：先用 polyfit 求一个低阶近似的参数作为初始值，再传给 leastsq；或者改用 scipy.optimize.curve_fit，它对参数绑定的检查更友好，报错信息也更明确。

5.5 残差分析指标选错，对比失去意义

现象：直接比较不同拟合方式的“误差平方和”，得出高次多项式永远更好的结论。

原因：残差平方和会随着多项式次数升高单调下降——次数越高，拟合能力越强，9 次多项式甚至能精确穿过所有数据点，但预测能力反而最差。这就是过拟合。报告里的三个指标只能衡量“拟合程度”，不能衡量“预测能力”，所以得到立方曲线最好这个结论本身就有前提限制。

解决：要做真正的模型对比，需要用交叉验证或留一法，把数据分成训练集和测试集，看测试集上的预测误差。后面改造章节会给一个最小实现。

6. 二次开发思路：把验证型脚本改成可复用的拟合工具

前四个实验脚本的价值在于演示，但如果要在自己的数值分析作业或真实项目里复用，代码还差点意思——所有数据都硬编码在数组里、模型函数写死、评估逻辑散落在主流程中。我一般会先做一次小重构，把“拟合”和“评估”抽成通用函数。

核心思路是封装一个拟合器类，支持任意多项式次数，并自动完成残差指标计算：

import numpy as np from scipy.optimize import leastsq class PolyFitter: def __init__(self, x, y, degree): self.x = np.asarray(x, dtype=float) self.y = np.asarray(y, dtype=float) self.degree = degree self.params = None def _model(self, coef, x): return np.polyval(coef, x) def _residual(self, coef): return self._model(coef, self.x) - self.y def fit(self, init=None): if init is None: init = np.zeros(self.degree + 1) # 初始参数全零 self.params = leastsq(self._residual, init)[0] return self def predict(self, x): return np.polyval(self.params, x) def metrics(self): res = self._residual(self.params) return { "l2_sq": (res**2).sum(), "std": res.std(), "abs_mean": np.abs(res).mean(), } x = np.arange(1994, 2004, dtype=float) y = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777]) fitter = PolyFitter(x - 1994, y, degree=3).fit() print(fitter.params) print(fitter.metrics()) print("2010 预测：", fitter.predict(2010 - 1994))

逻辑说明：init接收原始数据，fit 里用 leastsq 优化多项式系数，predict 用 polyval 做预测，metrics 一次性算出三个评估指标。这样改完，换新数据、换多项式次数都只需要改一行代码。

再进一步，如果想把数据源换成文件而不是写死的数组，就加一个 csv 读取函数；如果想做模型对比，就把不同 degree 的 PolyFitter 放进一个循环里，分别算指标并输出对比表。最后一件事是验证：对石油产量这类时间序列，我建议至少做一次留一法交叉验证——每次拿 9 个点拟合，预测剩下的 1 个点，统计预测误差。你会发现结果和三指标排名很可能不一样。从那以后，我每次拿到拟合任务都会强制走一遍“画散点图 → 试不同次数 → 交叉验证”的流程再下结论，这个习惯帮我避开了不少“拟合效果好但预测全错”的坑。希望帮到你。

本文还有配套的精品资源，点击获取