简介资源内容为一维伽辽金型无网格法的MATLAB编程实现面向正在学习无网格法、伽辽金法或需要编写数值计算代码的高校学生与科研人员。程序主体是一个可运行的m脚本文件配合一个rar压缩包共两个文件整体仅8KB非常轻量便于快速下载使用。通过这份代码读者可以完整理解一维问题中无网格法的主要流程节点影响域划分、形函数构造、刚度矩阵与荷载向量组装、位移边界条件处理等关键环节并能对照理论推导逐行检查具体实现细节。代码结构清晰模块划分明确适合作为课程实践、本科毕业设计或科研起步的参考模板同时在一维示例基础上读者也可以进一步扩展为二维或非线性问题。目前已有646人学习下载这一热度说明其对于快速上手无网格法编程具有较好的实用参考价值。1. 一维伽辽金无网格法MATLAB程序比有限元多算三步却能绕开网格有限元的前处理大头在网格划分一维伽辽金型无网格法MATLAB程序把这块网格依赖整个去掉节点照布背景单元照积分但形函数改用移动最小二乘在每个积分点上现算。省掉了网格生成代价是多算三步——每个 Gauss 点都要重新搜索支持节点、求矩数组逆、重构形函数导数还要用罚函数或拉格朗日乘子去“补”本质边界条件。这篇笔记围绕这套一维程序讲清楚MLS 形函数怎么算、刚度矩阵怎么组装、影响域半径和积分点数怎么选、边界条件和后处理有哪些坑最终给出一份能在 MATLAB 里直接跑通并验证收敛阶的最小实现。适合看完有限元想转无网格法的新手也适合准备扩展二维前后先回来打基础的人。2. 伽辽金型无网格法的离散原理为什么形函数不插值还要用2.1 MLS形函数用局部最小二乘拼出全域近似伽辽金型无网格法最常对应的叫法是无单元伽辽金法也就是 EFG。它和有限元共用同一套弱形式差别全在形函数上。一维程序里最基础也最关键的就是移动最小二乘形函数。把未知场写成叠加形式u^h(x) Σ φ_I(x) · u_I对固定求值点 x先在影响域内找到 n 个节点比如编号 i₁ 到 iₙ。在这一点附近用一个线性基 p(y) [1, y]ᵀ 去逼近真实函数逼近系数 a(x) 由局部加权最小二乘确定J(a) Σ w(x - xᵢⱼ) · [pᵀ(xᵢⱼ)a - uᵢⱼ]²对 a 求极小得到 a A⁻¹ B u于是φ(x) pᵀ(x) · A⁻¹(x) · B(x)其中 A Σ wⱼ p(xⱼ)pᵀ(xⱼ)B [w₁p(x₁), …, wₙp(xₙ)]。注意 A、B 里的 p 用的是节点坐标 xⱼ不是求值点 x。A 是 2×2 矩阵B 是 2×n 矩阵。这里最反直觉的一点φᵢ(xⱼ) 不等于 δᵢⱼ。你把第 J 个自由度取 1其余取 0得到的函数在 x_J 处并不等于 1这就是“形函数不插值”。它带来两个直接后果。第一边界上的 u(0)0 不能写成一维自由度 u₁0因为边界点的真实数值由整个自由度向量共同决定。第二后处理看位移时不能把 u_I 当节点位移直接画曲线。这两个问题会在后面专门处理。为什么非要用 MLS 而不是多项式插值因为 MLS 天然能随着求值点位置改变参与拟合的节点集合影响域外节点自动被权函数筛掉形函数光滑性由权函数保证只要 A 可逆形函数总是定义的。代价就是每个求值点都要动态组装一次 A、B计算量比有限元形函数大一个量级。2.2 权函数与影响域半径α取2.0还是3.0结果差一个量级权函数决定每个节点对拟合的贡献强度。一维下我常用三次样条w(s) 1 - 3s² 2s³s |x - x_I| / r_Is ≥ 1 时取 0。它满足 w(0)1、w(1)0、w(0)w(1)0边界处函数值和导数都连续对线性基来说够光滑代码也短。追求更高光滑性可以换四次样条 w(s) 1 - 6s² 8s³ - 3s⁴但导数表达式多两项前期调试没必要。如果求解三阶以上方程建议换高斯型权函数并适当缩小影响域否则形函数光滑度不够会直接影响收敛阶。影响域半径 r_I α·hh 是节点间距。α 取多少是无网格法里最像“玄学”的一个参数。一维线性基下我的经验值α 取值每个积分点覆盖节点数表现 1.5可能只有 1 个A 矩阵奇异程序直接报错1.8 ~ 2.52 ~ 3 个精度高带宽小首选区间3.0 左右4 ~ 5 个形函数变宽精度略降带宽变大 4.06 个以上结果过于光滑局部特征被抹平对均匀一维网格α2.4 时绝大多数 Gauss 点覆盖 3 个节点个别靠近边界的点覆盖 2 个刚好满足线性基最少节点数又不至于让刚度矩阵带宽过大。跑收敛性测试前先固定 α2.4程序通了再在 2.0~2.8 之间扫一遍看误差曲线找最优。α 增大并不总是更好覆盖节点多了形函数彼此重叠加深精度会被过度平均拖下来。2.3 从弱形式到刚度矩阵一维方程的组装路径用最经典的一维模型问题演示-u u fu(0)u(1)0。乘检验函数 v(x) 并在 (0,1) 上积分分部积分后得到弱形式∫(uv uv) dx ∫ f v dx把 u^h Σφ_I u_I、v φ_J 代入得到线性方程组 K u F其中K_IJ ∫(φ_I φ_J φ_I φ_J) dxF_I ∫ f φ_I dx这个方程组形式和有限元一模一样差别只在 φ_I 每个积分点上都要重新动态构造。一维实现里不需要单元形函数库流程变成三层循环外层遍历背景单元中层遍历 Gauss 点内层调用 MLS 形函数函数算完 φ 和 φ 后累加 K 和 F。背景单元只是积分载体和自由度没有任何拓扑绑定。这个结构也是后续二维 EFG 程序的基本骨架。3. 一维程序的模块拆解节点、背景积分和边界条件各自怎么落3.1 节点离散与背景积分单元网格不在明面仍在后台无网格法是“无单元”不是“无积分”。被积函数是形函数及其导数必须落在某个积分结构上。常见做法是直接用节点坐标把求解区间切成 N-1 段背景单元每段上配 Gauss 点。背景单元与有限元网格的区别在于它不承载形函数定义只是积分载体可以任意加密而不影响自由度数目。节点怎么布一维最简单是 linspace(0,1,N)。无网格法对非均匀节点的容忍度比有限元高因为形函数是局部重拟合而不是固定单元但节点间距突变过大时建议每个节点用自己的局部间距 h_I 定半径即 r_I α·h_I否则稀节点处覆盖过少节点。本程序用均匀节点r_i 统一。背景单元上 Gauss 点数要单独说。每个单元 2 点 Gauss 在均匀节点、线性基下勉强能算但 MLS 形函数在单元内变化不像有限元多项式那样平滑权函数在节点附近快速衰减积分分辨率不够时容易欠积分。实际使用 3 点或 4 点本程序用 3 点对边界层等剧烈变化问题可加到 5 点。每增加一个 Gauss 点所有积分点的支持节点搜索和 2×2 矩阵求逆都要重来一遍计算量线性上涨所以不是越高越好够用即可。3.2 本质边界条件的处理罚函数法与拉格朗日乘子法怎么选因为形函数不插值本质边界条件不能像有限元那样把某一行改成单位向量。边界节点自由度 u_b 不是边界点位移 u^h(x_b)直接置零相当于给了一个错误约束解出来边界一塌糊涂。罚函数法的思路是在弱形式里加边界罚项。一维情况下对每个边界点 x_b给刚度矩阵加一项 β·φ(x_b)·φᵀ(x_b)给载荷向量加 β·g_b·φ(x_b)。离散形式K β φ(x_b) φᵀ(x_b) F β g_b φ(x_b)β 典型值取 1e6 ~ 1e8。β 太小边界约束没压实边界残差大β 太大导致刚度矩阵条件数恶化 10 次方量级甚至把主对角元淹没。罚函数法是有“后悔药”的——β 取错了解出来一眼能看出来边界残差明显偏大调大 β 重算一遍就好。拉格朗日乘子法把约束写进系统[[K, Cᵀ], [C, 0]] · [[u], [λ]] [[F], [g]]其中 C_IJ φ_I(x_bJ)。优点是约束被精确满足、不引入病态缺点是矩阵变成不定矩阵自由度混入 λ求解器不能随便用常规消元。一维入门建议先用罚函数跑通后想追求边界精度再换拉格朗日。对于一维问题罚函数造成的边界误差在工程范围内可接受只要 β 足够大。3.3 后处理重构u_I不是节点位移直接画图会翻车把无网格法自由度当有限元节点位移用是新手最常见的事故。解出 u 数组后直接 plot(xNode, u)结果曲线和解析解“有点像但又不对”边界尤其怪。这不是程序有 bug是 u_I 本来就不是位移值。要看真实位移必须在评价点集合上重新组装形函数u^h(x_eval) Σ φ_I(x_eval) · u_I代码层面就是对每个绘图点调用一次 MLS 形函数函数得到形函数向量后点乘自由度向量。绘图点取 200~500 个不要在节点上画。应力应变恢复同理只不过要多算一个导数形函数 dφ。这个习惯要带到二维云图显示前必须把自由度映射到场值不能用节点色阶直接画。4. 用MATLAB跑通一维伽辽金无网格法可复现代码与参数说明4.1 主程序 efg1d.m组装、求解、后处理一条龙程序按三个函数文件加一个测试脚本组织efg1d.m 是主程序MLS1D.m 负责形函数wfun3.m 是权函数test_convergence.m 做收敛性验证。主程序如下function [u, xNode] efg1d(N, alpha, beta) % 一维伽辽金型无网格法EFG主程序 % 求解 -uu f, x in (0,1), u(0)u(1)0 % f (pi^21)*sin(pi*x)解析解 u sin(pi*x) % 输入N 节点数alpha 影响域半径系数beta 罚因子 % 输出u 广义自由度注意不是节点位移xNode 节点坐标 xNode linspace(0, 1, N); h xNode(2) - xNode(1); r_i alpha * h * ones(N, 1); % 每个节点的影响域半径 % 三点 Gauss 积分作用在 [-1,1] gx [-sqrt(3/5), 0, sqrt(3/5)]; gw [5/9, 8/9, 5/9]; K zeros(N, N); F zeros(N, 1); f (x) (pi^2 1) * sin(pi * x); for e 1 : N-1 xa xNode(e); xb xNode(e1); Jc 0.5 * (xb - xa); % 坐标变换雅可比 for q 1 : 3 xq 0.5*(xaxb) Jc*gx(q); wq Jc * gw(q); [phi, dphi] MLS1D(xq, xNode, r_i); K K wq * (dphi*dphi. phi*phi.); F F wq * (f(xq) * phi); end end % 罚函数法施加本质边界条件 u(0)u(1)0 [phi0, ~] MLS1D(xNode(1), xNode, r_i); [phi1, ~] MLS1D(xNode(N), xNode, r_i); K K beta * (phi0*phi0. phi1*phi1.); % 边界值都是 0F 不需要加修正项 u K \ F; end三个 Gauss 点对应的 gx、gw 可以直接背下来不需要额外调 lgwt 工具函数。组装刚度矩阵时写成 dphidphi. 和 phiphi.因为 dphi、phi 都是 N 维列向量外积得到 N×N 矩阵正好对应 φ_Iφ_J 和 φ_Iφ_J 对所有 I、J 的组合。MLS1D 返回的 phi、dphi 是与节点等长的列向量不在支持域内的位置自动置 0所以累加时不需要做编号映射这是和一维有限元最大的编程区别。罚函数边界修改里phi0 和 phi1 分别代表 x0 和 x1 两点的形函数向量。由于 α 在 2.0~2.5 时两边支持域不会互相重叠两块修正互不干扰。边界值都是 0所以 F 不需要动。如果你换个边界值非零的问题一定记得在 F 里加上 β·g·φ(x_b)。4.2 MLS1D 与权函数形函数和导数在哪算、为什么这么写MLS 形函数是这套程序的核心把它单独拎出来看function [phi, dphi] MLS1D(xq, xNode, r_i) % 一维线性基 MLS 形函数及导数的完整实现 % xq 求值点标量xNode 节点坐标列向量r_i 影响域半径列向量 tol 1e-12; idx find(abs(xq - xNode) r_i tol); % 影响域内节点 n length(idx); if n 2 error(xq%.3f 处支持节点数不足 2增大 alpha 或加密节点, xq); end A zeros(2, 2); dA zeros(2, 2); B zeros(2, n); dB zeros(2, n); for j 1 : n xj xNode(idx(j)); rj r_i(idx(j)); s abs(xq - xj) / rj; [w, dwdx] wfun3(s, rj, xq - xj); pj [1; xj]; % 在节点 xj 处求值的基向量 A A w * (pj * pj); dA dA dwdx * (pj * pj); B(:, j) w * pj; dB(:, j) dwdx * pj; end Ai inv(A); % 2x2 矩阵逆 dAi -Ai * dA * Ai; % (A^{-1}) pe [1; xq]; dpe [0; 1]; phi (pe * Ai * B); % 1xn dphi dpe * Ai * B ... pe * dAi * B ... pe * Ai * dB; % 1xn % 扩展回 N 维与节点编号对齐 phiOut zeros(size(xNode)); dphiOut zeros(size(xNode)); phiOut(idx) phi; dphiOut(idx) dphi; end两个地方特别容易写错。第一A、B 是用节点坐标 xj 算的不是用求值点 xq 算的第二dA 不能漏因为 A 里含权函数 w而 w 随求值点 x 变化A 对 x 的导数必须通过权函数导数 dwdx 传入。很多人写无网格法翻车就是只对 B 求导、漏了 dA 这一项导致导数形函数完全错误收敛阶直接崩掉。权函数子程序function [w, dwdx] wfun3(s, r, dx) % 三次样条权函数及对 x 的导数 % s |xq - xj| / rr 影响域半径dx xq - xj if s 1 w 0; dwdx 0; else w 1 - 3*s^2 2*s^3; % w(0)1, w(1)0 dwds -6*s 6*s^2; % dw/ds dwdx dwds * sign(dx) / r; % 链式求导 end enddwdx 在 s0 处有符号不定问题但 dwds 在 s0 处正好是 0所以 dx0 时 sign(0)0整个表达式为 0不存在 NaN 风险。这个权函数在 s≥1 的返回值是 0而主程序传入的求值点都在影响域内兜住边界容差引起的异常情况。4.3 收敛性测试换N、换alpha确认程序没写错程序写完先别急着算实际问题跑一遍收敛性测试% test_convergence.m alphas [2.0, 2.4, 2.8]; Ns [11, 21, 41, 81]; for a alphas fprintf(alpha %.1f\n, a); fprintf(%4s %12s %12s %10s\n, N, L2err, H1err, L2阶); errL2_prev inf; for N Ns beta 1e7; [u, xNode] efg1d(N, a, beta); h xNode(2) - xNode(1); xE linspace(0, 1, 200); uH zeros(200, 1); duH zeros(200, 1); for k 1 : 200 [phi, dphi] MLS1D(xE(k), xNode, a*h*ones(N,1)); uH(k) phi * u; duH(k) dphi * u; end uex sin(pi*xE); duex pi * cos(pi*xE); errL2 sqrt(trapz(xE, (uH-uex).^2)); errH1 sqrt(trapz(xE, (duH-duex).^2) errL2^2); order log(errL2_prev / errL2) / log(2); fprintf(%4d %12.3e %12.3e %10.2f\n, N, errL2, errH1, order); errL2_prev errL2; end end收敛阶按节点数翻倍、误差减半的方式评估。一维线性基 MLS 在积分充分时L2 误差期望接近二阶H1 误差接近一阶。如果你跑出来 L2 阶只有 0.9优先查 dA 项再查是不是罚因子太小。这套程序跑通后就可以开始改边界条件、换权函数、换基函数了。5. 避坑一维伽辽金无网格法最常见的五个翻车点这一章列的问题按踩中频率排序绝大多数不是数学问题是实现细节。如果你跑出来的结果不收敛不要先怀疑理论按清单排查。5.1 现象刚度矩阵条件数爆炸现象K\F 直接得到 NaN或者 cond(K) 超过 1e16。原因某个 Gauss 点处 MLS 支持域内有效节点数少于基函数个数。线性基最少需要 2 个节点但实际要求 A 矩阵良态最好 3 个以上。常见于 α 取得太小或者非均匀节点局部间距过大。解决在 MLS1D 里先数 idx 长度小于 2 就报错提示增大 α。均匀节点从 α2.4 起步非均匀节点改成每个节点独立半径 r_I α·h_Ih_I 取该节点与最近邻居的距离避免稀节点区域覆盖节点数不足。5.2 现象解出来是锯齿波现象解曲线在节点之间振荡看起来像有限元没做稳定化。原因背景积分点数太少。MLS 形函数每个 Gauss 点上重新构造被积函数在背景单元内不是多项式2 点 Gauss 的分辨率不足以覆盖影响域重叠区的变化刚度矩阵欠积分出现额外零空间模式。解决每个背景单元从 2 点 Gauss 提到 3 点或 4 点同时把 α 调到 2.4 以上。锯齿依然在的话把 Gauss 点数和 α 一起检查一个管积分分辨率一个管形函数重叠宽度两者合起来决定有效积分采样密度。5.3 现象边界处结果明显偏离解析解现象内部解很好x0 和 x1 附近误差突然变大边界重构值 u^h(x_b) 与给定边界值差得远。原因罚函数法的罚因子太小边界约束没压实或者你直接给边界自由度赋值了。因为 φ_I(x_J)≠δ_IJ直接赋 u(0)0 等于是给节点自由度一个错误的约束。解决确认施加的是罚函数形式而不是直接赋值。β 先用 1e7然后看边界重构残差 u^h(x_b) - g。残差小于 1e-6 说明 β 够用否则继续推高 β。想彻底摆脱 β 选择换拉格朗日乘子法。5.4 现象把u_I直接当节点位移画图现象后处理曲线在节点处“对不上”以为程序算错了其实画法错了。原因u_I 是广义自由度不是 u 在 x_I 处的值。MLS 形函数不插值u^h(x_I)≠u_I。直接 plot(xNode, u) 画的是系数向量不是位移场。解决统一用评价点后处理。生成 xE linspace(0,1,200)对每个点算 φ * u再画图。这个习惯同样要带到二维云图显示前必须把自由度映射到场值不能用节点色阶直接画。5.5 现象二维程序一写就卡死或带宽失控现象把一维逻辑直接搬到二维每个 Gauss 点对全部节点搜索支持域N1000 时程序跑几分钟刚度矩阵满阵存储内存直接爆掉。原因无网格法的形函数支持域重叠在二维下每个点覆盖的节点数随 α 平方增长且没有单元拓扑可以提前限定影响域。一旦每个积分点都对全部节点计算距离总计算量是 N × GaussPts × N 量级。解决二维版本必须做空间搜索常见做法是分块网格或 k-d 树先筛出候选节点再在候选节点上组装 A、B刚度矩阵按带状或稀疏存储因为每个点的影响域宽度约 2α·h矩阵实际带宽有限。一维程序练熟后二维扩展最先改这三处节点搜索、稀疏组装、背景网格积分。6. 收敛阶验证与二维EFG扩展的入场准备程序写完后最重要的一件事是验证收敛阶这比看任何文档都诚实。误差数字不会骗人L2 阶低于 1.5 的程序一定还有潜伏问题。三个必查指标检查项预期值异常时优先排查L2 误差收敛阶≈ 2.0dA 项、Gauss 点数、罚因子 βH1 误差收敛阶≈ 1.0形函数导数、权函数连续性边界重构残差 1e-6β 太小或边界装配写错跑收敛测试时不要只测一组 N。把 N 从 11 跳到 81每次翻倍误差按 log-log 斜率看。如果 α 固定在 2.4L2 阶稳定在 2.0 附近程序基本可信。再用 α2.0 和 2.8 各跑一遍观察误差曲线变化能帮你对这个参数建立感觉。往二维扩展之前先确认一维程序的三个习惯已经刻进肌肉第一所有评价点都走 MLS1D不直接使用自由度当节点值第二本质边界用罚函数或拉格朗日乘子不直接改矩阵行第三背景积分点数宁可多不要少。二维 EFG 的形函数推导和组装逻辑一维完全一致只是 A 矩阵变 3×3 或 6×6B 变宽节点搜索必须用空间数据结构。我现在拿到任何无网格法程序第一件事永远是跑一遍收敛阶表格而不是先看云图。误差数字不会骗人L2 阶低于 1.5 的程序一定还有潜伏问题。这套一维程序你能跑出二阶收敛再往二维三维走才有底气。希望帮到你。本文还有配套的精品资源点击获取