企业官网建设流程全解析

简介：一份围绕有限元弱形式主题的 doc 学习文档，面向有限元初学者、物理仿真研究人员及工程分析人员。内容系统梳理 PDE 弱形式的基本概念、三种物理问题描述方式（偏微分方程、能量最小化形式、弱形式）之间的联系，并结合弹性静力学 Navier 方程与弹性能量表达式，展示从 PDE 到泛函变分再到有限元离散求解的完整思路。文档还特别说明弱形式在非线性多物理场问题中的优势，以及 COMSOL、SOL Multiphysics、ANSYS 等软件中的应用场景，可帮助读者理解有限元底层数学基础，并为使用商业软件的自定义 PDE 建模提供参考。资料压缩包共 1 个 doc 文件，大小 649KB，内容集中、便于快速通读。已有 128 人学习，适合作为研究生课程补充或工程师入门有限元弱形式的精炼材料。

1. 有限元的弱形式：从“偏微分方程解不出来”到“矩阵方程”

很多做数值计算的人第一次被迫面对“弱形式”这三个字，是在一份名为《有限元的弱形式.doc》的讲义里：前几页还在规规矩矩推导偏微分方程，翻过一页突然冒出两个积分和一堆带下标的函数空间，再往后就是刚度矩阵、载荷向量和网格剖分。如果你也卡在这个位置上，先把最反直觉的结论说出来：弱形式不是把方程“变弱”，而是把偏微分方程的逐点约束换成积分约束。换完之后，解的光滑性要求从“二阶连续可导”降到“一阶导数平方可积”，代价是方程从单点严格满足变成在任意测试函数的加权意义下满足。这个交换的直接收益是：分片多项式现在可以进入解空间，而分片多项式正是有限元刚度矩阵能组装起来的前提。做matlab有限元编程求解实例时，第一步永远是从一维泊松方程开始写弱形式，再把组装出来的矩阵打印出来和手算结果核对。这篇文章就把这条路完整走一遍：先推弱形式，再写最小可运行代码，最后用误差范数验证程序没写错。适合刚接触有限元的工程师，也适合边界条件总处理不干净的从业者。

2. 强形式到弱形式：分部积分与测试函数怎么选

2.1 先写出强形式，再乘一个测试函数

以一维稳态热传导问题为例，控制方程是

-d/dx(k du/dx) = f，x ∈ (0, 1)

边界条件取 u(0) = u0，k du/dx |_{x=1} = g。

所谓“强形式”，字面意思是方程在定义域内每一个点 x 上都要严格成立。这要求 u 在闭区间内二阶连续可导。对实际工程中的温度场、位移场来说，这个要求相当苛刻：两种材料交界面处温度连续但温度梯度可以突变，集中载荷作用点附近位移的导数也会间断。强形式在数学上仍然成立，但直接求解析解或者构造数值格式都不容易。

有限元的思路是：不要求方程逐点满足，而是先选一个测试函数 v(x)，把强形式两边同时乘 v，再对全区间积分。测试函数可以理解成一把尺子，方程不再逐点测量，而是用无数把尺子去量“总体偏差”。当 v 取遍所有足够光滑且在边界上满足特定条件的函数时，加权等式与原方程等价。这一步之后，微分方程变成了积分方程，允许被积函数出现有限个跳跃点，因为跳跃点测度为零，不影响积分值。

2.2 分部积分：导数转移给测试函数

