企业官网建设流程全解析

简介：资源内容为一维伽辽金型无网格法的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_I

s ≥ 1 时取 0。它满足 w(0)=1、w(1)=0、w'(0)=w'(1)=0，边界处函数值和导数都连续，对线性基来说够光滑，代码也短。追求更高光滑性可以换四次样条 w(s) = 1 - 6s² + 8s³ - 3s⁴，但导数表达式多两项，前期调试没必要。如果求解三阶以上方程，建议换高斯型权函数并适当缩小影响域，否则形函数光滑度不够会直接影响收敛阶。

影响域半径 r_I = α·h，h 是节点间距。α 取多少，是无网格法里最像“玄学”的一个参数。一维线性基下我的经验值：

对均匀一维网格，α=2.4 时绝大多数 Gauss 点覆盖 3 个节点，个别靠近边界的点覆盖 2 个，刚好满足线性基最少节点数，又不至于让刚度矩阵带宽过大。跑收敛性测试前先固定 α=2.4，程序通了再在 2.0~2.8 之间扫一遍，看误差曲线找最优。α 增大并不总是更好，覆盖节点多了，形函数彼此重叠加深，精度会被过度平均拖下来。

2.3 从弱形式到刚度矩阵：一维方程的组装路径

用最经典的一维模型问题演示：-u'' + u = f，u(0)=u(1)=0。乘检验函数 v(x) 并在 (0,1) 上积分，分部积分后得到弱形式：

∫(u'v' + uv) dx = ∫ f v dx

把 u^h = Σφ_I u_I、v = φ_J 代入，得到线性方程组 K u = F，其中：

K_IJ = ∫(φ_I' φ_J' + φ_I φ_J) dx，F_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）主程序 % 求解 -u''+u = f, x in (0,1), u(0)=u(1)=0 % f = (pi^2+1)*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(e+1); Jc = 0.5 * (xb - xa); % 坐标变换雅可比 for q = 1 : 3 xq = 0.5*(xa+xb) + 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.'); % 边界值都是 0，F 不需要加修正项 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 分别代表 x=0 和 x=1 两点的形函数向量。由于 α 在 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| / r，r = 影响域半径，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 end

dwdx 在 s=0 处有符号不定问题，但 dwds 在 s=0 处正好是 0，所以 dx=0 时 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_I，h_I 取该节点与最近邻居的距离，避免稀节点区域覆盖节点数不足。

5.2 现象：解出来是锯齿波

现象：解曲线在节点之间振荡，看起来像有限元没做稳定化。

原因：背景积分点数太少。MLS 形函数每个 Gauss 点上重新构造，被积函数在背景单元内不是多项式，2 点 Gauss 的分辨率不足以覆盖影响域重叠区的变化，刚度矩阵欠积分，出现额外零空间模式。

解决：每个背景单元从 2 点 Gauss 提到 3 点或 4 点，同时把 α 调到 2.4 以上。锯齿依然在的话，把 Gauss 点数和 α 一起检查：一个管积分分辨率，一个管形函数重叠宽度，两者合起来决定有效积分采样密度。

5.3 现象：边界处结果明显偏离解析解

现象：内部解很好，x=0 和 x=1 附近误差突然变大，边界重构值 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 点对全部节点搜索支持域，N=1000 时程序跑几分钟，刚度矩阵满阵存储，内存直接爆掉。

原因：无网格法的形函数支持域重叠，在二维下每个点覆盖的节点数随 α 平方增长，且没有单元拓扑可以提前限定影响域。一旦每个积分点都对全部节点计算距离，总计算量是 N × GaussPts × N 量级。

解决：二维版本必须做空间搜索，常见做法是分块网格或 k-d 树，先筛出候选节点，再在候选节点上组装 A、B；刚度矩阵按带状或稀疏存储，因为每个点的影响域宽度约 2α·h，矩阵实际带宽有限。一维程序练熟后，二维扩展最先改这三处：节点搜索、稀疏组装、背景网格积分。

6. 收敛阶验证与二维EFG扩展的入场准备

程序写完后最重要的一件事是验证收敛阶，这比看任何文档都诚实。误差数字不会骗人，L2 阶低于 1.5 的程序一定还有潜伏问题。三个必查指标：

跑收敛测试时不要只测一组 N。把 N 从 11 跳到 81，每次翻倍，误差按 log-log 斜率看。如果 α 固定在 2.4，L2 阶稳定在 2.0 附近，程序基本可信。再用 α=2.0 和 2.8 各跑一遍，观察误差曲线变化，能帮你对这个参数建立感觉。

往二维扩展之前，先确认一维程序的三个习惯已经刻进肌肉：第一，所有评价点都走 MLS1D，不直接使用自由度当节点值；第二，本质边界用罚函数或拉格朗日乘子，不直接改矩阵行；第三，背景积分点数宁可多不要少。二维 EFG 的形函数推导和组装逻辑一维完全一致，只是 A 矩阵变 3×3 或 6×6，B 变宽，节点搜索必须用空间数据结构。

我现在拿到任何无网格法程序，第一件事永远是跑一遍收敛阶表格，而不是先看云图。误差数字不会骗人，L2 阶低于 1.5 的程序一定还有潜伏问题。这套一维程序你能跑出二阶收敛，再往二维三维走才有底气。希望帮到你。

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