有限元弱形式推导与MATLAB编程实现：从微分方程到线性方程组

简介一份关于有限元弱形式的doc文档面向学习有限元方法的学生、科研人员及工程应用者系统梳理弱形式的核心知识点。内容从物理问题的三种描述方式偏微分方程、能量最小化形式、弱形式入手讲解三者的等价关系与各自适用范围并以弹性力学中的Navier方程和总势能泛函为例展示从PDE到泛函变分的推导过程。文档强调弱形式作为积分形式对积分变量连续性要求更低特别适合非线性及多物理场问题——这也是COMSOL Multiphysics等软件选择弱形式作为底层求解基础的重要原因有助于读者理解软件高级设置并拓展到传热、流体N-S方程等标准模块难以覆盖的场景。资源为单文件doc格式体积仅649KB内容紧凑便于离线阅读和快速查阅。目前已有128人学习下载适合作为有限元入门与进阶的补充笔记。1. 有限元弱形式从微分方程到线性方程组的转换关口接触有限元的第一道坎多半是“弱形式”。一个反直觉的事实是商业软件和自编程序真正求解的从来不是原始微分方程逐点成立的强形式而是它的积分形式。弱形式把“每个点都必须满足方程”的强约束放松成“在测试函数加权意义下整体满足”换来的是对解连续性要求的显著降低——二阶问题只用到一次导数分片线性单元就够用。下面的内容按我平时调程序的顺序推进先亲手推导从强形式到弱形式再做离散化和矩阵组装给出一个可直接运行的 MATLAB 有限元求解实例最后用制造解方法验证程序没有写错。适合写过一段有限元代码但对变分推导含糊的工程师也适合正在啃计算力学教材的研究生。2. 从强形式到弱形式加权残值与分部积分推导2.1 强形式的连续性要求二阶导数为什么难伺候一维问题的强形式写出来很紧凑− d/dx(c(x) du/dx) f(x)x ∈ (0, L)。加上边界条件后它要求在定义域内每一点都恒等成立。这隐含一个很高的门槛c(x)u′(x) 必须可导也就是 u 至少要有二阶连续导数u ∈ C²。对光滑单一材料这个要求不算过分可一旦遇到复合材料界面、层状结构或几何突变弹模 c(x) 不连续界面上通量 c u′ 连续而 u′ 发生跳跃强形式在这些点并不成立。真实工程问题里这种“不光滑”才是常态。弱形式换一个问法不要求残差处处为零只要求残差在一族测试函数权函数的加权积分意义下为零。二阶导数通过分部积分转移到测试函数上对解的要求降为“一阶导数平方可积”。分片线性函数的一阶导数是分片常数自然满足这个要求——这就是线性有限元成立的数学根基。多数教材从结论讲起初学者真正缺的是亲手推导一遍的过程。2.2 加权残值与 Galerkin 弱形式推导推导分三步。第一步把强形式写成残差形式 R(x) −(c u′)′ − f要求它与任意测试函数 v 的内积为零∫₀ᴸ R(x)v(x) dx 0。第二步对 ∫₀ᴸ(c u′)′v dx 做分部积分。这一步是整个弱形式的灵魂二阶导数从 u 身上“卸载”到 v 身上同时蹦出一个边界项∫₀ᴸ c u′v′ dx ∫₀ᴸ f v dx [c u′v]₀ᴸ第三步把边界项按边界条件归类通量 c u′ 被指定的边界直接代入u 被指定的边界令测试函数 v 在该处取零边界项自动消失。整理后得到标准弱形式求 u ∈ V使对任意 v ∈ W 有 a(u, v) F(v)其中 a(u,v) ∫ c u′v′ dxF(v) ∫ f v dx 边界通量项。分部积分这行最容易抄错可以用符号计算验一遍。下面这段 Python 用 sympy 随机取一组光滑函数验证恒等式 ∫(c u′)′v dx [c u′v] − ∫ c u′v′ dx 两边严格相等import sympy as sp x sp.symbols(x) c 1 x**2 # 变系数材料 u sp.sin(x) # 假定的解 v x**3 # 任取的测试函数 a, b 0, 1 # 积分区间 lhs sp.integrate(sp.diff(c*sp.diff(u, x), x) * v, (x, a, b)) rhs (c*sp.diff(u, x)*v).subs(x, b) - (c*sp.diff(u, x)*v).subs(x, a) \ - sp.integrate(c*sp.diff(u, x)*sp.diff(v, x), (x, a, b)) print(sp.simplify(lhs - rhs)) # 输出 0 即恒等式成立代码逻辑lhs 和 rhs 分别计算恒等式两侧subs代入上、下界得到边界项数值simplify对差值做化简。把 c、u、v 换成其他光滑函数结论不变。实际调试中如果这行符号检验过不去说明强形式或弱形式某一步写错了后面的矩阵组装做得再对也是白做。Galerkin 选择是让测试函数空间与试解空间用同一套基函数。直接收益是双线性形式 a(u,v) 对称离散后刚度矩阵对称既省存储又能用共轭梯度这类高效求解器。物理问题若本身不对称例如有对流项矩阵不对称但 Galerkin 处理流程不变这一点在第 5 章对拍时会再次遇到。2.3 本质边界条件与自然边界条件两类边界谁进右端项边界条件处理是弱形式最容易踩坑的地方先记住结论出现在分部积分边界项里的通量 c u′ 属于自然边界条件指定它就直接代入边界项解函数 u 本身的指定值属于本质边界条件必须从试解空间里强制扣除不能靠代入边界项实现。边界条件强形式写法弱形式处理典型物理量本质Dirichletu(x₀) ū试解空间压缩测试函数在边界取零固定位移、给定温度自然Neumannc u′(x₀) ḡ边界项替换为 ḡ·v(x₀)端部力、热流密度为什么 Dirichlet 边界不能“代入求解”因为弱形式是积分方程它只约束解的加权积分行为不约束点值u(x₀) ū 这种点约束在积分里被“抹平”了必须显式地从解的展开式中把对应自由度替换掉。实现上对应第 3 章的消去/置一法和第 4 章程序里的边界自由度处理。记住这条规则排错时能少一半的困惑。3. 弱形式的离散化基函数、刚度矩阵与全局组装3.1 有限维子空间与 Lagrange P1 基函数弱形式仍是无穷维问题v 可以是任意平方可积函数。要变成线性方程组必须把函数空间截断成有限维子空间最常见是分片线性多项式即 P1 单元。把 [0, L] 剖成 n 个单元节点坐标 x₁ … xₙ₊₁。每个节点 i 对应一个帽子函数 φᵢ(x)节点 i 上取值 1其余节点上取值 0单元内线性插值。这个“节点值为 1、其余为零”的插值性保证离散解在节点处的值就是未知向量 u 的分量。试解写成 u_h(x) Σᵢ uᵢφᵢ(x)未知量从“一个函数”变成“一个 Rⁿ⁺¹ 向量”。P1 基函数在每个单元上只有两个非零左端 N₁(ξ) 1 − ξ右端 N₂(ξ) ξξ (x − xₑ)/hₑ 是单元自然坐标hₑ 是单元长度。导数很干净N₁′ −1/hₑN₂′ 1/hₑ。“导数恰好是常数”这一点让单元刚度矩阵的积分变成一次算术。3.2 单元刚度矩阵一次积分算出 2×2 矩阵把 u_h 代入弱形式左侧、令 v 取基函数得到每个单元的贡献。单元 e 只有两个自由度单元刚度矩阵是 2×2kᵢⱼ⁽ᵉ⁾ ∫ₑ c(x)Nᵢ′(x)Nⱼ′(x) dx。当 c 在单元内为常数 cₑ 时代入两个常数导数就得到Kₑ (cₑ/hₑ) [[1, −1], [−1, 1]]这是有限元里最常见的矩阵。它的行和为零说明它只感受相对位移——刚体平移不产生应变能与物理直觉一致。右端等效节点载荷同样积分得到 Fᵢ⁽ᵉ⁾ ∫ₑ fNᵢ dxf 为常数 fₑ 时 Fₑ fₑhₑ/2 [1;1]相当于把单元总载荷平均分给两个节点。两个细节值得注意。其一c 在单元内变化梯度材料时不能提前面那个常数必须上数值积分其二f 是坐标函数时等效载荷的精度取决于积分规则第 4 章用一张表专门说明。还有一点容易被忽略这里积分是在物理坐标 x 上做的实际程序常在标准区间 [−1,1] 上做多出一个雅可比因子 hₑ/2后面随时要盯住它。3.3 全局组装局部自由度到全局自由度的映射单元矩阵算完剩下的就是“把局部矩阵放进全局矩阵”——这是组装唯一要做的事。一维情形映射规则很简单单元 e 的左节点对应全局自由度 e右节点对应全局自由度 e1相邻单元共享节点全局矩阵在该位置累加两个单元的贡献。单元编号 e局部自由度 1局部自由度 2全局索引数组112[1, 2]223[2, 3]eee1[e, e1]映射落在代码里就是索引数组 idx [e, e1]。MATLAB 组装循环如下它是第 4 章完整程序的核心K zeros(n1, n1); % 预分配全局刚度矩阵 F zeros(n1, 1); for e 1:n idx [e, e1]; % 单元 e 的两个全局自由度 he x(e1) - x(e); % 单元长度非均匀网格各自计算 Ke ce/he * [1, -1; -1, 1]; % ce 是单元材料参数 Fe fe * he/2 * [1; 1]; % fe 是单元载荷值 K(idx, idx) K(idx, idx) Ke; % 组装核心按索引累加 F(idx) F(idx) Fe; end参数说明idx 是单元到全局自由度的桥直接决定矩阵元素落位he 在非均匀网格下每个单元不同必须在循环内取值ce、fe 对常系数问题与单元无关但接口上按单元传值方便后续扩展为随空间变化的材料组装必须用累加而不是赋值因为内部节点同时属于两个单元。预分配零矩阵是必须的MATLAB 里逐次增长矩阵会反复分配内存节点数过万时明显拖慢程序。3.4 本质边界条件的处理消去法与罚函数法组装完的 K 是奇异的行和为零刚体位移还没被约束。本质边界条件施加后矩阵才可逆。最常用的是消去法把已知自由度从方程组里去掉只对未知自由度求解。实现上常见两种等价写法删除对应行列或把对应行列置零、对角线置 1、右端项置已知值。后者保持矩阵规模不变实现直观。另一种是罚函数法把本质边界条件当成大刚度弹簧加到对角元上例如 u(x₁) ū 时 K(1,1) 加一个大数 λ右端加 λū。优点是改动最小缺点是 λ 影响精度和条件数太大病态太小约束不牢。个人建议学习阶段用置一法工程代码用消去法罚函数法留给接触这类边界会动态变化的场景。提示置一法要求边界处理前 K 已经对称如果先处理边界再组装顺序颠倒会得到不对称矩阵排查时先确认组装循环完整跑完。4. MATLAB 有限元编程求解实例从弱形式到线性系统4.1 最小可运行程序一维 Poisson 问题把第 3 章的组装和第 2 章的弱形式串成完整程序。求解 −u″ 1x ∈ (0,1)边界条件 u(0) 0、u′(1) 0。解析解 u x − x²/2直接对拍。% 1D Poisson: -u 1, u(0) 0, u(1) 0 n 10; % 单元数 x linspace(0, 1, n1); % 节点坐标列向量 K zeros(n1, n1); F zeros(n1, 1); for e 1:n idx [e, e1]; % 单元自由度映射 he x(e1) - x(e); Ke (1/he) * [1, -1; -1, 1]; % c 1 Fe (he/2) * [1; 1]; % f 1 K(idx, idx) K(idx, idx) Ke; F(idx) F(idx) Fe; end % 本质边界条件 u(0)0置一法 K(1,:) 0; K(:,1) 0; K(1,1) 1; F(1) 0; u K \ F; % 单元内两点 Gauss 积分计算 L2 误差 xi [-1/sqrt(3), 1/sqrt(3)]; wt [1, 1]; L2 0; for e 1:n he x(e1) - x(e); idx [e, e1]; for q 1:2 xq x(e) (xi(q)1)/2 * he; % 积分点全局坐标 N [(1-xi(q))/2, (1xi(q))/2]; % 两点形状函数 uh N * u(idx); % 数值解在 xq 处的值 ue xq - xq^2/2; % 解析解 L2 L2 wt(q) * (uh - ue)^2 * (he/2); end end fprintf(L2 误差 %.3e\n, sqrt(L2));逻辑说明前半段与第 3 章组装完全一致边界处理用置一法把第 1 行第 1 列置零、对角线置 1、右端置 0等价于强制 u₁ 0处理后矩阵保持对称u K\F 是 MATLAB 线性求解入口对对称正定系统自动选 Cholesky。误差计算用单元内两点 Gaussxq 把标准积分点映射回物理坐标N*u(idx) 是有限元解在积分点的值he/2 是坐标变换的雅可比因子。这里有个容易误解的点对一维常系数问题线性单元的节点值有超收敛性质节点处误差接近机器精度所以不要用节点上的 max 范数判断精度一定要算单元内部的 L2 范数。上面的程序 n10 时 L2 误差约在 1e-3 量级加密到 n20 误差约降为四分之一。4.2 非均匀网格与数值积分Gauss 积分点数怎么设上面程序能跑通因为 c 和 f 在单元内是常数。真实问题的 c(x)、f(x) 是坐标函数单元刚度 ∫ cNᵢ′Nⱼ′ dx 和等效载荷 ∫ fNᵢ dx 都要数值积分。最常用是一维 Gauss-Legendre取若干积分点用加权和代替积分。规则是 N 点 Gauss 可精确积分 2N−1 次多项式。被积函数形式多项式最高阶次最少积分点数典型情形c 为常数刚度项01常截面杆单元f 为常数载荷项11均布载荷f 为线性函数22线性分布载荷c 为线性函数22变截面杆两点 Gauss 对三次多项式也精确因此工程惯例是刚度、载荷统一用两点代价几乎为零换来对任意光滑 f、c 的一致精度。用两点 Gauss 改写载荷计算的代码Fe [0; 0]; for q 1:2 xq x(e) (xi(q)1)/2 * he; % 标准点映射到物理坐标 N [(1-xi(q))/2, (1xi(q))/2]; % 单元上两点形状函数 Fe Fe wt(q) * f(xq) * N * (he/2); % 雅可比因子 he/2 end参数说明xi、wt 是标准区间 [−1,1] 的积分节点和权重xq 是物理坐标f 在 xq 处取值N 是形状函数行向量N′ 转置成列向量乘标量后得到 2×1 的贡献向量he/2 是 dx (he/2)dξ 的雅可比。自测方法把 f 取常数 1这段代码的结果应当与 Fe he/2[1;1] 完全一致。提示被积函数超过三次多项式时才需要三点以上的 Gauss。看到“积分点越多越准”的说法要分辨语境——增加积分点只影响求积精度不影响有限元解的收敛阶。4.3 调试技巧对称性检查与常见错误程序写完后先做三个快速检查再增加功能。第一检查刚度矩阵对称性。弱形式双线性形式对称离散后边界处理前的 K 必须对称。MATLAB 里运行 norm(K − K′,′fro′) 应为 0。不对称几乎都是组装索引写错典型是把累加写成赋值K(idx,idx) Ke 覆盖了相邻单元的贡献。第二配合何晓明老师公开课那套 MATLAB 有限元程序框架对照调试。其结构就是“组装、处理边界、求解”三步先跑通最小例程再逐步加功能。手写一遍的关键不是背代码而是记住每步对应的弱形式操作组装对应双线性形式 a(·,·)载荷对应线性泛函 F(·)边界处理对应本质边界条件的压缩。第三用精细网格对拍解析解。n 从 10 加密到 40L2 误差应大致下降一个量级线性单元收敛阶为 2。误差不降问题几乎总在边界条件处理最常见的是把自然边界条件又手动设了一次值导致约束重复。5. 弱形式程序的验证技巧制造解与收敛率对拍5.1 制造解方法的基本流程解析解只存在于教科书问题。制造解方法MMS的思路是反着用弱形式先任意选一个光滑函数 u*(x) 当“真解”代回强形式反解出载荷 f(x)再带着这个 f 去跑程序比较数值解与 u*。u* 不需要有物理意义只要足够光滑、边界形式与你的问题一致。对 −(cu′)′ f制造解载荷是 f −(cu*)″正好用符号计算避免手推错误。常见做法是选三角函数和指数函数的组合让解不是低阶多项式从而真正检验积分和组装。5.2 变系数问题的符号生成与对拍取变系数 c(x) 1 x制造解 u* sin(πx)边界取两端 Dirichlet u(0) u(1) 0。生成载荷的代码import sympy as sp x sp.symbols(x) c 1 x ustar sp.sin(sp.pi * x) f sp.simplify(-sp.diff(c * sp.diff(ustar, x), x)) print(f) # pi**2*(x 1)*sin(pi*x) - pi*cos(pi*x)sp.diff 两次求导对应强形式左端得到的符号表达式可以直接写成 MATLAB 函数句柄f (x) pi^2*(x1).sin(pix) − picos(pix)然后跑第 4 章程序。注意此时 c 是坐标的函数刚度矩阵也要用两点 Gauss不能用常数公式。这个对拍同时验证了变系数刚度项和载荷项两处积分。5.3 收敛率对拍用误差比率代替肉眼看曲线量化验证的最后一环是收敛率网格每次加密一倍计算 L2 误差 e ‖u_h − u*‖线性单元理论收敛阶为 2误差应近似缩小 4 倍。列出几个网格的误差和相邻比率比看一眼对数坐标曲线的斜率更严格。单元数 nL2 误差误差比率10≈5.2e-3—20≈1.3e-34.040≈3.2e-44.080≈8.0e-54.0比率稳定在 4 附近说明弱形式程序实现正确比率接近 2 时先查积分点数再查本质边界是否被重复约束比率不下降时从组装索引查起。实际排错中按“边界约束 → 积分规则 → 组装索引”这个顺序排查是我遇到过的最高效路径。本文还有配套的精品资源点击获取