24行代码教你学会有限元编程
MATLAB 有限元 Poisson 方程 三角形 元 稀疏矩阵 向量化装配 数值 PDE
挑战最短有限元代码：24 行 MATLAB 写完一个二维 P1 FEM
这篇文章不准备比谁更会把代码挤成一团。
如果允许
a=...;b=...;c=...;d=...;
这样无限往一行里塞语句，那么任何 MATLAB 程序理论上都可以写成“一行”，行数也就失去了意义。
所以先把规则说清楚。
代码规则
解一个二维有限元问题：
取精确解
因此右端项为
挑战代码需要同时完成下面这些事情：
必须是二维有限元。 使用三角形 元。 网格必须在脚本里生成。 不允许读取提前准备好的 mesh 文件。 必须真正组装有限元矩阵。 不能调用 PDE Toolbox，也不能调用外部 FEM 软件包。 不能调用自己提前写好的辅助函数。 一个 .m文件复 制后即可运行。必须包含完整流程。 网格、单元几何、刚度矩阵、载荷向量、Dirichlet 边界、线性方程求解都不能省。 必须有结果检查。 至少计算一次与精确解的误差。 必须有可视化。 空行和注释不计入行数，但不允许靠一行塞多个独立语句作弊。 代码按 MATLAB R2016b 及以后版本书写，允许使用隐式数组扩展。
先看完整代码
下面这段直接保存为 fem24.m 即可运行。
n=30;
[X,Y]=meshgrid(linspace(0,1,n+1));
p=[X(:),Y(:)];
t=delaunay(p(:,1),p(:,2));
N=size(p,1);
x=p(:,1);
y=p(:,2);
xt=x(t);
yt=y(t);
A=abs((xt(:,2)-xt(:,1)).*(yt(:,3)-yt(:,1))-(xt(:,3)-xt(:,1)).*(yt(:,2)-yt(:,1)))/2;
b=[yt(:,2)-yt(:,3),yt(:,3)-yt(:,1),yt(:,1)-yt(:,2)];
c=[xt(:,3)-xt(:,2),xt(:,1)-xt(:,3),xt(:,2)-xt(:,1)];
I=kron(t,ones(1,3));
J=repmat(t,1,3);
V=(kron(b,ones(1,3)).*repmat(b,1,3)+kron(c,ones(1,3)).*repmat(c,1,3))./(4*A);
K=sparse(I(:),J(:),V(:),N,N);
q=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;
f=2*pi^2*sin(pi*q(:,1)).*sin(pi*q(:,2));
F=accumarray(t(:),repmat(A.*f/3,3,1),[N,1]);
free=find(x>0 & x<1 & y>0 & y<1);
u=zeros(N,1);
u(free)=K(free,free)\F(free);
trisurf(t,x,y,u,'EdgeColor','none');
title(sprintf('max nodal error = %.3e',norm(u-sin(pi*x).*sin(pi*y),inf)));
24 行里，包含了每个个有限元最基本的环节。
当然，这不是工程代码。它更像是把一个有限元程序剥到只剩骨架，然后看看那些我们平时习惯写成几十个函数、几百行代码的东西，究竟哪些才是不可删除的核心。
24 行代码解读
有限元离散从弱形式开始。
对于任意 ，Poisson 方程写成
取三角剖分 ，令
离散问题就是求 ，使
最后得到熟悉的线性系统
24 行代码核心要做的事情，无非就是把这里的 和 真正造出来。
第一部分：网格其实只需要两个数组
[X,Y]=meshgrid(linspace(0,1,n+1));
p=[X(:),Y(:)];
t=delaunay(p(:,1),p(:,2));
有限元网格最核心的数据结构只有两样：
p：节点坐标；t：每个三角形由哪三个节点组成。
如果 t(k,:)=[17,18,49]，它的意思就是第 个单元的三个顶点分别是全局节点 17、18、49。
后面所有局部到全局的装配，本质上都在使用这个映射。
这一点很值得品味。
有限元程序看起来有形函数、积分、矩阵、边界条件一大堆概念，但真正进入计算机以后，拓扑关系最终就是一个整数矩阵 t。
第二部分：三角形 P1 元为什么特别适合凝练代码
设一个三角形 的三个顶点为
定义
以及
对线性三角形元，三个局部基函数的梯度在单元内部是常数，并且
所以局部刚度矩阵直接就是
代码里的
A=abs((xt(:,2)-xt(:,1)).*(yt(:,3)-yt(:,1))-(xt(:,3)-xt(:,1)).*(yt(:,2)-yt(:,1)))/2;
b=[yt(:,2)-yt(:,3),yt(:,3)-yt(:,1),yt(:,1)-yt(:,2)];
c=[xt(:,3)-xt(:,2),xt(:,1)-xt(:,3),xt(:,2)-xt(:,1)];
就是把所有三角形的 、 、 一次算完。
这里没有显式写出任何形函数。不是因为形函数被忽略了，而是因为对于 三角形，形函数梯度的全部信息已经压缩进了 b 和 c。
❝短代码最有意思的地方，不是少敲几个字符，而是找到一个数据表达，使数学结构自动包含在数组里。
第三部分：真正有意思的是这三行装配
整段代码里，最重要的应该是下面三行：
I=kron(t,ones(1,3));
J=repmat(t,1,3);
V=(kron(b,ones(1,3)).*repmat(b,1,3)+kron(c,ones(1,3)).*repmat(c,1,3))./(4*A);
随后一句
K=sparse(I(:),J(:),V(:),N,N);
完成全局刚度矩阵装配。
为什么？
每个三角形有三个自由度，所以一个局部刚度矩阵有 个条目。
假设某个单元的全局编号是
那么这 9 个条目需要加到全局位置
假设一个三角形单元的三个全局节点是：
它的 局部刚度矩阵要装到全局位置：
MATLAB 中 kron(A,B) 表示 Kronecker 积。这里I=kron(t,ones(1,3));会把每个节点编号重复 3 次。对 t=[5 8 12]，得到 I=[5,5,5,8,8,8,12,12,12]。 而J=repmat(t,1,3);得到 J=[5,8,12,5,8,12,5,8,12]。 于是 (I,J) 正好列出局部 矩阵的 9 个全局坐标。
对于一般的单元数据t,kron(t,ones(1,3)) 负责生成所有行指标，其具体函数功能参见克罗内克张量积。repmat(t,1,3) 负责生成所有列指标。V 则一次性生成对应的所有局部刚度条目。
最关键的一点是，sparse 在遇到重复的 (I,J) 时会自动把对应数值相加。 而有限元所谓的 assembly，不就是“不同单元共享同一个全局自由度时，把贡献加起来”吗？ 于是原本教科书里需要专门讲一节的局部到全局装配，在 MATLAB 里几乎变成了一次稀疏矩阵构造。
编者认为：这比单纯少写几行更值得看。它体现的是 MATLAB 这种数组语言和有限元数据结构之间非常自然的一次对接。是 MATLAB 向量化编程范式中非常具有美感的一个经典例子。
第四部分：载荷向量也没有逐单元循环
右端项为
这里为了保持代码简洁，在每个三角形上使用重心一点积分。
设 是三角形重心，则近似有
代码先算所有单元重心：
q=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;
再计算右端函数：
f=2*pi^2*sin(pi*q(:,1)).*sin(pi*q(:,2));
最后：
F=accumarray(t(:),repmat(A.*f/3,3,1),[N,1]);
accumarray 在这里干的事情和刚才 sparse 的重复索引累加非常像，只不过这里是一维的。
每个三角形给自己的三个顶点各贡献一份 ，共享节点收到的所有单元贡献自动相加。 所以刚度矩阵用 sparse 装，右端载荷向量用 accumarray 装。下面的动态演示图直观地展示了 accumarray 的运行逻辑：
两个 MATLAB 原生函数，正好对应有限元里最频繁出现的“局部贡献汇总”。
第五部分：Dirichlet 边界只保留自由自由度
对于齐次 Dirichlet 条件，边界节点上的解就是零。
于是没有必要先修改矩阵行列，再往对角线上塞 1。
直接找内部节点：
free=find(x>0 & x<1 & y>0 & y<1);
然后解
u(free)=K(free,free)\F(free);
就结束了。
数学上对应的其实就是把线性系统限制到自由自由度子空间。
这也是短代码里另一个很好的取舍：如果问题本身已经是齐次边界，就不要为了“通用写法”引入不需要的代码。
24 行为什么还能保持完整
把代码重新拆开看，它其实就是五件事：
p t | ||
A b, c | ||
I J, V, sparse | ||
q f, accumarray | ||
free K\F |
平时我们使用 FreeFEM、FEniCS、deal.II、MFEM，或者自己维护一个成熟 FEM 框架时，有限元的基本步骤很容易被封装层遮住。求解实际问题是方便了，但是对于初学者理解代码做了什么就比较困难。
把程序压缩到二十多行以后，反而看得很清楚：
有限元最基本的计算骨架就是这些。后面更复杂的方程，更复杂的有限元空间编程基本都是这个思路，万变不离其宗。
还能更短？！
首先声明代码短并不自动等于代码好。 比如下面几种做法当然还能继续减少“物理行数”：
用 deal把多个赋值合并；把很多独立语句用分号塞到同一行； 把右端项、边界和网格全部硬编码； 省掉误差检查； 省掉画图； 调用一个提前写好的 assemble()；直接写出结构化网格对应的五点差分矩阵。
这些都能让代码看更短。但最后一条甚至已经不太算是在写一个真正的三角形有限元程序了，更像是“作弊”。
这次挑战想保留的是：
数学结构不丢，程序骨架不丢，然后再看 MATLAB 自己能帮我们省掉多少机械工作。所以我更愿意把 sparse、accumarray、向量化单元几何这些叫作“压缩”，而不是简单把十条命令挤到一行。
为什么不用单元循环
很多有限元入门代码都会写成
for elem=1:NT
...
Ke=...
K(nodes,nodes)=K(nodes,nodes)+Ke;
end
这个写法非常直观，教学上完全没有问题，甚至概念更加清晰便于理解。
但 MATLAB 里直接在循环中反复修改一个大稀疏矩阵，通常不是最舒服的做法。
更自然的思路是：
一次算出所有单元的 9 个矩阵条目； 一次生成所有全局行列索引； 最后只调用一次 sparse。
也就是说，不要把思路写成
“第一个单元装一次，第二个单元再装一次……”
而是改成
“把所有局部条目先列出来，然后一次汇总。”
这不只是为了少几行代码。它实际上把程序从“单元驱动的逐次更新”改成了“批量生成 COO 数据，再构造 CSR/CSC 稀疏矩阵”的思路。 规模真正上去以后，这种数据组织方式也更接近高性能有限元装配的基本思想。
这里代码没有考虑的事情
需要指出，这份程序没有处理：
一般多边形区域； 非齐次 Dirichlet 条件； Neumann 或 Robin 边界； 高阶积分公式； 变系数扩散； 反应项、对流项； 或更高阶有限元； 自适应网格； 后验误差估计； 网格质量检查； 大规模迭代求解器与预条件器。
尤其是载荷项这里只用了三角形重心的一点积分。对于当前演示足够，但如果要写通用有限元程序，数值积分必须被认真设计。
成熟有限元软件之所以长，是因为它要处理通用性、鲁棒性、效率、可维护性和各种边界情况。 短代码的价值，是帮我们确认这些复杂性到底建立在什么最小核心之上。
可以思考以下推广
加入非齐次 Dirichlet 条件｜最少增加几行
把向量化装配改回单元循环｜比较速度
把 P1 换成 P2｜看看代码会在哪里变长
最后
好的教学代码，不该把数学藏起来。
有限元代码一旦写到上万行，最容易忘掉的是矩阵里每一个非零元到底从哪里来。而在这 24 行里，这件事几乎没有可以删去的地方。 三角形给出局部梯度，局部梯度给出 9 个刚度条目，t 把局部编号映射成全局编号，sparse 把共享节点的贡献加起来。然后解方程。就这么多。
也正是在这种时候，才更容易体会到计算数学的魅力：那些写在公式里的抽象结构，并不是停留在纸面上的符号，它们可以一步步落到代码里，最后变成一个真正可以计算、可以观察、可以验证的数值结果。所以老师们总强调程序基本功，大概也不只是为了“把程序写出来”。更重要的是，当代码足够熟练之后，我们才有机会透过实现本身，重新看见背后的数学。
参考文献
J. Alberty, C. Carstensen, S. A. Funken, Remarks around 50 lines of Matlab: short finite element implementation, Numerical Algorithms, 20, 117–137, 1999.
C. Carstensen, D. Gallistl, J. Hu, A discrete Helmholtz decomposition with Morley finite element functions and the optimality of adaptive finite element schemes, Computers & Mathematics with Applications, 68(12), 2167–2181, 2014. DOI: 10.1016/j.camwa.2014.07.019.
O. J. Sutton, The Virtual Element Method in 50 lines of MATLAB, Numerical Algorithms, 75, 1141–1159, 2017.
S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Springer, 2008.
