简介广东工业大学数值分析课程实验四的Python实验报告完整覆盖超定方程组最小二乘解与最小二乘拟合两大核心实验。报告为docx文档共1个文件、112KB包含题目三、题目六及石油产量预测三道典型问题均附有可运行的Python源码、算法介绍和结果分析。题目三基于最小二乘解的充要条件构造正规方程组ATAxATb借助numpy的线性代数求解器高效算出两组超定方程的解题目六以1和1/x为基函数建立线性组合模型调用sklearn的LinearRegression拟合观察数据得到形如y2.01431.00931/x的拟合函数石油产量预测则演示了将数值方法应用于实际数据建模的完整流程。读者可从中学习最小二乘问题的矩阵推导、正规方程构造、Python库调用及结果分析方法适合正在学习数值分析、需要参考完整实验报告或复习相关算法的高校学生。目前已有212人学习该资源。1. 某高校Python数值分析实验四为什么代码三分钟跑完作业却要写一整晚数值分析这门课真正动手做起来和听课时完全是两个世界。某高校的Python数值分析实验四表面上是插值和拟合两个经典主题——拉格朗日插值、牛顿插值、最小二乘拟合代码用NumPy就能写完跑通不过三分钟。但真正消耗时间的是把课本上的公式变成能交代得清楚的程序把结果画成能放进实验报告的图表以及处理那些课本没写、IDE也不会告诉你的边界情况。这篇笔记按我实际完成实验四的顺序来写从任务拆解、环境准备、算法实现到避坑排查每一步都给可复现的代码和参数说明。适合第一次做数值分析实验的学生也适合想快速拾起插值与拟合实现的工程师。2. 拆开实验四的作业要求任务清单、环境准备与报告结构2.1 先搞清实验四考什么两个主线任务和三类题型数值分析实验四在多数高校的教学安排里主题都落在“插值法与曲线拟合”这一章。拆开来看主线任务基本是固定的用Python实现拉格朗日插值和牛顿插值再用最小二乘法做曲线拟合最后把结果和理论做对比。不同年份、不同老师的题目形式会有变化但底层要求逃不出三类题型。第一类是给定函数或数据表要求计算指定点的插值结果第二类是给一组观测数据要求用多项式拟合出趋势第三类是综合题比如对某个病态函数做插值并讨论误差或者用拟合预测下一组数据点。先确认自己拿到的是哪类题再决定代码怎么组织能省下大量返工时间。我见过不少同学上来就写拉格朗日的循环写完才发现题目要求的是拟合。更常见的坑是题目给了十个数据点但要求分别用三次、五次、七次多项式做拟合这时候如果没提前把拟合逻辑封装成函数后面每一个次数都要复制粘贴改参数实验报告的代码部分会显得很乱查错也麻烦。所以第一步不是写代码是把题目要求翻译成一个清单要实现的算法有几个、输入输出是什么、要求画几张图、误差分析要算哪些量。2.2 实验环境准备conda创建虚拟环境与依赖清单实验四的代码依赖面不宽但版本不一致会在后面图表上翻车。我一般会先建一个干净的虚拟环境避免和机器上其他项目共用NumPy和SciPy。用conda创建虚拟环境的命令如下conda create -n numeric-analysis python3.11 conda activate numeric-analysis pip install numpy matplotlib scipy选Python 3.11而不是最新版本是因为NumPy和SciPy对最新Python版本的支持往往滞后几个月装完才发现某个依赖编译失败再回退更浪费时间。NumPy负责数组计算和多项式运算SciPy在拟合部分会用到curve_fitMatplotlib负责画图。这三个库覆盖面足够pandas在这个实验里用不上装不装看个人习惯。依赖安装完以后建议顺手验证一下版本号是否匹配尤其是NumPy和SciPy之间的版本兼容。我遇到过NumPy 1.x和SciPy 1.11搭配时某个接口行为变化导致polyfit返回结果格式不同。验证方式很简单在Python交互环境里分别import并打印版本号确认能正常加载再开始写实验代码。如果用的是机房预装的环境版本往往偏旧建议先确认能不能装新包不能装就基于现有版本调整写法别到最后一刻才暴露问题。2.3 实验报告的结构算法原理、程序实现、结果分析三位一体实验报告怎么写决定了这份作业的最终成绩。数值分析实验的评分点通常有三个算法原理是否说清楚、代码是否能跑出正确结果、结果分析是否有深度。三者缺一不可但很多同学只重视代码原理部分抄课本分析部分写“图像符合预期”这等于把最容易拿分的部分放弃了。我的习惯是先在报告里画一个整体框架再按照框架去补内容。第一部分写算法原理包括公式推导和算法步骤不需要长篇大论但要写清楚公式中每个符号的含义比如插值节点、基函数、差商。第二部分写程序实现代码按功能拆分关键函数要给注释。第三部分写结果分析这一部分要包含数值结果和图像的对比以及误差分析——插值法算出的值和理论值差多少、拟合的残差有多大、为什么会有这个误差。整份报告的逻辑是读者只看报告不看你跑代码也能复现你的结果这才是一份合格的实验报告。3. 拉格朗日插值与牛顿插值从公式到Python实现3.1 拉格朗日插值法的数学原理与代码实现插值问题的本质是已知n1个数据点找一条通过所有这些点的曲线。拉格朗日插值的思想最直白构造一组基函数每个基函数在某个数据点上取值为1在其他数据点上取值为0然后把所有基函数加权叠加。权重就是对应数据点的函数值。这样做的好处是公式对称、好理解代码实现也不需要额外存储中间量。import numpy as np def lagrange_interp(x_points, y_points, x_eval): 拉格朗日插值计算单个点 x_eval 处的插值结果 x_points: 插值节点数组例如 [1, 2, 3] y_points: 对应的函数值数组例如 [1, 4, 9] x_eval: 需要求插值的点可以是标量或数组 n len(x_points) if n 0: raise ValueError(插值节点不能为空) total 0.0 for i in range(n): term y_points[i] for j in range(n): if j ! i: term * (x_eval - x_points[j]) / (x_points[i] - x_points[j]) total term return total这段代码的逻辑是两层循环外层遍历每个数据点内层累乘计算对应基函数在当前求值点x_eval处的值。内层循环里分母是x_points[i] - x_points[j]当x_eval恰好等于某个插值节点时对应项的分母会精确抵消成一其他项因为分子出现零而整个为0最终返回该节点的函数值这是拉格朗日插值最基本的性质。代码里加了个空节点检查防止传入空数组时报出难以理解的numpy错误。用的时候注意两个参数语义x_points是插值节点而非插值区间端点它决定多项式的次数n个节点对应n-1次多项式。x_eval可以是单个数值也可以是numpy数组函数内部会自动做广播。实测中对数组求值时这个函数会在每个元素上独立执行两层循环数据量上千点时速度会明显下降所以如果只需要在特定几个点求值不要传一个很长的密集数组这属于代码效率和精度无关但和体验相关的问题。3.2 牛顿插值法差商表与递推实现牛顿插值和拉格朗日插值数学上等价但实现思路不同。牛顿插值通过构造差商表来逐步构建多项式它的优势在于当新增一个插值节点时不需要重写所有基函数只需在原有多项式上加一项。这个性质在实验里体现得最明显的是题目要求先取前三个点插值看结果再取前四个点插值看结果用牛顿插值可以复用之前的计算。def newton_divided_differences(x_points, y_points): 计算差商表返回对角线上的差商系数数组 n len(x_points) coef np.array(y_points, dtypefloat) for j in range(1, n): for i in range(n - 1, j - 1, -1): coef[i] (coef[i] - coef[i - 1]) / (x_points[i] - x_points[i - j]) return coef def newton_interp(coef, x_points, x_eval): 用差商系数计算牛顿插值多项式在 x_eval 处的值 coef: newton_divided_differences 的返回值 n len(coef) result coef[0] for i in range(1, n): term coef[i] for j in range(i): term * (x_eval - x_points[j]) result term return result差商表计算是原地更新的coef数组在每一轮外层循环中被覆盖最终coef[i]存放的是第i阶差商。注意内层循环从后向前遍历这保证计算新值时用到的是上一轮的结果而不是本轮已经覆盖的值顺序反了整个差商表就是错的。牛顿插值的求值过程用到了嵌套乘法从最高次项开始逐层提取公因子这样比直接展开多项式运算量更小也更符合数值分析教材里对秦九韶算法的介绍。实操中差商表函数只在构造阶段调用一次之后对多个x_eval求值时只需要反复调用newton_interp。如果你要插值的点有几百上千个这个做法的优势非常明显——差商构造是O(n²)而每个点求值是O(n)。拉格朗日插值对每个求值点都要重新计算所有基函数同样是O(n²)但乘在点的数量上数据量一大差距就拉开了。3.3 用matplotlib画出插值效果对比图实验报告里最直观的内容是一张插值效果图原始数据点、拉格朗日插值曲线、牛顿插值曲线在同一坐标系中展示。画图代码不长但有两个细节要注意插值曲线的横坐标要足够密不然曲线会显得很生硬纵坐标范围要留出余量避免曲线在边缘处超出画布。import matplotlib.pyplot as plt x_points np.array([0.0, 1.0, 2.0, 3.0, 4.0]) y_points np.array([1.0, 2.0, 1.5, 3.0, 2.8]) x_dense np.linspace(0, 4, 200) fig, ax plt.subplots(figsize(8, 5)) ax.plot(x_points, y_points, o, label原始数据点) ax.plot(x_dense, lagrange_interp(x_points, y_points, x_dense), -, label拉格朗日插值) coef newton_divided_differences(x_points, y_points) ax.plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛顿插值) ax.set_xlabel(x) ax.set_ylabel(y) ax.legend() ax.grid(True, linestyle--, alpha0.6) plt.show()画图时x_dense用np.linspace(0, 4, 200)生成200个均匀分布的点具体点数取决于插值多项式的次数和区间跨度一般200到500足够。如果两条插值曲线完全重合说明两个算法在精度上一致这是正常的如果出现细微差异优先检查差商表的迭代顺序是否正确。对比图存在的意义是验证两个算法结果一致顺带检查曲线的形状是否符合直觉。4. 最小二乘曲线拟合把公式拟合到数据点上的完整流程4.1 最小二乘原理与numpy.polyfit最小实现拟合和插值的核心区别在于插值要求曲线严格穿过每个数据点拟合只要求曲线在整体上最接近数据点。最小二乘的思想是找一组参数让模型计算值与观测值的平方误差之和最小。平方误差比绝对误差更敏感而且方便求导所以成为最常用的拟合准则。import numpy as np x_data np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]) y_data np.array([0.1, 0.8, 1.8, 2.9, 5.1, 7.2, 10.0]) degree 2 coeffs np.polyfit(x_data, y_data, degree) p np.poly1d(coeffs) print(拟合系数, coeffs) print(拟合多项式, p)np.polyfit返回的是多项式系数从最高次到常数项排列系数顺序和poly1d的输入是对应的。degree参数是最关键的调试项取1时是线性拟合取2时是抛物线拟合取大了会出现过拟合——曲线为了贴近个别数据点而剧烈抖动尤其是在数据本身不是高阶多项式的情况下。判断次数是否合适的标准很简单看拟合曲线走势是否符合数据的物理含义再看残差是否随机分布而不是呈现系统性偏差。poly1d是把系数数组包装成可以直接调用的多项式对象调用p(x_eval)就能批量求值。这里有个小细节polyfit内部用的是基于SVD的稳定算法对于次数不超过8的多项式拟合很少出现数值问题但超过10次就要警惕矩阵病态结果可能不可信。所以实验题里如果要求拟合次数特别高建议检查一下拟合系数是否出现数量级异常的大值那通常是条件数过大导致的数值性翻车。4.2 自定义函数拟合scipy.optimize.curve_fit的使用实验题里有时候不要求多项式拟合而是指定一个函数形式比如指数衰减模型或者带偏移的正弦函数。np.polyfit只支持多项式这时候要换用scipy的curve_fit它可以拟合任意你写出的Python函数。核心用法是先定义一个函数形式函数的第一个参数是自变量后面是待拟合参数。from scipy.optimize import curve_fit def exp_model(x, a, b, c): 指数衰减加常数偏移的拟合模型 return a * np.exp(-b * x) c popt, pcov curve_fit(exp_model, x_data, y_data, p0[1.0, 1.0, 0.0]) a_opt, b_opt, c_opt popt print(f拟合参数: a{a_opt:.4f}, b{b_opt:.4f}, c{c_opt:.4f}) x_smooth np.linspace(x_data.min(), x_data.max(), 100) y_fit exp_model(x_smooth, *popt)p0是参数初值这个参数在很多时候决定了拟合能不能收敛。初值差得太远算法可能陷入局部极小值甚至直接报错。给p0的经验法则是先看数据点的取值范围大致估计一下参数的量级。比如数据y值在0到10之间c的初值可以设为数据中点或最小值附近a的初值设为y的跨度b的初值根据曲线衰减的快慢估一个正数。如果你发现拟合结果明显不合理最先该改的是p0而不是模型本身。curve_fit返回的pcov是参数协方差矩阵对角线元素开根号是参数的标准差代表拟合参数的不确定度。实验报告里如果想体现严谨性可以把参数误差也打印出来。当协方差矩阵里的值大到离谱时往往说明模型参数之间存在强相关性或者数据量不足以确定这么多参数这时候的拟合系数虽然能算出数但不能当作物理解释使用。这个点算是一个隐藏的知识点写在报告的分析部分会很加分。4.3 拟合效果的量化评估残差、R方与可视化拟合做完不评估就写报告等于黑匣子输出直接交差。评估拟合质量最常用的是三个量残差、均方根误差和R方。残差是每个数据点的观测值减拟合值R方表示拟合曲线解释了数据方差的百分比越接近1越好。residuals y_data - exp_model(x_data, *popt) ss_res np.sum(residuals**2) ss_tot np.sum((y_data - np.mean(y_data))**2) r_squared 1 - ss_res / ss_tot rmse np.sqrt(np.mean(residuals**2)) print(fR方 {r_squared:.4f}) print(f均方根误差 RMSE {rmse:.4f})注意R方有一个陷阱当模型形式比较复杂、参数很多时R方会虚高甚至接近1但这不代表模型真的好。一个十次多项式拟合十个数据点R方必然是1但那完全没有预测意义。所以评估时要同时看残差图——如果残差随x增大而增大说明模型形式不对如果残差呈现波浪形说明漏掉了周期性分量。这些判断是数值分析实验最有价值的训练比单纯算对拉格朗日插值结果重要得多。把残差图画出来和拟合效果图画在一起报告的分析部分立刻有了深度。5. 实验四避坑指南5个最常见的翻车点与排查方法5.1 拉格朗日插值次数过高导致龙格现象数据越多结果越离谱现象用拉格朗日插值去逼近某些连续函数插值节点增加到十几个甚至更多插值曲线在区间两端剧烈震荡误差反而比少用几个节点时更大。原因这是数值分析教材里反复强调的龙格现象。拉格朗日插值多项式在等距节点下对龙格函数、绝对值函数这类目标做高次逼近时插值余项在端点附近会放大多项式次数越高震荡越明显。解决明确实验四的题目通常只要求4到6个节点的插值不需要追求高次。如果题目要求对比不同节点数下的误差把插值点数从3逐个增加到10画误差曲线让读者直观看到误差先降后升的过程这本身就是一份合格的现象分析。如果实际工程中遇到类似需求应该改用分段插值或样条插值而不是硬上高次多项式。5.2 图表的横纵坐标出现方块字中文标注全部乱码现象matplotlib画图时标题和图例里的中文显示成一个个小方块英文和数字正常实验报告截图里一眼就能看出问题。原因matplotlib默认字体不包含中文字形需要显式指定中文字体。Windows系统一般有SimHei或Microsoft YaHeiLinux和macOS的字体名不同用系统自带字体名逐个尝试。解决在画图代码最前面加上字体设置让matplotlib全局使用中文字体如下所示import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, PingFang SC] plt.rcParams[axes.unicode_minus] False第二行去掉坐标轴负号的特殊处理是因为切换了中文字体之后负号有可能显示成一个奇怪的方框。这个问题属于程序能跑但报告不能看的隐形翻车查错成本高提前设置好就能完全避免。5.3 差商表计算顺序写反牛顿插值结果对不上拉格朗日现象牛顿插值和拉格朗日插值结果不一致差距可能很小也可能很大越小规模的数据差距越小很多同学因此没发现。原因newton_divided_differences函数里内层循环的方向错了。差商表必须从后向前原地更新如果写成从前向后计算新的差商时用的是已经被本轮覆盖过的值导致差商全部错误。数据点少时误差不明显数据点多时完全对不上。解决内层循环固定使用for i in range(n - 1, j - 1, -1)从后往前遍历。排查思路也很简单用三个点的手工差商计算结果和Python输出对比如果差商对不上问题就在循环方向。这类错误的特点是代码能跑、结果贴近但不对最容易在实验报告的误差分析部分暴露靠直觉很难发现。5.4 拟合结果出现巨大系数多项式次数太高导致矩阵病态现象polyfit指定degree9或更高时拟合系数出现10的6次方以上的大数拟合曲线在数据点之间剧烈波动部分区间甚至超出数据范围好几倍。原因高次多项式拟合的等效设计矩阵条件数随次数指数增长数据点x跨度大时尤其严重。矩阵接近奇异SVD分解得到的解被微小舍入误差放大系数自然不稳定。解决首先判断题目是否真的要求高次拟合。如果要求做不同次数拟合的对比把次数范围控制在1到6之间重点观察拟合误差变化趋势。如果必须拟合高次多项式考虑先对x数据做归一化把数据平移到零均值并缩放到单位方差再拟合最后把系数还原回原始坐标。归一化这步在实验报告里是加分项能体现你对数值稳定性的理解。5.5 求值点落在插值节点范围之外结果完全失控现象用插值函数计算区间外的点计算结果非常大和真实函数值差出几个数量级有的同学以为是代码精度不够。原因插值多项式只在节点覆盖的区间内有可靠逼近能力。区间外没有任何数据点约束多项式形状高次项起主导作用结果必然外溢。这是插值方法的固有属性不是bug。解决求值前先检查x_eval是否在x_points的范围内。如果题目故意要求计算外推值那么在报告里要明确记录这个现象用误差对比说明外推不可靠这属于正面的实验结果。如果在做拟合预测用polyfit拟合出的多项式做短距离预测是可行的做远距离预测同样会失效。边界意识是从数值分析实验里能带走的最重要工程习惯。6. 用一份可视化脚本把实验结果整理成答辩图表实验四答辩或提交报告时图表质量往往直接影响印象分。同样是解决同一个问题一个用默认绘图风格、图例位置遮挡曲线、坐标轴没标单位和另一个图面干净、曲线粗细分明、每个子图都有标题和注释给阅卷人的感受完全不同。整理图表这一步建议独立写一个脚本不要和算法代码混在一起输出为高分辨率图片再插入报告而不是直接截图。import matplotlib.pyplot as plt import numpy as np fig, axes plt.subplots(1, 3, figsize(14, 4)) # 第一个子图插值结果对比 axes[0].plot(x_points, y_points, o, label数据点) axes[0].plot(x_dense, lagrange_interp(x_points, y_points, x_dense), label拉格朗日) coef newton_divided_differences(x_points, y_points) axes[0].plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛顿) axes[0].set_title(插值结果对比) axes[0].legend() # 第二个子图拟合效果 axes[1].plot(x_data, y_data, o, label观测数据) axes[1].plot(x_smooth, y_fit, -, label拟合曲线) axes[1].set_title(最小二乘拟合) # 第三个子图残差分布 axes[2].plot(x_data, residuals, o-, label残差) axes[2].axhline(0, colorgray, linestyle--) axes[2].set_title(拟合残差分布) axes[2].legend() fig.tight_layout() plt.savefig(experiment4_summary.png, dpi200, bbox_inchestight) plt.show()一处细节是plt.savefig里的参数dpi至少200才够打印清晰bbox_inchestight可以把空白边缘裁掉避免图片四周留白过多。如果你是准备提交电子版报告300dpi更好图片文件大一点没关系打印效果匀称。另一个细节是每个子图的标题要能自解释也就是说读者不看正文只看图就知道这张图在表达什么这是实验报告图表规范的基本要求多数同学没意识到所以认真做反而会留下印象。我自己的习惯是在画图脚本开头放一段注释记录这份图表对应的题目要求、数据文件、生成日期方便复查时快速定位。数值分析实验四做完回头看最大的收获往往不是那几个算法本身而是理解了“数据点给出的信息量是有限的任何试图从有限数据外推的做法都必须说明假设和边界”。这个意识在后面的课程和工程项目里反复被验证。希望这篇笔记能帮你把实验四从头到尾走通拿到一个能进报告、能讲清楚、也能经得起提问的结果。本文还有配套的精品资源点击获取 Python数值分析实验四：拉格朗日插值、牛顿插值与最小二乘拟合实践指南
简介广东工业大学数值分析课程实验四的Python实验报告完整覆盖超定方程组最小二乘解与最小二乘拟合两大核心实验。报告为docx文档共1个文件、112KB包含题目三、题目六及石油产量预测三道典型问题均附有可运行的Python源码、算法介绍和结果分析。题目三基于最小二乘解的充要条件构造正规方程组ATAxATb借助numpy的线性代数求解器高效算出两组超定方程的解题目六以1和1/x为基函数建立线性组合模型调用sklearn的LinearRegression拟合观察数据得到形如y2.01431.00931/x的拟合函数石油产量预测则演示了将数值方法应用于实际数据建模的完整流程。读者可从中学习最小二乘问题的矩阵推导、正规方程构造、Python库调用及结果分析方法适合正在学习数值分析、需要参考完整实验报告或复习相关算法的高校学生。目前已有212人学习该资源。1. 某高校Python数值分析实验四为什么代码三分钟跑完作业却要写一整晚数值分析这门课真正动手做起来和听课时完全是两个世界。某高校的Python数值分析实验四表面上是插值和拟合两个经典主题——拉格朗日插值、牛顿插值、最小二乘拟合代码用NumPy就能写完跑通不过三分钟。但真正消耗时间的是把课本上的公式变成能交代得清楚的程序把结果画成能放进实验报告的图表以及处理那些课本没写、IDE也不会告诉你的边界情况。这篇笔记按我实际完成实验四的顺序来写从任务拆解、环境准备、算法实现到避坑排查每一步都给可复现的代码和参数说明。适合第一次做数值分析实验的学生也适合想快速拾起插值与拟合实现的工程师。2. 拆开实验四的作业要求任务清单、环境准备与报告结构2.1 先搞清实验四考什么两个主线任务和三类题型数值分析实验四在多数高校的教学安排里主题都落在“插值法与曲线拟合”这一章。拆开来看主线任务基本是固定的用Python实现拉格朗日插值和牛顿插值再用最小二乘法做曲线拟合最后把结果和理论做对比。不同年份、不同老师的题目形式会有变化但底层要求逃不出三类题型。第一类是给定函数或数据表要求计算指定点的插值结果第二类是给一组观测数据要求用多项式拟合出趋势第三类是综合题比如对某个病态函数做插值并讨论误差或者用拟合预测下一组数据点。先确认自己拿到的是哪类题再决定代码怎么组织能省下大量返工时间。我见过不少同学上来就写拉格朗日的循环写完才发现题目要求的是拟合。更常见的坑是题目给了十个数据点但要求分别用三次、五次、七次多项式做拟合这时候如果没提前把拟合逻辑封装成函数后面每一个次数都要复制粘贴改参数实验报告的代码部分会显得很乱查错也麻烦。所以第一步不是写代码是把题目要求翻译成一个清单要实现的算法有几个、输入输出是什么、要求画几张图、误差分析要算哪些量。2.2 实验环境准备conda创建虚拟环境与依赖清单实验四的代码依赖面不宽但版本不一致会在后面图表上翻车。我一般会先建一个干净的虚拟环境避免和机器上其他项目共用NumPy和SciPy。用conda创建虚拟环境的命令如下conda create -n numeric-analysis python3.11 conda activate numeric-analysis pip install numpy matplotlib scipy选Python 3.11而不是最新版本是因为NumPy和SciPy对最新Python版本的支持往往滞后几个月装完才发现某个依赖编译失败再回退更浪费时间。NumPy负责数组计算和多项式运算SciPy在拟合部分会用到curve_fitMatplotlib负责画图。这三个库覆盖面足够pandas在这个实验里用不上装不装看个人习惯。依赖安装完以后建议顺手验证一下版本号是否匹配尤其是NumPy和SciPy之间的版本兼容。我遇到过NumPy 1.x和SciPy 1.11搭配时某个接口行为变化导致polyfit返回结果格式不同。验证方式很简单在Python交互环境里分别import并打印版本号确认能正常加载再开始写实验代码。如果用的是机房预装的环境版本往往偏旧建议先确认能不能装新包不能装就基于现有版本调整写法别到最后一刻才暴露问题。2.3 实验报告的结构算法原理、程序实现、结果分析三位一体实验报告怎么写决定了这份作业的最终成绩。数值分析实验的评分点通常有三个算法原理是否说清楚、代码是否能跑出正确结果、结果分析是否有深度。三者缺一不可但很多同学只重视代码原理部分抄课本分析部分写“图像符合预期”这等于把最容易拿分的部分放弃了。我的习惯是先在报告里画一个整体框架再按照框架去补内容。第一部分写算法原理包括公式推导和算法步骤不需要长篇大论但要写清楚公式中每个符号的含义比如插值节点、基函数、差商。第二部分写程序实现代码按功能拆分关键函数要给注释。第三部分写结果分析这一部分要包含数值结果和图像的对比以及误差分析——插值法算出的值和理论值差多少、拟合的残差有多大、为什么会有这个误差。整份报告的逻辑是读者只看报告不看你跑代码也能复现你的结果这才是一份合格的实验报告。3. 拉格朗日插值与牛顿插值从公式到Python实现3.1 拉格朗日插值法的数学原理与代码实现插值问题的本质是已知n1个数据点找一条通过所有这些点的曲线。拉格朗日插值的思想最直白构造一组基函数每个基函数在某个数据点上取值为1在其他数据点上取值为0然后把所有基函数加权叠加。权重就是对应数据点的函数值。这样做的好处是公式对称、好理解代码实现也不需要额外存储中间量。import numpy as np def lagrange_interp(x_points, y_points, x_eval): 拉格朗日插值计算单个点 x_eval 处的插值结果 x_points: 插值节点数组例如 [1, 2, 3] y_points: 对应的函数值数组例如 [1, 4, 9] x_eval: 需要求插值的点可以是标量或数组 n len(x_points) if n 0: raise ValueError(插值节点不能为空) total 0.0 for i in range(n): term y_points[i] for j in range(n): if j ! i: term * (x_eval - x_points[j]) / (x_points[i] - x_points[j]) total term return total这段代码的逻辑是两层循环外层遍历每个数据点内层累乘计算对应基函数在当前求值点x_eval处的值。内层循环里分母是x_points[i] - x_points[j]当x_eval恰好等于某个插值节点时对应项的分母会精确抵消成一其他项因为分子出现零而整个为0最终返回该节点的函数值这是拉格朗日插值最基本的性质。代码里加了个空节点检查防止传入空数组时报出难以理解的numpy错误。用的时候注意两个参数语义x_points是插值节点而非插值区间端点它决定多项式的次数n个节点对应n-1次多项式。x_eval可以是单个数值也可以是numpy数组函数内部会自动做广播。实测中对数组求值时这个函数会在每个元素上独立执行两层循环数据量上千点时速度会明显下降所以如果只需要在特定几个点求值不要传一个很长的密集数组这属于代码效率和精度无关但和体验相关的问题。3.2 牛顿插值法差商表与递推实现牛顿插值和拉格朗日插值数学上等价但实现思路不同。牛顿插值通过构造差商表来逐步构建多项式它的优势在于当新增一个插值节点时不需要重写所有基函数只需在原有多项式上加一项。这个性质在实验里体现得最明显的是题目要求先取前三个点插值看结果再取前四个点插值看结果用牛顿插值可以复用之前的计算。def newton_divided_differences(x_points, y_points): 计算差商表返回对角线上的差商系数数组 n len(x_points) coef np.array(y_points, dtypefloat) for j in range(1, n): for i in range(n - 1, j - 1, -1): coef[i] (coef[i] - coef[i - 1]) / (x_points[i] - x_points[i - j]) return coef def newton_interp(coef, x_points, x_eval): 用差商系数计算牛顿插值多项式在 x_eval 处的值 coef: newton_divided_differences 的返回值 n len(coef) result coef[0] for i in range(1, n): term coef[i] for j in range(i): term * (x_eval - x_points[j]) result term return result差商表计算是原地更新的coef数组在每一轮外层循环中被覆盖最终coef[i]存放的是第i阶差商。注意内层循环从后向前遍历这保证计算新值时用到的是上一轮的结果而不是本轮已经覆盖的值顺序反了整个差商表就是错的。牛顿插值的求值过程用到了嵌套乘法从最高次项开始逐层提取公因子这样比直接展开多项式运算量更小也更符合数值分析教材里对秦九韶算法的介绍。实操中差商表函数只在构造阶段调用一次之后对多个x_eval求值时只需要反复调用newton_interp。如果你要插值的点有几百上千个这个做法的优势非常明显——差商构造是O(n²)而每个点求值是O(n)。拉格朗日插值对每个求值点都要重新计算所有基函数同样是O(n²)但乘在点的数量上数据量一大差距就拉开了。3.3 用matplotlib画出插值效果对比图实验报告里最直观的内容是一张插值效果图原始数据点、拉格朗日插值曲线、牛顿插值曲线在同一坐标系中展示。画图代码不长但有两个细节要注意插值曲线的横坐标要足够密不然曲线会显得很生硬纵坐标范围要留出余量避免曲线在边缘处超出画布。import matplotlib.pyplot as plt x_points np.array([0.0, 1.0, 2.0, 3.0, 4.0]) y_points np.array([1.0, 2.0, 1.5, 3.0, 2.8]) x_dense np.linspace(0, 4, 200) fig, ax plt.subplots(figsize(8, 5)) ax.plot(x_points, y_points, o, label原始数据点) ax.plot(x_dense, lagrange_interp(x_points, y_points, x_dense), -, label拉格朗日插值) coef newton_divided_differences(x_points, y_points) ax.plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛顿插值) ax.set_xlabel(x) ax.set_ylabel(y) ax.legend() ax.grid(True, linestyle--, alpha0.6) plt.show()画图时x_dense用np.linspace(0, 4, 200)生成200个均匀分布的点具体点数取决于插值多项式的次数和区间跨度一般200到500足够。如果两条插值曲线完全重合说明两个算法在精度上一致这是正常的如果出现细微差异优先检查差商表的迭代顺序是否正确。对比图存在的意义是验证两个算法结果一致顺带检查曲线的形状是否符合直觉。4. 最小二乘曲线拟合把公式拟合到数据点上的完整流程4.1 最小二乘原理与numpy.polyfit最小实现拟合和插值的核心区别在于插值要求曲线严格穿过每个数据点拟合只要求曲线在整体上最接近数据点。最小二乘的思想是找一组参数让模型计算值与观测值的平方误差之和最小。平方误差比绝对误差更敏感而且方便求导所以成为最常用的拟合准则。import numpy as np x_data np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]) y_data np.array([0.1, 0.8, 1.8, 2.9, 5.1, 7.2, 10.0]) degree 2 coeffs np.polyfit(x_data, y_data, degree) p np.poly1d(coeffs) print(拟合系数, coeffs) print(拟合多项式, p)np.polyfit返回的是多项式系数从最高次到常数项排列系数顺序和poly1d的输入是对应的。degree参数是最关键的调试项取1时是线性拟合取2时是抛物线拟合取大了会出现过拟合——曲线为了贴近个别数据点而剧烈抖动尤其是在数据本身不是高阶多项式的情况下。判断次数是否合适的标准很简单看拟合曲线走势是否符合数据的物理含义再看残差是否随机分布而不是呈现系统性偏差。poly1d是把系数数组包装成可以直接调用的多项式对象调用p(x_eval)就能批量求值。这里有个小细节polyfit内部用的是基于SVD的稳定算法对于次数不超过8的多项式拟合很少出现数值问题但超过10次就要警惕矩阵病态结果可能不可信。所以实验题里如果要求拟合次数特别高建议检查一下拟合系数是否出现数量级异常的大值那通常是条件数过大导致的数值性翻车。4.2 自定义函数拟合scipy.optimize.curve_fit的使用实验题里有时候不要求多项式拟合而是指定一个函数形式比如指数衰减模型或者带偏移的正弦函数。np.polyfit只支持多项式这时候要换用scipy的curve_fit它可以拟合任意你写出的Python函数。核心用法是先定义一个函数形式函数的第一个参数是自变量后面是待拟合参数。from scipy.optimize import curve_fit def exp_model(x, a, b, c): 指数衰减加常数偏移的拟合模型 return a * np.exp(-b * x) c popt, pcov curve_fit(exp_model, x_data, y_data, p0[1.0, 1.0, 0.0]) a_opt, b_opt, c_opt popt print(f拟合参数: a{a_opt:.4f}, b{b_opt:.4f}, c{c_opt:.4f}) x_smooth np.linspace(x_data.min(), x_data.max(), 100) y_fit exp_model(x_smooth, *popt)p0是参数初值这个参数在很多时候决定了拟合能不能收敛。初值差得太远算法可能陷入局部极小值甚至直接报错。给p0的经验法则是先看数据点的取值范围大致估计一下参数的量级。比如数据y值在0到10之间c的初值可以设为数据中点或最小值附近a的初值设为y的跨度b的初值根据曲线衰减的快慢估一个正数。如果你发现拟合结果明显不合理最先该改的是p0而不是模型本身。curve_fit返回的pcov是参数协方差矩阵对角线元素开根号是参数的标准差代表拟合参数的不确定度。实验报告里如果想体现严谨性可以把参数误差也打印出来。当协方差矩阵里的值大到离谱时往往说明模型参数之间存在强相关性或者数据量不足以确定这么多参数这时候的拟合系数虽然能算出数但不能当作物理解释使用。这个点算是一个隐藏的知识点写在报告的分析部分会很加分。4.3 拟合效果的量化评估残差、R方与可视化拟合做完不评估就写报告等于黑匣子输出直接交差。评估拟合质量最常用的是三个量残差、均方根误差和R方。残差是每个数据点的观测值减拟合值R方表示拟合曲线解释了数据方差的百分比越接近1越好。residuals y_data - exp_model(x_data, *popt) ss_res np.sum(residuals**2) ss_tot np.sum((y_data - np.mean(y_data))**2) r_squared 1 - ss_res / ss_tot rmse np.sqrt(np.mean(residuals**2)) print(fR方 {r_squared:.4f}) print(f均方根误差 RMSE {rmse:.4f})注意R方有一个陷阱当模型形式比较复杂、参数很多时R方会虚高甚至接近1但这不代表模型真的好。一个十次多项式拟合十个数据点R方必然是1但那完全没有预测意义。所以评估时要同时看残差图——如果残差随x增大而增大说明模型形式不对如果残差呈现波浪形说明漏掉了周期性分量。这些判断是数值分析实验最有价值的训练比单纯算对拉格朗日插值结果重要得多。把残差图画出来和拟合效果图画在一起报告的分析部分立刻有了深度。5. 实验四避坑指南5个最常见的翻车点与排查方法5.1 拉格朗日插值次数过高导致龙格现象数据越多结果越离谱现象用拉格朗日插值去逼近某些连续函数插值节点增加到十几个甚至更多插值曲线在区间两端剧烈震荡误差反而比少用几个节点时更大。原因这是数值分析教材里反复强调的龙格现象。拉格朗日插值多项式在等距节点下对龙格函数、绝对值函数这类目标做高次逼近时插值余项在端点附近会放大多项式次数越高震荡越明显。解决明确实验四的题目通常只要求4到6个节点的插值不需要追求高次。如果题目要求对比不同节点数下的误差把插值点数从3逐个增加到10画误差曲线让读者直观看到误差先降后升的过程这本身就是一份合格的现象分析。如果实际工程中遇到类似需求应该改用分段插值或样条插值而不是硬上高次多项式。5.2 图表的横纵坐标出现方块字中文标注全部乱码现象matplotlib画图时标题和图例里的中文显示成一个个小方块英文和数字正常实验报告截图里一眼就能看出问题。原因matplotlib默认字体不包含中文字形需要显式指定中文字体。Windows系统一般有SimHei或Microsoft YaHeiLinux和macOS的字体名不同用系统自带字体名逐个尝试。解决在画图代码最前面加上字体设置让matplotlib全局使用中文字体如下所示import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, PingFang SC] plt.rcParams[axes.unicode_minus] False第二行去掉坐标轴负号的特殊处理是因为切换了中文字体之后负号有可能显示成一个奇怪的方框。这个问题属于程序能跑但报告不能看的隐形翻车查错成本高提前设置好就能完全避免。5.3 差商表计算顺序写反牛顿插值结果对不上拉格朗日现象牛顿插值和拉格朗日插值结果不一致差距可能很小也可能很大越小规模的数据差距越小很多同学因此没发现。原因newton_divided_differences函数里内层循环的方向错了。差商表必须从后向前原地更新如果写成从前向后计算新的差商时用的是已经被本轮覆盖过的值导致差商全部错误。数据点少时误差不明显数据点多时完全对不上。解决内层循环固定使用for i in range(n - 1, j - 1, -1)从后往前遍历。排查思路也很简单用三个点的手工差商计算结果和Python输出对比如果差商对不上问题就在循环方向。这类错误的特点是代码能跑、结果贴近但不对最容易在实验报告的误差分析部分暴露靠直觉很难发现。5.4 拟合结果出现巨大系数多项式次数太高导致矩阵病态现象polyfit指定degree9或更高时拟合系数出现10的6次方以上的大数拟合曲线在数据点之间剧烈波动部分区间甚至超出数据范围好几倍。原因高次多项式拟合的等效设计矩阵条件数随次数指数增长数据点x跨度大时尤其严重。矩阵接近奇异SVD分解得到的解被微小舍入误差放大系数自然不稳定。解决首先判断题目是否真的要求高次拟合。如果要求做不同次数拟合的对比把次数范围控制在1到6之间重点观察拟合误差变化趋势。如果必须拟合高次多项式考虑先对x数据做归一化把数据平移到零均值并缩放到单位方差再拟合最后把系数还原回原始坐标。归一化这步在实验报告里是加分项能体现你对数值稳定性的理解。5.5 求值点落在插值节点范围之外结果完全失控现象用插值函数计算区间外的点计算结果非常大和真实函数值差出几个数量级有的同学以为是代码精度不够。原因插值多项式只在节点覆盖的区间内有可靠逼近能力。区间外没有任何数据点约束多项式形状高次项起主导作用结果必然外溢。这是插值方法的固有属性不是bug。解决求值前先检查x_eval是否在x_points的范围内。如果题目故意要求计算外推值那么在报告里要明确记录这个现象用误差对比说明外推不可靠这属于正面的实验结果。如果在做拟合预测用polyfit拟合出的多项式做短距离预测是可行的做远距离预测同样会失效。边界意识是从数值分析实验里能带走的最重要工程习惯。6. 用一份可视化脚本把实验结果整理成答辩图表实验四答辩或提交报告时图表质量往往直接影响印象分。同样是解决同一个问题一个用默认绘图风格、图例位置遮挡曲线、坐标轴没标单位和另一个图面干净、曲线粗细分明、每个子图都有标题和注释给阅卷人的感受完全不同。整理图表这一步建议独立写一个脚本不要和算法代码混在一起输出为高分辨率图片再插入报告而不是直接截图。import matplotlib.pyplot as plt import numpy as np fig, axes plt.subplots(1, 3, figsize(14, 4)) # 第一个子图插值结果对比 axes[0].plot(x_points, y_points, o, label数据点) axes[0].plot(x_dense, lagrange_interp(x_points, y_points, x_dense), label拉格朗日) coef newton_divided_differences(x_points, y_points) axes[0].plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛顿) axes[0].set_title(插值结果对比) axes[0].legend() # 第二个子图拟合效果 axes[1].plot(x_data, y_data, o, label观测数据) axes[1].plot(x_smooth, y_fit, -, label拟合曲线) axes[1].set_title(最小二乘拟合) # 第三个子图残差分布 axes[2].plot(x_data, residuals, o-, label残差) axes[2].axhline(0, colorgray, linestyle--) axes[2].set_title(拟合残差分布) axes[2].legend() fig.tight_layout() plt.savefig(experiment4_summary.png, dpi200, bbox_inchestight) plt.show()一处细节是plt.savefig里的参数dpi至少200才够打印清晰bbox_inchestight可以把空白边缘裁掉避免图片四周留白过多。如果你是准备提交电子版报告300dpi更好图片文件大一点没关系打印效果匀称。另一个细节是每个子图的标题要能自解释也就是说读者不看正文只看图就知道这张图在表达什么这是实验报告图表规范的基本要求多数同学没意识到所以认真做反而会留下印象。我自己的习惯是在画图脚本开头放一段注释记录这份图表对应的题目要求、数据文件、生成日期方便复查时快速定位。数值分析实验四做完回头看最大的收获往往不是那几个算法本身而是理解了“数据点给出的信息量是有限的任何试图从有限数据外推的做法都必须说明假设和边界”。这个意识在后面的课程和工程项目里反复被验证。希望这篇笔记能帮你把实验四从头到尾走通拿到一个能进报告、能讲清楚、也能经得起提问的结果。本文还有配套的精品资源点击获取