简介这是一份面向数学、物理及工程领域学习者的偏微分方程数值解入门资料内容从椭圆型、抛物型、双曲型三类典型方程出发系统讲解泊松方程、热传导方程和波动方程的定解问题并重点展示如何在MATLAB中通过有限差分法、pdepe等工具进行数值求解适合需要结合代码理解偏微分方程数值方法的本科高年级学生或科研人员。压缩包内包含1个PDF文件体积仅1.67MB便于下载后随时查阅。目前已有222人学习浏览内容结构清晰完整覆盖网格剖分、差分格式构造、边值问题处理、解的收敛性与误差分析等关键环节同时给出差分方法的实现思路与MATLAB函数使用要点对初值和边界条件的分类也做了细致梳理。整体而言这是一份兼具理论梳理和实操参考价值的紧凑型学习手册。 手边这本从网上下载的《偏微分方程-matlab.pdf》封面已经有点旧了但我还是愿意把它放在“常翻”文件夹里。原因很简单偏微分方程的数值求解MATLAB 仍然是我见过最省心的环境。不管你是做课程作业、毕业论文还是工程项目里要快速验证某个扩散、波动或传热过程用 MATLAB 搭数值模型基本都能在很短时间里得到能直接分析的结果。这篇文章不打算堆数学公式而是把我在使用这份 PDF 资料和实际项目过程中积累的求解方法、代码模板、以及踩过的坑整理出来给同样被偏微分方程折磨过的人一份可以直接参考的实操记录。1. 拿到《偏微分方程-matlab.pdf》后我重新梳理的求解路线很多人把偏微分方程当成一门纯数学课来学动辄翻公式翻到怀疑人生。但如果你手里有 MATLAB其实可以换一种思路先把它当成工具跑通一个结果再回头补理论。我最早看这本 PDF 的时候就是被里面“用数值方法求解 PDE本质上是把无限维问题变成有限维问题”这句话点醒的。从那以后我遇到实际的偏微分方程问题基本都按同一条流程走先判断方程类型再确认定解条件然后选求解路径最后写代码、调精度、验证结果。1.1 为什么 MATLAB 值得用来啃偏微分方程MATLAB 在这类问题上有三个天然优势。第一矩阵运算是它的看家本领而有限差分、有限元、谱方法最终都落在矩阵运算上第二内置可视化工具非常成熟求解完直接就能画温度场、波场、浓度场不需要额外装库第三工具箱和内置函数覆盖了大量常用场景比如 pdepe、PDE Toolbox以及更底层的 ode 求解器用来处理半离散化的 PDE 问题也足够顺手。我在项目里经常要快速验证一个传热模型在 Python 里要用 numpy、scipy、matplotlib 搭一堆环境在 MATLAB 里只需要一个脚本文件就能从定义方程到出图全部完成。当然Python 也很优秀但如果你手上的资料就是《偏微分方程-matlab.pdf》这种教材按它的思路走是最省力的。1.2 从连续到离散贯穿所有 PDE 求解的一条主线偏微分方程描述的是连续变量之间的关系但计算机只能处理有限个数字。所以数值求解的核心思路就是把连续问题离散成代数方程组。比如一维热传导方程[ \frac{\partial u}{\partial t} \kappa \frac{\partial^2 u}{\partial x^2} ]计算机并不知道这个偏导数怎么算我们得把空间 (x) 和时间 (t) 切成很多小段在网格点上用差分公式近似导数。这个过程叫离散化。离散化之后原来的 PDE 就变成了一系列关于网格点上的值 (u_{i}^{n}) 的代数关系式。MATLAB 做的无非就是帮你高效地解这些代数关系式。理解这条主线比记住某一类方程的特定解法重要得多。因为不管是 pdepe、PDE Toolbox还是自己写有限差分本质上都在做这件事。1.3 拿到方程先做类型判断再决定用哪个求解器不是所有偏微分方程都能用同一个函数解。拿到一个 PDE我一般先问三个问题方程是什么类型椭圆、抛物还是双曲自变量是哪些维度有多高边界条件和初始条件给的是完整这三个问题直接决定选型。我把常用场景整理成了下面的表格方便对照方程类型经典例子数学特征MATLAB 常用路径椭圆方程拉普拉斯方程、泊松方程稳态问题没有时间项PDE Toolbox或自写迭代法抛物方程热传导方程、扩散方程含一阶时间导数扩散特性pdepe或自写隐式差分双曲方程波动方程、对流方程含二阶时间导数或一阶时间导数强对流项pdepe 部分适用复杂情况自写差分/FVM我见的比较多的坑是有人拿椭圆方程的求解器去解双曲方程结果波峰处产生剧烈振荡还以为是参数设置问题。其实第一步方程类型没判断对后面全是白费。判断类型的方法不复杂就看最高阶导数全是空间二阶导、没有时间项的基本是椭圆含时间一阶导和空间二阶导的基本是抛物含时间二阶导或者有强一阶对流项的基本是双曲。2. 动手前必须想清楚的事定解条件、离散化与边界处理我见过不少同学方程拿到手就急着写代码结果算出来一张漂亮的图拿去和解析解一比完全对不上。原因往往不是差分格式不对而是边界条件或者初始条件没理清楚。定解条件就像题目的“额外规则”少了它方程的解不唯一错了它整个数值解就是一场大型幻觉。2.1 三类边界条件必须熟到条件反射偏微分方程的边界条件常见有三种MATLAB 里 pdepe 的写法也基于这三类第一类边界条件Dirichlet直接给定边界上的函数值。比如杆的一端温度固定为 100 度就是 (u(0,t)100)。第二类边界条件Neumann给定边界上的导数或通量。比如杆的末端绝热即传热通量为 0也就是 (\frac{\partial u}{\partial x}0)。第三类边界条件Robin函数值和导数的线性组合常见于对流换热。比如边界上有热交换可写成 (\alpha u \beta \frac{\partial u}{\partial x} \gamma)。在 pdepe 里面边界条件统一写成 (p q \cdot f 0) 的形式其中 (p) 和 (q) 是我们需要返回的函数值(f) 是 PDE 方程里的通量项。一开始接触会有点绕但多写几次就顺了。2.2 离散化基础前向、后向、中心差分离散化是所有数值解法的地基。以一维网格点 (x_i) 上某个函数 (u(x_i)) 为例常见的差分近似有前向差分(\frac{du}{dx} \approx \frac{u_{i1} - u_i}{\Delta x})后向差分(\frac{du}{dx} \approx \frac{u_i - u_{i-1}}{\Delta x})中心差分(\frac{du}{dx} \approx \frac{u_{i1} - u_{i-1}}{2\Delta x})二阶导 (\frac{d^2u}{dx^2} \approx \frac{u_{i1} - 2u_i u_{i-1}}{\Delta x^2})其中中心差分的精度通常更高而前向或后向在显式时间推进时实现更简单。我个人习惯空间二阶导一律用中心差分时间一阶导根据情况选择显式欧拉或隐式格式。2.3 在 MATLAB 里“定义微分方程”的标准姿势不少人搜过“matlab 中定义微分方程”这类词其实在 MATLAB 里定义微分方程的关键不是像写数学公式一样写出一整行 (dt/dx)而是把方程拆成标准形式再写成函数。以 pdepe 为例它要求方程写成[ c(x,t,u,\frac{\partial u}{\partial x}) \cdot \frac{\partial u}{\partial t} \frac{\partial}{\partial x} \left( f(x,t,u,\frac{\partial u}{\partial x}) \right) s(x,t,u,\frac{\partial u}{\partial x}) ]这里三个函数分别叫 (c)、(f)、(s)然后在代码里对应定义三个子函数。一开始我觉得这个格式很死板后来发现它把“扩散项、源项、时间系数”拆得清清楚楚反而有利于排查错误。2.4 步长和网格设计的黄金法则网格步长不是越小越好因为步长太小会显著增加计算量而步长太大则会导致精度不足甚至数值发散。对一个抛物型方程如果使用显式格式推进时间时间步长 (\Delta t) 和空间步长 (\Delta x) 之间必须满足稳定性限制典型条件形如[ \frac{\kappa \Delta t}{\Delta x^2} \le \frac{1}{2} ]这个条件叫 CFL 条件或扩散稳定性条件。我第一次写显式差分程序时完全没在意这个限制结果温度场越算越奇怪甚至出现了负温度最后查了很久才发现是步长匹配的问题。所以网格设计要综合考虑稳定性、精度和计算资源不能盲目加密。3. 手把手跑通一个模型热传导方程的 pdepe 实现看完前面那些基础我们来点能直接上手的。以经典的一维热传导方程为例目标区间是 (x \in [0,1])时间从 0 到 0.5扩散系数 (\kappa 1)初始温度分布为 (u(x,0) \sin(\pi x))两端边界温度固定为 0。这是一个教科书式的抛物型方程非常适合用 pdepe 跑通。3.1 三个子函数把数学方程翻译成 MATLAB 代码先给整段代码后面我再逐行解释关键点。新建一个脚本命名为heat_conduction_pdepe.m内容如下function heat_conduction_pdepe m 0; x linspace(0, 1, 101); t linspace(0, 0.5, 101); sol pdepe(m, pdefun, icfun, bcfun, x, t); u sol(:, :, 1); surf(x, t, u) xlabel(x) ylabel(t) zlabel(u) title(一维热传导方程的 pdepe 求解结果) end function [c, f, s] pdefun(x, t, u, DuDx) kappa 1; c 1; f kappa * DuDx; s 0; end function u0 icfun(x) u0 sin(pi * x); end function [pl, ql, pr, qr] bcfun(xl, ul, xr, ur, t) pl ul; ql 0; pr ur; qr 0; end3.2 逐段解释pdepe 到底在解什么pdepe的第一个参数m代表坐标系类型0 代表直角坐标1 代表柱坐标2 代表球坐标。我们这个问题是一维直角坐标所以填 0。接着传进去的四个函数句柄分别对应 PDE 系数、初始条件、边界条件。最后的x和t是网格点返回的sol是一个三维数组第一维是时间点第二维是空间点第三维是方程解的个数。pdefun里的写法非常直观方程的标准形式 (c \cdot u_t \partial_x(f) s)这里 (c1)(f \kappa u_x)(s0)。icfun就是初始条件返回一个值sin(pi*x)。bcfun里返回的四项分别代表左右边界的 (p) 和 (q)因为边界条件是 (p q \cdot f 0)对于固定温度边界 (u0)可以直接令 (pu)(q0)这样边界条件就变成了 (u0)。3.3 结果怎么验证解析解是最快的照妖镜这个热传导方程的解析解很容易求[ u(x,t) \sin(\pi x) \cdot e^{-\pi^2 \kappa t} ]我通常会在代码里加一个对照部分把解析解算出来然后算最大绝对误差u_analytical sin(pi * x) * exp(-pi^2 * t); error max(max(abs(u - u_analytical))); disp([最大绝对误差: , num2str(error)])实测下来pdepe 在默认容差下精度很不错最大误差通常在 (10^{-4}) 量级。这一步非常关键一旦代码有小错误对比解析解能第一时间发现而不是等画出来的图变得稀奇古怪才去排查。4. 踩过的坑与调试笔记从数值振荡到矩阵运算代码跑通只是第一步真正花时间的是调试。这里我把遇到过的几个典型问题展开讲讲全是实际操作中容易踩的地方。4.1 显式差分振荡与稳定性条件用 pdepe 很少遇到显式差分那种“炸场”但自己写显式格式时稳定性问题几乎是绕不开的。我之前写过一段求解扩散方程的显式差分代码核心迭代就三行u(2:end-1) u(2:end-1) ... kappa * dt / dx^2 * (u(3:end) - 2*u(2:end-1) u(1:end-2));如果 ( \kappa dt / dx^2 0.5 )数值解就会开始出现高频振荡幅度越来越大最后全是 NaN。这不是代码逻辑错误而是格式本身就发散。解决办法有两条把时间步长压缩到稳定范围内或者改用隐式格式。隐式格式没有这个稳定性限制是实际工程中更稳妥的选择。4.2 点乘和直接乘的区别一个不起眼但致命的错误很多刚接触 MATLAB 的人会搜“matlab 中点乘和直接乘有啥区别”。在偏微分方程数值解中这个问题特别容易跳出来咬你一口。比如在有限差分迭代里如果某个系数向量要和网格上的值逐点相乘必须用.*如果误写成*MATLAB 就会按照矩阵乘法去计算轻则报维度错误重则矩阵形状碰巧匹配但结果完全错了。我印象最深的一次是在求解带变系数的扩散方程时把kappa .* (u(i1)-u(i))写成了kappa * (u(i1)-u(i))而kappa恰好是一个向量另一个是向量MATLAB 没报错因为两者的维度正好匹配但计算结果完全不对整个浓度场都扭曲了。从那次以后我在涉及逐点运算的代码里都会特意检查乘号前面有没有点。4.3 pdepe 报错和边界条件函数的分工pdepe 报错常见的有几类。一类是“Spatial discretization has failed”大概率是边界条件函数返回的pl、ql、pr、qr某些项写错了导致离散后的 Jacobian 矩阵奇异。还有一类是求解到中途出现NaN一般是系数函数pdefun里传入了u和DuDx结果你在里面做了类似1/u的操作而某些时刻u会正好等于 0。调试边界条件时我习惯先把左右边界条件都简化成 Dirichlet 形式跑通之后再逐步加复杂条件。这样能把问题范围缩小不会一上来就是一团乱麻。4.4 验证数值解的几板斧写偏微分方程求解代码我基本都会做这几步验证找解析解对比没有解析解就想办法构造一个带有已知解的“人造源项”。检查守恒量。比如热传导方程无源无热交换时若边界是绝热的总热量应该不随时间变化。做网格收敛性测试步长减半看误差是否按预期比例下降。这三个方法在绝大多数情况下都能把隐蔽的 bug 挖出来。我自己的体会是越早做验证后面返工越少。5. 进阶用 PDE Toolbox 和自写有限差分应对更复杂的问题pdepe 虽然好用但它主要解决一维问题。遇到二维、三维或者不规则几何区域时就要升级工具了。5.1 PDE Toolbox 适合的场景PDE Toolbox 是 MATLAB 专门用来求解 PDE 的图形化工具箱支持有限元方法能处理二维和三维区域上的椭圆方程、抛物方程、双曲方程也支持多个方程联立。对于传热、结构力学、电磁场这类问题我通常直接用工具箱里预定义的方程模块来建模型设置好几何区域、材料参数、边界条件和初始条件然后生成网格并求解。尤其当几何区域复杂时手写网格生成和组装有限元矩阵是非常费时的事情PDE Toolbox 能把这部分时间省下来。5.2 什么时候值得自写有限差分自写差分最大的优势是透明、可控适合教学演示也适合处理一些工具箱覆盖不到的非标准方程。比如带上强非线性源项或者需要人为控制离散格式细节的时候写自己的迭代代码往往比强行套工具箱更灵活。对于一维问题我习惯先写一个显式模板做快速试探确认方程形态没问题后再改写成隐式格式用于正式计算。5.3 一维扩散方程的隐式差分模板这里给一个最常用的隐式差分模板使用 Crank-Nicolson 格式兼顾精度和稳定性function u implicit_diffusion_1d(kappa, L, T, Nx, Nt) dx L / (Nx - 1); dt T / Nt; u zeros(Nt, Nx); u(1, :) sin(pi * linspace(0, L, Nx)); r kappa * dt / dx^2; A diag((1 r) * ones(Nx, 1)) ... diag(-r/2 * ones(Nx-1, 1), 1) ... diag(-r/2 * ones(Nx-1, 1), -1); A(1, :) 0; A(1, 1) 1; A(end, :) 0; A(end, end) 1; for n 1:Nt-1 b u(n, :).; b(2:end-1) b(2:end-1) ... r/2 * (u(n, 3:end) - 2*u(n, 2:end-1) u(n, 1:end-2)).; b(1) 0; b(end) 0; u(n1, :) (A \ b).; end end这个模板虽然不复杂但已经可以处理很多实际扩散问题了。如果你要处理非线性项只需要在组装A或计算右端项时加入对应的线性化处理。5.4 后续可以继续扩展的方向偏微分方程本身不是终点算出来的结果经常还要对接后续分析。比如用lsqcurvefit根据实验观测数据反演方程中的未知扩散系数或者把 PDE 模型的输出和卷积、系统辨识等方法结合用于更复杂的动态系统建模。这些方向需要单独写文章来讲但核心还是先把 PDE 数值解做稳后面才有底气做反向题。我自己这些年用 MATLAB 解 PDE最大的体会是“先跑通再优化”。第一版代码往往不是最高效的但它能让你快速知道方程行为、参数范围和潜在问题。等一切都在掌控之中再去考虑更精细的格式、更快的算法或者更复杂的模型完全来得及。本文还有配套的精品资源点击获取