对等式 ∫(-d/dx(k u')) v dx = ∫ f v dx 的左端做分部积分：

∫ k u' v' dx - [k u' v]₀¹ = ∫ f v dx

边界项展开是 k u'(1) v(1) - k u'(0) v(0)。此时要区分两类边界信息：右端 x=1 给出的是 k u'(1) = g，属于已知自然边界条件，可以直接把 g 代入；左端 x=0 给出的是 u(0) = u0，但 k u'(0) 是未知量，这一项必须被消掉。

消掉的办法是约束测试函数 v(0)=0。本质边界条件上的测试函数取零，是弱形式理论里最容易被忽略的规则。取了 v(0)=0 后，弱形式的最终形式是：

求 u，使得 u(0)=u0 且 u、u' 平方可积；对所有满足 v(0)=0 且 v、v' 平方可积的 v，有

∫ k u' v' dx = ∫ f v dx + g v(1)

右端最后一项 g v(1) 就是 Neumann 边界条件进入有限元方程的入口。如果 g=0，这一项自然消失，不需要做任何额外处理。这就是“自然”两个字的含义：它自己掉进弱形式里，而不是被人为塞进去的。

2.3 “弱”在哪里：函数空间与解的光滑性

强形式要求 u ∈ C²，弱形式只要求 u ∈ H¹。H¹ 是 Sobolev 空间，包含所有函数本身平方可积、一阶导数也平方可积的函数。分片线性函数正好落在 H¹ 里：它在每个单元内是线性函数，导数在单元内为常数，在单元边界处跳跃，但跳跃值有限，平方积分有限。

“弱”不代表计算精度差。它放宽的是逐点要求，换来的是数值可操作性。你选一个分片线性的 v_h，代入弱形式，得到的就是关于节点值的代数方程组。理解这个层次后，那些带下标的函数空间定义就不再是障碍，而只是说明了“允许哪些函数进来”。

2.4 用具体函数验证分部积分边界项

弱形式推导容易在边界项符号上出错。常见做法是用一组具体函数核对恒等式，在 MATLAB 里跑通后再写正式程序：

% 用具体函数验证分部积分恒等式 % ∫ -u'' v dx = ∫ u' v' dx - [u'v]_0^1 syms x u = x^2; % 测试函数 u v = x * (1 - x); % 满足 v(0)=0, v(1)=0 的测试函数 L = 1; I1 = int(-diff(u, x, 2) * v, x, 0, L); % 原始弱形式左端 I2 = int(diff(u, x) * diff(v, x), x, 0, L) - ... (subs(diff(u, x) * v, x, L) - subs(diff(u, x) * v, x, 0)); simplify(I1 - I2) % 输出 0 说明分部积分推导正确

逻辑说明：I1 是原方程左端乘测试函数后的积分，I2 是分部积分后去掉边界项的结果。对满足 v(0)=v(1)=0 的测试函数，边界项subs(diff(u,x)*v, x, L)和subs(diff(u,x)*v, x, 0)都应该为 0，但代码里保留完整形式，是为了验证符号运算本身没出错。把 u、v 换成分段函数也能做同样验证，只是符号积分可能变慢。这个习惯值得保留：任何弱形式推导，先用具体函数验边界项，再进离散。

3. 有限元离散：线性基函数下的刚度矩阵与载荷向量

3.1 从弱形式到线性方程组

弱形式仍然是无限维问题，因为函数空间 V 中有无穷多个函数。有限元离散做的事情是把 V 换成一个有限维子空间 V_h，比如所有在节点 0 到 N 上取值、在单元内部是一次多项式的连续分片线性函数。这个空间里的任意函数完全由 N+1 个节点值 U₀, U₁, ..., U_N 决定。

用形状函数 φ_i(x) 描述第 i 个节点上的“尖塔”：φ_i(x_i)=1，在相邻两个单元内线性降到 0，其他位置全为 0。V_h 中任意函数都能写成 u_h(x) = Σ U_j φ_j(x)，代入弱形式并依次取 v = φ_i，得到 N+1 个代数方程：

K_ij = ∫ k φ_i' φ_j' dx，F_i = ∫ f φ_i dx + g φ_i(1)

这就是 K U = F 的来源。刚度矩阵 K 是稀疏的，因为 φ_i 和 φ_j 的支撑区间只在 i、j 相同或相邻时相交，每个单元只对相邻自由度有贡献。这个局部性决定了有限元组装的“单元循环”方案：逐个单元计算局部矩阵，再按节点编号投放进全局矩阵。

3.2 单元刚度矩阵：为什么要记住 1 -1 -1 1

在单个单元 e 上，单元长度为 h，局部节点 1 和 2 对应全局节点 e 和 e+1。两个线性形状函数的导数为 φ₁' = -1/h，φ₂' = 1/h，代入 K_ij 的积分，k 取常数时得到局部刚度矩阵：

ke = (k/h) * [[1, -1], [-1, 1]]

这个 2×2 矩阵是整个一维有限元最常用的公式。物理意义可以从两个角度看：每行元素之和为零，对应刚体平移模式下单元内没有应力；对角线元素为正，保证系统正定。组装时，局部节点 1 的贡献累加到全局第 e 行、第 e 列，局部节点 2 的贡献累加到第 e+1 行、第 e+1 列。

载荷向量的单元贡献取决于 f 的积分。最简单的近似是两点梯形积分：fe = h/2 * (f(x_e) + f(x_{e+1}))。f 是线性函数时这个公式精确；f 是非线性函数时它是一阶近似。先把梯形积分跑通框架，再升级成高斯积分，是更稳的推进节奏。

3.3 matlab有限元编程求解实例：组装与求解的最小代码

下面是一段完整的一维 matlab有限元编程求解实例代码，对应方程 -u'' = π² sin(πx)，边界条件 u(0)=u(1)=0，解析解 u = sin(πx)。k 取 1，g 取 0。

% 一维泊松方程 -u'' = f, u(0)=u(1)=0 % 解析解 u = sin(pi*x) f = @(x) pi^2 .* sin(pi * x); N = 20; % 单元数 x = linspace(0, 1, N + 1)'; % 节点坐标 h = 1 / N; % 单元长度 K = sparse(N + 1, N + 1); % 全局刚度矩阵 F = zeros(N + 1, 1); % 全局载荷向量 for e = 1:N n = [e, e + 1]; % 当前单元对应的全局节点 ke = (1 / h) * [1, -1; -1, 1]; % 单元刚度矩阵，k=1 fe = h / 2 * [f(x(e)); f(x(e + 1))]; % 单元载荷向量，梯形积分 K(n, n) = K(n, n) + ke; % 投放进全局矩阵 F(n) = F(n) + fe; % 投放进全局向量 end % 本质边界条件 u(0)=0, u(1)=0，只求解内部自由节点 free = 2:N; U = zeros(N + 1, 1); U(free) = K(free, free) \ F(free); % 与解析解对比前 4 个节点 ue = sin(pi * x); disp([x(1:4), U(1:4), ue(1:4)]);

逻辑说明：每个单元独立计算局部刚度矩阵和局部载荷向量，再依据全局节点编号把贡献“累加”进 K 和 F。sparse在 N 增大时节省大量内存，避免全稠密矩阵。free = 2:N排除左端节点 1 和右端节点 N+1，因为这两个节点值已知为零；只对内部节点求解，解出来后与端节点上的零合并成完整 U。

参数说明：N 控制网格密度，N 增大时近似解更接近真解，但二维问题自由度会平方增长；f 是右端热源项，以函数句柄传入可以让代码不绑死具体表达式；h 是网格尺寸，初始化后先打印确认等于 1/N，能排查网格划分错误。U(free) = K(free, free) \ F(free)用的是 MATLAB 内置稀疏直接求解器，一维问题 N 在 10 万以内都很快。

3.4 载荷向量升级：两点高斯积分

梯形积分对光滑但非线性的 f 会有可感知的误差。常见做法是换成两点高斯-勒让德积分，积分点取 ±1/√3，权重取 1。把参考区间 [-1, 1] 映射到物理单元 [x_e, x_{e+1}] 后，单元中点为 x_m = (x_e + x_{e+1})/2，高斯点位于 x_m ± h/(2√3)。

替换前面代码里的 fe 行：

xm = (x(e) + x(e + 1)) / 2; s = h / (2 * sqrt(3)); fe = h / 2 * [f(xm - s); f(xm + s)];

逻辑说明：两点高斯积分对三次多项式以内的 f 精确，梯形积分只对一次多项式精确。高斯点上的 f 值经过雅可比因子 h/2 映射到物理单元。对大多数光滑右端项，两点高斯足够；如果 f 含强烈变化，先加密网格，再考虑增加高斯点数量，而不是只靠单元内加密。

4. 弱形式中的边界条件：Dirichlet与Neumann的两张面孔

4.1 本质边界条件与自然边界条件的自由度差异

回到弱形式，边界条件分两类：一类直接给定 u 的值，叫本质边界条件或 Dirichlet 条件；一类给定导数组合，叫自然边界条件或 Neumann 条件。本质边界条件必须进入解空间：u_h 在对应节点上取给定值，测试函数在对应节点上取 0。自然边界条件不需要改解空间，它通过弱形式中的边界积分直接进入右端项。

数值上经常出错的点有两个：一是在 Neumann 边界上错误地强制节点为零，二是在 Dirichlet 边界上忘记做任何处理。判断方法很简单：看弱形式边界积分项里 u 是否已知。u 已知，是本质边界，要改矩阵；u 未知但它的导数组合已知，是自然边界，只加右端项。

4.2 直接消去法：把已知节点值移到右边

非齐次 Dirichlet 边界 u(0)=a、u(1)=b 是最常见的情况。直接消去法的思路是：先把已知节点值写入解向量，再把已知节点对未知节点方程的影响移到右端，最后只对自由节点求解。

bnodes = [1; N + 1]; % 本质边界对应的全局节点编号 bvals = [a; b]; % 边界节点值 U = zeros(N + 1, 1); U(bnodes) = bvals; % 写入已知节点值 F2 = F - K * U; % 已知节点贡献移到右端 free = setdiff(1:N + 1, bnodes); U(free) = K(free, free) \ F2(free);

逻辑说明：完整方程是 K U = F，但 U 中部分分量已知。把 K 与已知 U 相乘，得到的是边界节点对自由节点方程贡献的“力”，移到右端得到 F2。如果只取 K(free, free) 和 F(free)，边界节点影响被完全丢掉，结果会明显偏离真解。setdiff适合教学和小规模问题；大规模运算建议在组装前完成自由度编号，避免每次构造 free 索引。

直接消去法的优点是精确：边界值以机器精度进入解。缺点是维护自由度索引稍微麻烦。自适应加密下节点编号会变化，边界节点和自由节点必须分开记录。

4.3 罚函数法：一个参数决定收敛性

罚函数法把 Dirichlet 条件作为强惩罚项放进弱形式，相当于在边界上加了很大的弹簧。离散后的修改非常简单：

C = 1e8; % 惩罚系数，经验值取刚度矩阵最大对角元的 1e6 到 1e8 倍 for i = 1:numel(bnodes) K(bnodes(i), bnodes(i)) = K(bnodes(i), bnodes(i)) + C; F(bnodes(i)) = F(bnodes(i)) + C * bvals(i); end U = K \ F;

惩罚系数 C 选小了，边界约束不足，表现为边界节点值偏离给定值；选大了，矩阵条件数迅速恶化，迭代求解器可能不收敛。经验准则是让 C 与刚度矩阵对角最大元素保持至少 6 个数量级差距，然后检查abs(U(bnodes) - bvals)是否在可接受范围。实际工程里直接消去法更可控，罚函数法适合快速原型实现。

4.4 自然边界条件：往右端项加值的顺序不能错

右端 Neumann 边界 k u'(1)=g，弱形式载荷向量最后一个分量需要加 g：F(end) = F(end) + g。原因是最末端节点的形状函数 φ_N 在该点取值为 1，边界积分 g v(1) 在测试函数取 φ_N 时直接落到最后一个自由度上。

有一个顺序问题容易忽略：如果同一节点上同时存在 Dirichlet 和 Neumann 条件，必须先加自然边界贡献，再做本质边界消去。先消去本质边界会把 Neumann 贡献错误地带进自由节点方程，结果整个右端项都会偏。推荐把“先自然、后本质”作为任何有限元程序里的固定顺序，写进注释里。

5. 用误差范数验证弱形式实现：MATLAB算例与网格收敛

5.1 构造有解析解的光滑问题

验证弱形式程序是否写对的可靠方式，是设计一个已知解析解的问题。取 u(x) = sin(πx)，则 -u'' = π² sin(πx)。在代码里设 f = @(x) pi^2 .* sin(pi*x)，k=1，边界条件 u(0)=u(1)=0。这个解足够光滑，没有角点奇异性干扰，误差只来自网格离散和数值积分，方便观察收敛阶。

5.2 用单元高斯积分计算误差范数

误差计算不要用 interp1 对数值解重新采样，那样会混入插值误差。更干净的做法是在每个单元上用两点高斯积分直接计算精确解和数值解之差：

err_L2 = 0; err_grad = 0; for e = 1:N xm = (x(e) + x(e + 1)) / 2; s = h / (2 * sqrt(3)); xg = [xm - s, xm + s]; % 高斯点 ug = [U(e), U(e + 1)]; % 数值解在高斯点的线性插值 ue = sin(pi * xg); % 精确解 ue_x = pi * cos(pi * xg); % 精确解的导数 grad_uh = (U(e + 1) - U(e)) / h; % 单元内数值解导数为常数 err_L2 = err_L2 + h / 2 * sum((ug - ue).^2); err_grad = err_grad + h / 2 * sum((grad_uh - ue_x).^2); end err_L2 = sqrt(err_L2); err_grad = sqrt(err_grad); err_H1 = sqrt(err_L2^2 + err_grad^2);

逻辑说明：高斯点上的数值解由该单元两个节点值线性插值得到，精确解直接用解析公式。误差平方乘上雅可比因子 h/2 后累加，最后开方得到 L2 范数。梯度误差部分利用了线性元导数为常数的特性，数值导数不需要额外差分。线性元的理论收敛阶为 L2 误差约 O(h²)，H1 误差约 O(h)。加密网格时可以按下面表格核对：

5.3 程序能跑通的三个快速检查

第一个检查刚度矩阵对称性：norm(K - K', 'fro')应该为 0。组装索引写错时，经常会丢失非对角线的一半贡献，这个语句能立刻暴露问题。第二个检查刚体模式：在纯 Neumann 问题里把 U 全部置 1，计算K * ones(N+1,1)应该接近 0；带 Dirichlet 边界时内部行的结果仍接近 0，边界行会有响应。第三个检查求解残差：res = norm(F - K * U)，如果消去边界条件时只是把行清零而不是缩减矩阵，残差里会出现边界行的巨大数值。

N 超过 10 万后，一维三对角刚度矩阵的直接求解仍然很快，但推广到二维三角形网格时带宽增大，内存开销会明显上升。常见做法是切换成共轭梯度法加不完全 Cholesky 预条件：pcg(K(free,free), F2(free), 1e-8, 400, ichol(K(free,free)))。此时如果还在用罚函数法处理 Dirichlet 边界，巨大对角线会破坏预条件子，收敛会非常慢。弱形式理论看到这里，你会明白直接消去法虽然代码多一点，却是最可预测的方案。把 K 和 F 的维度与网格节点数对齐打印出来，能快速定位组装漏单元的常见问题；下一步把一维代码推广到二维三角形单元时，弱形式本身不动，变的只是单元矩阵和积分规则。

本文还有配套的精品资源，点击获取