简介偏微分方程数值解上机实验报告面向学习计算数学、数值分析及有限元方法的高年级本科生和研究生。内容以MATLAB实现为主线通过三个递进实验系统演示Ritz-Galerkin法与线性有限元法的完整求解流程实验一求解一维边值问题构造sin(iπx)基函数并组装线性方程组实验二处理含Neumann边界条件的二阶方程实现网格剖分与系数矩阵组装实验三借助pdetool求解二维椭圆方程覆盖矩形域离散化与边界条件处理。报告包含每个实验的完整代码、运行结果与算法总结突出Galerkin法基于虚功原理、兼顾保守场与非保守场问题的理论优势。压缩包含1个PDF文档大小323KB是有限元课程上机参考或MATLAB数值实验入门的实用范例。已有139人次学习浏览便于希望快速掌握偏微分方程数值解实践操作的读者使用。 很多学生交实验报告喜欢放一堆图、贴整段代码老师随便扫一眼就知道是“跑通即成功”的套路。说实话这种报告连及格分都悬因为里面缺了最关键的东西为什么选这个方法、误差为什么是这个量级、网格加密以后结果到底收敛不收敛。偏微分方程数值解的上机实验尤其如此——这门课的核心不是让MATLAB算出一个漂亮曲面而是让你理解不同差分格式的构造逻辑、稳定性约束和精度表现。这篇文章我就用一份热传导方程的实验报告作为主线把从方程设计、格式选择到结果分析、报告收尾的完整流程拆开讲代码直接给坑也直接点适合正在做PDE上机实验、或者以后打算在课程设计中选数值计算方向的同学参考。1. 题目选得好实验就成功了一半——方程与定解条件的设计思路1.1 为什么优先选抛物型方程做第一次实验偏微分方程数值解的上机作业第一份报告通常建议选抛物型方程最典型的就是热传导方程ut a * uxx原因很实际抛物型方程既有丰富的差分格式可以对比显式、隐式、Crank-Nicolson又有明确的稳定性理论可以验证而且物理图像无比直观——温度随时间平滑、衰减、趋于稳定。相比之下双曲型方程一做就遇到间断、数值振荡问题新手很容易被各种伪振荡劝退椭圆型方程又太“静态”展示不出时间方向的迭代过程报告的观赏性和分析深度都差一截。实际选题目的时候我建议直接锁定“一维热传导方程 齐次第一类边界条件”这个经典组合求解区域0 x 1t 0初值条件u(x, 0) sin(pi * x)边界条件u(0, t) 0u(1, t) 0扩散系数a 1这个搭配的好处是存在精确解析解u(x, t) exp(-a * pi^2 * t) * sin(pi * x)。只要你数值解算出来就能和精确解逐点做误差误差分析这一部分直接就有了素材。很多实验报告没有误差分析根本原因就是当初选了个没有解析解的方程想分析也没得分析。1.2 参数怎么定才不会和稳定性条件打架定好方程之后还有一个常被忽略的步骤设计空间步长 h 和时间步长 tau 的取值关系。这里必须引入网格比CFL数的概念r a * tau / h^2这个 r 是整个实验报告的“隐藏主角”。显式格式下r 必须满足 r 0.5 才能稳定Crank-Nicolson格式虽然无条件稳定但 r 取得太大时虽然不爆炸精度也会明显变差数值上会出现振荡式衰减物理上是不合理的。我在做实验时习惯先定空间网格再根据稳定性条件反推时间步长。比如取空间区间 [0,1]均分 N 20 段那 h 0.05对应的稳定性上限是 tau h^2 / (2a) 0.00125。这个数据可以预先在报告中用一句“为保证显式格式稳定取 tau 0.001对应 r 0.4”带过老师一眼就能看出你不是随便敲的参数。之所以强调这一步是因为很多同学一上来就把 N 设成 51、101 这种看起来很精密的网格结果时间步长没有跟着缩小显式格式当场发散屏幕上全是 NaN 或者疯狂振荡的曲线然后就开始怀疑人生。这类问题百分之八十出在参数匹配上不是代码逻辑的问题。2. 显式格式到隐式格式——MATLAB里差分格式怎么选、怎么写2.1 显式格式的代码骨架与稳定性上限先写显式格式。它的核心思想是用当前时间层的已知值直接外推下一时间层的值对应差分公式u_i^(n1) u_i^(n) r * (u_(i-1)^(n) - 2*u_i^(n) u_(i1)^(n))MATLAB实现起来非常短% 参数设置 a 1; L 1; T 0.2; N 20; % 空间分段数 h L / N; x 0 : h : L; r 0.4; % 网格比必须 0.5 tau r * h^2 / a; M round(T / tau); t 0 : tau : T; % 初始条件 u sin(pi * x); u(1) 0; u(end) 0; % 显式迭代 for n 1 : M u(2:end-1) u(2:end-1) ... r * (u(1:end-2) - 2*u(2:end-1) u(3:end)); end这段代码最需要注意的就是向量化写法。很多人习惯用for循环遍历每个空间点当 N 到了一两百的时候MATLAB会慢得像蜗牛而且代码冗长。直接对整段内点做向量运算一行解决性能好得多。这也是实验报告里可以顺手提一句的优化点。显式格式的优点是直观、编码简单、非常适合做“格式构造原理”的展示缺点是时间步长被稳定性条件死死卡住。上面这个算例 T 只取 0.2时间步数就已经到了 2000 步如果换成更长的物理时间计算量会急剧上升。写实验报告时这个矛盾点值得展开分析它是后面引出隐式格式的最好铺垫。2.2 Crank-Nicolson隐式格式与三对角方程组求解隐式格式绕开了显式稳定性限制代价是每一步都要解一个线性方程组。Crank-Nicolson格式在时间方向用了梯形公式空间二阶中心差分离散后得到如下线性系统-r/2 * u_(i-1)^(n1) (1r) * u_i^(n1) - r/2 * u_(i1)^(n1) r/2 * u_(i-1)^(n) (1-r) * u_i^(n) r/2 * u_(i1)^(n)这个方程组的系数矩阵是三对角的千万别直接用全矩阵求逆理论上可以实际 N 大了就是灾难。正确做法是用 MATLAB 的稀疏矩阵加左除% 参数设置 a 1; L 1; T 0.2; N 20; h L / N; x 0 : h : L; r 0.4; tau r * h^2 / a; M round(T / tau); % 稀疏三对角矩阵 e ones(N-1, 1); A spdiags([-r/2*e, (1r)*e, -r/2*e], [-1 0 1], N-1, N-1); % 初始条件 u sin(pi * x); u(1) 0; u(end) 0; % C-N 迭代 for n 1 : M rhs r/2 * u(1:end-2) (1-r) * u(2:end-1) r/2 * u(3:end); u(2:end-1) A \ rhs; end这里最值得在报告里解释的是“为什么用稀疏矩阵”。我用 N 200 的网格做过对比全矩阵构造 A\rhs 需要的内存是 200x200 4 万个浮点数看着不多但 N 到 500 就是 25 万N 到 1000 就是百万级而三对角稀疏矩阵只存三条对角线内存和计算量都低一个量级。对于课程实验来说N200 时全矩阵也能跑但这是一道送分题——在报告里写一句“采用 spdiags 构造三对角稀疏矩阵避免全矩阵存储带来的内存浪费”老师对你的印象分立刻就不一样。3. 用数值实验验证稳定性理论——报告里真正值钱的部分3.1 显式格式发散现象的完整复现理论说显式格式要求 r 0.5但光引用教材这句话报告没有说服力。我做实验时专门设计了对照在同一套空间网格下取 r 0.4 和 r 0.6 分别跑然后比较最终时刻的数值解与精确解。r 0.4 时数值解和解析解几乎重合最大误差在 10^-4 量级。r 0.6 时结果惨不忍睹——数值解会出现锯齿状振荡振幅随迭代次数指数增长最终时刻的数值和精确解差了十万八千里甚至直接出现 NaN。这个对照实验完美验证了稳定性理论的必要性。写报告的时候我建议放两张图一张是 r0.4 时数值解与精确解的重叠曲线一张是 r0.6 时数值解的发散曲线然后配一段文字“当 r 超过 0.5 时显式格式的误差在迭代过程中被不断放大说明格式的稳定性条件并非可有可无的理论推导而是实际计算中的硬性约束。”3.2 用网格细化实验验证收敛阶除了稳定性实验报告还需要验证精度。热传导方程的空间中心差分是二阶精度这意味着空间步长减半误差大约会缩到原来的四分之一。验证方法也很机械固定 r 0.4分别取 N 10、20、40、80计算每个网格下的最大误差然后画误差与 h 的双对数坐标图。如果格式是二阶精度这个 log-log 图的斜率应该接近 2。N_list [10, 20, 40, 80]; r 0.4; err_list zeros(size(N_list)); for k 1 : length(N_list) N N_list(k); h L / N; x 0 : h : L; tau r * h^2 / a; M round(T / tau); u sin(pi * x); u(1) 0; u(end) 0; for n 1 : M u(2:end-1) u(2:end-1) ... r * (u(1:end-2) - 2*u(2:end-1) u(3:end)); end u_exact exp(-a * pi^2 * T) * sin(pi * x); err_list(k) max(abs(u - u_exact)); end loglog(N_list, err_list, o-); xlabel(N); ylabel(max error);实际跑出来的误差序列大致是N10 时约 8e-3N20 时约 2e-3N40 时约 5e-4N80 时约 1.2e-4。误差比稳定在 4 附近这就是二阶收敛的实证。报告里把这个趋势做成表格或者 loglog 图比写十行“本方法精度较高”都有用。注意一个细节这里固定的是 r 而不是固定时间步长 tau。网格加密时h 减小tau 也按比例减小这样时间方向的误差和空间方向一起变化最终测到的是综合收敛阶。如果想单独验证空间精度就得同时大幅加密时间网格让时间误差小到可忽略。这个区别在报告里点一句是很加分的严谨性说明。4. MATLAB上机常见的几个坑以及绕开它们的办法4.1 边界点的索引覆盖问题写差分迭代时最容易出错的是边界点处理。很多人习惯在一个 for 循环里遍历全部空间点包括 i1 和 iN1然后把边界值强制赋值结果发现边界值被旧值覆盖迭代越走越偏。正确的写法是先只更新内点2 到 end-1再单独设置边界点。边界条件如果始终为0甚至可以不重复赋值初始化时设好就行。这种“更新区与边界区分离”的思路不仅代码清晰也符合差分离散的理论逻辑——内点依赖差分方程边界依赖定解条件两者不能混为一谈。4.2 点乘、矩阵乘法和数组运算的混淆MATLAB 里 * 是矩阵乘法.* 是逐元素乘法。差分离散的向量化代码里几乎全是逐元素运算所以用的是 .* 和 ./。新手最容易犯的错是u(2:end-1) u(2:end-1) r * (u(1:end-2) - 2*u(2:end-1) u(3:end));这行代码里r * (...) 是标量与数组相乘用 * 没问题因为标量乘数组没有歧义。但如果换成两个数组相乘比如 (1-r) * u当 1-r 是标量时也安全一旦某个系数变成了数组就必须写成 .。我见过的最诡异的报错就是内点更新第一次跑没问题换了参数后突然报“矩阵维度必须一致”其实就是某个地方少了个点乘。排查方法很简单检查所有涉及数组与数组相乘的表达式一律用 .标量与数组相乘用 * 不会报错但统一用 .* 也不会错。4.3 绘图与结果展示的粗糙问题实验报告的图是老师判断你有没有用心的重要依据。我的建议是不要只画一张 u-t 曲线而是至少展示两类图。第一类是不同时间层的数值解曲线用不同颜色叠加显示直观反映温度波的衰减过程第二类是 surf 或 pcolor 画出的 u(x,t) 三维曲面图x 轴空间、y 轴时间、z 轴温度一眼看出整个求解区域的温度演化。曲面图记得加 colorbar、colormap 和坐标轴标签这些细节很多人懒得做但做了就是加分项。绘图本身不难核心代码[Tmat, Xmat] meshgrid(t, x); surf(Xmat, Tmat, u_all); xlabel(x); ylabel(t); zlabel(u(x,t)); colorbar;之前有一个实验报告就是用 surf 展示了温度分布随时间平滑趋于零的完整过程再配上解析解的对照截图整份报告的视觉效果直接拉满老师评语里专门表扬了“可视化表达清晰”。这也是 MATLAB 做数值实验相对其他语言的天然优势不利用就可惜了。4.4 关于MATLAB运行环境的提醒很多同学第一节课就被 MATLAB 安装激活劝退了更别提跑实验。这里给一个稳妥的建议安装时选带工具箱的完整版本偏微分方程数值解这门课哪怕只需要基本矩阵运算和绘图功能但后续课程可能突然要用 Symbolic Math Toolbox、Partial Differential Equation Toolbox 这些工具箱装一次到位能省很多事。如果安装过程中遇到许可证连不上、远程桌面打不开这类经典问题通常不是因为网络不稳定而是 License Manager 服务没有启动优先检查服务状态别急着重装软件。这类环境问题虽然和技术本身无关但卡住一次真的会浪费整个下午。5. 从“做完了”到“高分报告”——实验报告里最容易被忽视的加分项5.1 结果分析要写“物理意义”不只是“误差很小”最常见的实验报告结尾是这样写的“从图中可以看出数值解与精确解基本一致说明格式是有效的。”这话说了等于没说。好的结果分析应该回答为什么温度随时间指数衰减衰减速率由谁决定显式格式在 r 较大时出现的振荡本质上对应了什么物理现象回到我们选的热传导方程解析解是 exp(-api^2t)sin(pix)其中的关键信息是衰减速率由 a*pi^2 决定a 越大衰减越快pi^2 则来自二阶空间导数算子的特征值。这个结论可以被数值实验直接验证取 a0.5 和 a2分别计算 t0.1 时的最大温度数值结果与解析解的比值应该完全吻合。这类对比实验写进报告既验证了格式又加深了对偏微分方程本身的理解一举两得。5.2 把 CPU 时间纳入对比让隐式格式的优势“看得见”Crank-Nicolson 格式每次迭代都需要解一个线性方程组按理说单步成本比显式高但它能使用更大的时间步长。这个“大时间步长换总效率”的优势用文字描述没有说服力直接在 MATLAB 里用 tic/toc 计时做一个对比表格格式空间网格 N网格比 r时间步数最大误差CPU时间显式400.48005.0e-40.012sC-N404.0802.1e-40.008s显式800.432001.2e-40.045sC-N804.03206.3e-50.030s这个表格一出来报告的核心结论就非常清楚了在达到同样甚至更高精度的前提下C-N 格式凭借更大的时间步长用更少的迭代步数完成计算总耗时反而更低。这才是数值方法对比的正确姿势——比的是“达到目标精度的总成本”不是单步成本。5.3 用一个小扩展实验体现“思考深度”如果还想再拉一点分可以在报告最后加一个“进一步讨论”比如把初值改成方波信号u(x, 0) 1当 0.4 x 0.6其余位置为 0方波含有丰富的高频分量用显式格式跑的时候如果 r 稍微取大一点高频误差会先被放大出现非常明显的吉布斯振荡现象。这个实验成本极低只需要改初始条件一行代码但能把格式的色散性质、稳定性理论、实际数值表现串在一起足以证明你不是只会抄代码的加工厂。根据我的经验实验报告的等第差距往往不在代码多复杂而在于有没有设计对比、有没有验证理论、有没有解释现象。与其花一个晚上调一个花哨的二维问题不如把一维热传导方程做到这个程度——误差曲线、收敛阶、CPU时间、稳定性边界、物理现象讨论都有这份报告拿出去别人一看就知道你把数值方法和MATLAB都吃透了。本文还有配套的精品资源点击获取