简介这是一份面向数值分析课程学习者的 MATLAB 上机作业资源内容对应东南大学出版教材中的典型编程题涵盖三次样条插值函数、重积分计算复化梯形公式与外推思想以及常微分方程初值问题的 RK4、AB4、AB4-AM4 及改进预测校正方法。资源包共 1 个 docx 文档体积约 566KB文档内包含三道上机题的算法说明、MATLAB 实现思路、计算结果、误差对比与编程感悟便于对照练习和复盘编程细节。该资源已有 904 人浏览学习适合正在修读数值分析、需要完成 MATLAB 上机作业或复习数值计算方法的本科生与研究生参考。读者可从中获得各题目的完整求解流程、关键公式与程序框架尤其能帮助理解不同数值方法的精度差异和实现要点。1. 数值分析上机题不是“算题”是验证算法行为的实验课某同学拿到一份上机题题目是“用Newton法求方程的一个根画出误差下降曲线”。他两分钟写完调了现成库输出一个根觉得完事了。成绩出来很低因为报告里没有误差表、没有收敛阶分析、没有初值影响的讨论。这不是个例——数值分析上机题真正考察的从来不是“能不能算出答案”而是你有没有亲眼见过算法在什么条件下失效、误差怎么涨怎么跌、哪些参数一拍脑袋就翻车。这门课的上机题通常落在插值、拟合、求根、数值积分、常微分方程初值问题这几大块适合正在修数值分析课程的学生也适合想补算法边界感的从业者。接下来我按自己交过多份作业、也帮人排查过的经验把一套能复现、能解释、能拿分的做法讲透。2. 上机题的题型地图与工具选型先看清题目在考什么拿到一份数值分析上机题先别急着敲代码。绝大多数题目表面是在“用某某方法算某某问题”实际考察的是三件事算法复现是否忠实、数值实验是否完整、结果分析是否严谨。搞清楚题目真正要什么比选工具更重要。2.1 上机题常见六大题型与真实结构数值分析上机题翻来覆去就那几个方向我把它们和考察点列成一张表题型常见载体真正的考察点非线性方程求根Newton法、弦截法、二分法初值敏感性、迭代收敛阶、导数计算是否正确插值Lagrange插值、Newton插值、样条插值Runge现象、节点分布对误差的影响拟合最小二乘多项式拟合、正交多项式拟合条件数、拟合阶数选择、过拟合数值积分梯形法、Simpson法、Romberg算法外推思想、精度阶、振荡与奇异被积函数常微分方程初值问题Euler法、RK4、Adams方法步长选择、稳定性、局部截断误差线性方程组迭代Jacobi、Gauss-Seidel、SOR收敛条件、谱半径、迭代矩阵构造一份典型上机题由三部分构成算法实现约三成、数值实验约四成、结果分析约三成。算法实现要求你不用现成库函数把方法写出来数值实验要求你对不同参数做对比比如改步长、改节点数、改初值结果分析要求你解释现象画误差图算收敛阶。理解了这三块你就明白为什么“调库算出答案”只能拿个及格分。2.2 MATLAB还是Python选型的标准不是“哪个好用”很多课程默认用MATLAB因为讲义里的示例代码是MATLAB写的老师改作业方便。但如果你有选择权我的建议是跟课程走别跟情绪走。课程指定了MATLAB就用MATLAB虽然要装商业软件但矩阵运算和绘图确实顺手课程没指定我用Python——numpy和matplotlib足够覆盖上机题的全部需求而且后续做数值实验、写报告、做自动化批处理更自由。这里有个容易踩的坑有的同学觉得“反正我用Python实现算法库里也有现成的”于是把牛顿时改成fsolve把积分改成quad。一旦题目要求“不使用内置求解函数”这种行为会被直接扣分。我自己的做法是核心算法一律自己写最多用numpy做向量计算内置函数只用来交叉验证结果对不对。这样既满足题目要求又能用已知解去核对实现有没有写错。2.3 用评分标准倒推作业结构拿到题目先列一个问题清单算法部分有没有要求不调库实验部分要求对比几个参数分析部分要求画什么图、算不算收敛阶报告有没有格式模板把评分点拆出来之后代码和文档的结构就清楚了。常见做法是先写一个“主程序”依次完成算法调用、实验循环、误差输出、画图最后用report模板把数据填进表格。如果题目明确要求“分析收敛阶”那你需要提前知道收敛阶怎么算对于求根问题观察误差比[e_{k1}/e_k^p]的极限对于积分和ODE问题通常做网格加密实验比较步长减半后误差的变化比例。后面第4章我会给一套可以直接套用的流程。3. 用 Python 跑通四个经典上机题求根、插值、积分与 ODE这一章是我最想让你照着敲的部分。四个题目覆盖最常见题型代码我都写过验证版关键参数旁边有注释跑完再对照说明调参数你会明显感觉到“原来上机题是这么设计的”。3.1 Newton-Raphson求根最小实现与收敛阶验证Newton法是最常出现在第一道题里的方法因为它短、好写、但坑也多。最小实现可以写成这样import numpy as np def newton(f, df, x0, tol1e-10, max_iter100): Newton法求根返回迭代序列和误差序列 x x0 seq [x] errors [] for k in range(max_iter): fx f(x) dfx df(x) if abs(dfx) 1e-14: # 导数接近零时迭代会爆炸 raise RuntimeError(导数过小请换初值或改用弦截法) x_new x - fx / dfx errors.append(abs(x_new - x)) x x_new seq.append(x) if errors[-1] tol: break return np.array(seq), np.array(errors) # 测试x^3 - 2x - 5 0根约为 2.0945... f lambda x: x**3 - 2*x - 5 df lambda x: 3*x**2 - 2 seq, errors newton(f, df, x02.0) print(f根 {seq[-1]:.12f}, 迭代次数 {len(seq)-1})逻辑说明每步用当前点切线近似函数公式是 x_new x - f(x)/f(x)误差序列记录相邻两步迭代值之差。这个“相邻差”在收敛时很接近真误差是判断停机的常用指标。参数里最关键的是x0和容差tolx0选在根附近且导数为零的区域迭代会慢或发散tol设太小会导致后期在机器精度附近震荡白白多迭代几十次。我一般先画函数曲线看根大概在哪再选x0这样比盲猜稳得多。收敛阶验证是这题的加分项。理论上Newton法对单根是二阶收敛也就是说误差每次大约平方递减。你可以这样验证把errors相邻项做除法观察比值是否趋近某个常数for k in range(1, len(errors)): if errors[k-1] 0: print(k, errors[k] / errors[k-1]**2)如果最后一组比值稳定在一个数量级说明程序实现了二阶收敛。如果比值越来越大多半是求导公式写错或者初值离根太远导致一开始没进入局部收敛区域。这个自查方法也适用于其它迭代法只是p值不同比如弦截法是1.618阶。3.2 Newton插值复现Runge现象才算真正看懂插值插值类上机题几乎必考Runge现象对区间上的连续函数做高次等距节点插值次数越高边界附近振荡越剧烈。这个现象不亲手画出来光看公式很难有体感。下面这段代码实现了Newton插值并在Runge函数上做对比import numpy as np def newton_interp(x_nodes, y_nodes, xx): Newton差商插值xx是待求点数组 n len(x_nodes) coef y_nodes.copy().astype(float) # 计算差商表 for j in range(1, n): for i in range(n-1, j-1, -1): coef[i] (coef[i] - coef[i-1]) / (x_nodes[i] - x_nodes[i-j]) # 用嵌套形式求值 res np.zeros_like(xx, dtypefloat) for k, xk in enumerate(xx): p coef[-1] for i in range(n-2, -1, -1): p p * (xk - x_nodes[i]) coef[i] res[k] p return res f lambda x: 1 / (1 25*x**2) # 经典Runge函数 a, b -1, 1 for n in [5, 10, 15]: nodes np.linspace(a, b, n1) # 等距节点 xx np.linspace(a, b, 200) yy_truth f(xx) yy_interp newton_interp(nodes, f(nodes), xx) err np.max(np.abs(yy_interp - yy_truth)) print(fn{n:2d}, 最大误差{err:.3e})逻辑说明差商表逐列计算嵌套求值是数值上比较稳的评估方式避免展开成单项式导致大数相消。等距节点下n5还算正常n10误差已经变大n15边界处会出现明显振荡。这就是Runge现象的实测版。判断插值质量不能只看在节点处对不对——节点上肯定对要看节点之间的最大偏差所以我在xx上取了200个密集点做误差统计。如果题目要求对比“等距节点”和“Chebyshev节点”把节点生成改成余弦节点nodes np.cos(np.pi * np.arange(n1) / n) # Chebyshev节点两端密中间疏其余代码不动。你会看到误差数量级立刻掉下来几个量级。这说明一个很重要的问题算法本身没变只是节点分布变了效果天差地别。上机题让你画这个对比目的就是别把“插值算法差”和“节点选得差”混为一谈。3.3 Romberg积分外推思想的代码表达数值积分上机题里Romberg算法是考察重点因为它把梯形法的序列通过外推公式变成高精度结果。先看实现def romberg(f, a, b, m10): Romberg积分m为外推行数返回最终近似值和递推表 R np.zeros((m, m)) h b - a R[0, 0] (f(a) f(b)) * h / 2 # 初始梯形值 for k in range(1, m): h h / 2 # 补算新增的中间点 s 0.0 for i in range(1, 2**k, 2): s f(a i * h) R[k, 0] 0.5 * R[k-1, 0] h * s # 梯形递推 # Richardson外推 for j in range(1, k1): R[k, j] R[k, j-1] (R[k, j-1] - R[k-1, j-1]) / (4**j - 1) return R[m-1, m-1], R # 测试∫0^1 x^3 dx 0.25 f lambda x: x**3 val, table romberg(f, 0, 1, m8) print(fRomberg结果 {val:.12f}, 误差 {abs(val - 0.25):.2e})逻辑说明R的第一列是逐步减半步长的梯形值每多一行步长减半一次第二列起把前一行的两个值组合消去误差主项。公式里的4**j - 1来自梯形法的误差按h^2展开每外推一次精度提高两阶。这段代码的要点是“按行推进而不是按列推进”——必须先算完第k行的梯形值才能做第k行的外推。参数m代表外推行数不是步长数。m8时实际上最大的梯形划分是2^7个区间。对光滑函数m10就能达到机器精度附近的水平。但要注意如果被积函数在端点附近有奇异性Romberg收敛会很慢甚至看起来不收敛这时候要先做变量替换或分解奇异项而不是一味增大m。这道题的观察点在于m从3增加到10误差的下降速度在变快这就是“外推加速”的实际效果。3.4 RK4求解初值问题步长是个要动手试的参数ODE初值问题是上机题里最容易“结果看着对、实际全错”的题目。RK4是经典方法实现不长def rk4(f, t0, y0, t_end, h): 经典四阶Runge-Kutta求解 y f(t, y) t, y t0, y0 ts, ys [t0], [y0] while t t_end - 1e-12: k1 f(t, y) k2 f(t h/2, y h/2 * k1) k3 f(t h/2, y h/2 * k2) k4 f(t h, y h * k3) y y h/6 * (k1 2*k2 2*k3 k4) t t h ts.append(t) ys.append(y) return np.array(ts), np.array(ys) # 测试y y, y(0)1, 真解 e^t f lambda t, y: y ts, ys rk4(f, 0.0, 1.0, 2.0, h0.1) print(ft2处近似值 {ys[-1]:.10f}, 真值 {np.e**2:.10f})逻辑说明k1到k4是区间内四个不同位置的斜率估计最后按1:2:2:1加权平均。这个“加权”不是拍脑袋而是保证局部截断误差为O(h^5)的组合方式。步长h决定精度和总步数h减半误差理论上降为原来的1/16但总步数翻倍舍入误差的累积也会增加。所以不是h越小越好——实际中你会看到误差随h先降后升这是正常的找那个“误差最低点”附近的h就是步长的合理取值。这道题的实验部分可以这样设计分别取h 0.2, 0.1, 0.05, 0.025在固定点t2处计算误差然后打印误差比。如果程序没错误差比应该接近16。如果比值只有10或者2说明代码实现有逻辑问题或者被积函数不够光滑导致低阶分量占主导。这个“四阶”验证本身就是上机题里最值钱的输出。4. 上机题实验设计制造已知解、收敛阶验证与误差表代码写完只是第一步上机题的大部分分值在实验设计和结果分析。很多同学在这个环节翻车不是算法错而是不知道怎么证明“我的实现是对的”。这里有一套通用实验设计流程可以直接套在任何题目上。4.1 制造已知解让算法用“标准答案”自检对积分和ODE问题最可靠的自检方法是构造一个解析解已知的测试问题。比如你想测ODE的RK4实现选一个yλy这样有精确解e^{λt}的方程就可以把数值解和精确解逐点对比。积分问题同样直接选∫_0^1 e^x dx真值是e-1一减就知道程序写没写对。这个方法常被称为“制造已知解”严格形式化之后也用在偏微分方程数值解的误差估计上。操作步骤是先选简单函数多项式、指数、三角都可以算好解析解再跑你的算法最后算最大误差。如果最大误差在预期精度范围附近说明实现基本可信如果误差大得离谱先查代码逻辑别急着改参数。固定一个测试用例还有个额外好处——后面调任何参数都能拿这个用例做回归检查改坏了立刻能发现。以RK4为例可以这样制造一组不同频率的测试问题for lam in [1, -1, 10, -10]: f lambda t, y: lam * y ts, ys rk4(f, 0.0, 1.0, 1.0, h0.05) true np.exp(lam * ts) err np.max(np.abs(ys - true)) print(flambda{lam:3d}, 最大误差{err:.3e})注意lambda-10时步长0.05已经不满足稳定性条件——这正是实验的目的之一。负特征值大而步长不够小数值解会发散这是ODE方法的稳定性边界亲眼看到它发生比背两页判据有用得多。4.2 收敛阶验证与loglog图斜率会告诉你代码对不对收敛阶验证是报告里最有说服力的部分。做法是把步长或节点数按倍数变化记录误差变化然后画双对数图。斜率的理论值就是方法的收敛阶。比如RK4在loglog图中斜率应为4梯形法为2Newton法如果画“误差-迭代次数”单对数图斜率会呈指数下降趋势。画图代码这样写import matplotlib.pyplot as plt hs [] errs [] for h in [0.2, 0.1, 0.05, 0.025, 0.0125]: ts, ys rk4(f, 0.0, 1.0, 2.0, h) err abs(ys[-1] - np.e**2) hs.append(h) errs.append(err) print(fh{h:.4f}, err{err:.3e}) plt.loglog(hs, errs, o-) plt.xlabel(步长 h) plt.ylabel(误差) plt.title(RK4 误差随步长变化) plt.grid(True, whichboth, ls--, alpha0.3) plt.show()如果你看到的数据点接近一条直线这条直线斜率大约是4可以放心写“实现了四阶精度”。如果斜率只有1或2大概率是代码里某个k项写错了或者输出数据时选错了观测点。有个常见误用很多人把相邻步进之间的差当误差但那个在收敛阶验证里不能替代真误差——只有和精确解比出来的误差才反映精度阶。4.3 误差表设计真误差、误差比与收敛阶怎么填报告里的误差表是评分老师第一眼扫的地方。设计要遵循一个原则让读者一眼看出误差在下降、下降速度符合预期。列的内容我建议固定为四列步长h、最大误差、误差比、推测收敛阶。h最大误差误差比推测阶0.22.6e-3--0.11.7e-415.3~3.90.051.1e-515.5~4.00.0256.8e-716.2~4.0误差比的计算方法是后一行的误差除以相邻前一行的误差。当h减半时如果方法是p阶误差比应接近2^p。看到误差比稳定在16附近就可以判断RK4代码写对了。注意第一行没有误差比这是正常的不用填“-”以外的内容。误差比小幅波动是正常的不必过分紧张。但如果误差比在某一档突然变小或变大很多先检查是不是误差受舍入误差控制了——当误差低于1e-12量级机器精度就会压过截断误差曲线会变平甚至抖动这在双精度浮点数下是物理规律不是bug。报告里如果能主动写一句“当h进一步缩小时误差受机器精度限制不再下降”反而是加分的分析。5. 上机题避坑六个把同学搞到半夜的细节这些坑我基本都踩过也帮人排查过。每一条都按“现象→原因→解决”写你遇到类似问题时可以直接对照。5.1 现象Newton法迭代几次就输出NaN原因初值选在了导数接近零的区域一步迭代把x推到很远或者导数公式写错符号反了。解决先把函数图像画出来看根在哪再选初值在代码里加导数绝对值检查碰到|f(x)|1e-14就主动报错换初值。还有一个常见偷懒做法是手算导数容易算错可以用差分近似求导但注意这会把收敛阶降到一阶题目要求“精确导数”时不适用。5.2 现象插值多项式次数越高边界误差反而越大原因这就是Runge现象等距节点下高次插值在边界处振荡次数越高振荡越厉害误差会随次数增长。解决换Chebyshev节点或者改用分段低次插值。两位同学的代码对比最能说清楚问题一个用等距节点n15画出的图在两端像波浪另一个用Chebyshev节点画出的图贴着真函数两张图放在一起谁理解了插值的边界条件一目了然。上机题如果只要求“比较两种节点”这个对比就是核心输出别只贴代码不贴图。5.3 现象数值积分“步长减半误差反而变大”原因被积函数可能在积分区间内有奇点或者函数振荡太快梯形公式和Simpson公式的误差估计前提不满足。另一个可能是误差已经降到机器精度附近后面的波动全是舍入噪声。解决先画出被积函数曲线自查有没有尖峰或振荡积分区间是否包含奇点——有的话要做奇异点分离或变量替换。误差低于1e-13时直接说明“已到机器精度”不再继续加密网格把这点写进报告反而是亮点。5.4 现象库函数和你手写算法的结果不一致原因内置函数用自适应算法容差设置不同比如quad的默认绝对容差是1e-8你手写的简单方法精度可能更高也可能更低这不能说明谁对谁错。解决用已知解析解做裁判。统一把问题换成有精确解的测试函数两边分别和真值比误差小的一方在“这一题、这一容差”下更接近准确答案。库函数不是检验标准真解才是。5.5 现象误差表看起来正常但收敛阶图一片模糊原因图的横轴步长用了线性刻度而误差跨了几个量级压缩到一个区域内几乎看不见趋势或者对收敛阶的纵轴用了线性坐标。解决用loglog图横纵轴都取对数。这是一个让收敛趋势立刻显现的简单操作。还要注意数据点的单位必须一致h是步长就用步长别忽然后面几组用了区间数——这在生成循环时容易写错。5.6 现象同样的代码在自己的机器上一个结果换台机器换个结果原因版本差异、浮点环境差异或者代码里依赖了未初始化的变量。解决报告里固定写明软件版本和运行环境例如“Python 3.10 with numpy 1.24”。数值实验要求可复现性一个没有版本信息的数值结果老师很难判断是你写的还是抄的。养成在代码开头打印版本信息的习惯对你自己排查也有帮助。6. 让上机题变成可复现实验的三个进阶习惯如果你想把上机题从“交作业”提升到“真正可控的数值实验”我建议养三个习惯全部在平时就能做到。第一个习惯是所有实验循环里用同一套测试问题固定解。不要在每一道题里临时想测试函数而是维护一个小函数库里面放三到五个有解析解的测试问题多项式积分、指数函数ODE、带振荡的Runge函数。每次算法改动都用这些测试跑一遍回归保证旧的正确性没有被破坏。这个习惯在很多工程项目里同样适用数值作业只是最早练手的地方。第二个习惯是导出数据而不是复制粘贴。写一个简单的辅助函数把误差表按Markdown或LaTeX表格格式输出到文件然后报告直接插入。手动抄数据容易抄错而且改一次参数就要重新抄一遍。代码量不大但省下的时间和避免的错误非常可观。我一般会在循环里直接构造二维数组然后一次性打印成表格。第三个习惯是记录每次实验的参数配置。步长、节点数、初值、容差、最大迭代次数这些值哪怕只改一个结果就不同。我在代码文件头部写一行注释记录当前实验用的参数组合跑完一次实验就复制一份结果图和数据文件命名里带上关键参数比如rk4_h0p05.png。这样回头写报告时不会面对一堆同名文件发懵。这三个习惯一开始会显得繁琐但等你同时做三四个上机题时它们能省下大量对账时间。最后的教训是数值分析上机题最坑人的地方从来不是数学公式而是你默认“程序跑出了数就算完成任务”。把每一步实验当成一次有计划的计算——先想清楚要验证什么再写代码最后用数据说话。整个过程里我最大的收获不是那些分数而是学会了对计算结果的怀疑态度不对答案不画趋势不算阶数就不敢说这个结果是可信的。希望帮到你。本文还有配套的精品资源点击获取