简介这是一份面向无网格法与计算力学初学者的MATLAB源码资源对应一维伽辽金型无网格法的基础编程实现可用于学习从形函数构造到方程求解的完整流程。压缩包内共2个文件其中m文件为主要程序涉及移动最小二乘形函数、伽辽金弱形式离散、刚度矩阵与荷载向量组装、边界条件处理等关键模块另一个rar文件为配套的算例或辅助文件整套资源仅8KB体量轻、结构紧凑便于逐行阅读和调试。借助该代码读者可以快速掌握一维无网格问题的编程框架理解基函数、权函数、节点影响域等概念对结果的影响并以此为基础扩展到二维或非线性问题。目前已有646人学习浏览适合本科生、研究生及MATLAB编程爱好者用于数值方法课程设计、科研入门或相关课题二次开发。1. 一维伽辽金型无网格法一个MATLAB程序包能帮你少走多少弯路先说结论无网格法不是不用网格而是不用“单元”这个概念。你手里这份资源是教材第三章配套的一维伽辽金型无网格法MATLAB程序解压出来就是一套能跑的、从节点生成到误差输出全流程的脚本。它解决的核心问题很直接当你不想碰有限元里那套网格剖分、想试试只用散乱节点做数值逼近时这个程序就是你从头复现无网格法的最短路径。适合谁刚接触无网格法的研究生、打算用无网格法写课程设计或小论文的初学者以及那些已经看过公式但卡在“代码该从哪写起”的从业者。你不需要先精通MATLAB能把脚本跑通、能改参数就已经比只看公式强很多。2. 一维伽辽金型无网格法的核心机制有限元形函数与MLS形函数的本质差异2.1 为什么从一维起步MLS形函数到底在算什么无网格法的种类很多伽辽金型Element-Free Galerkin MethodEFG是其中最成熟的一支它以移动最小二乘Moving Least Squares, MLS形函数作为近似函数。你可能会问既然最终都要解一个类似有限元的线性方程组那它跟有限元到底差在哪答案是形函数。有限元形函数是在单元内定义的形函数只依赖单元内的几何形状而且形函数在节点上满足克罗内克δ性质第I个形函数在节点I处取值为1在其他节点处取值为0。这意味着有限元的未知量就是节点处的真实位移值边界条件可以直接“硬塞”进方程。MLS形函数则完全不同。它基于一组离散节点在每个求解点x处做一次加权最小二乘拟合得到该点的形函数值。换个说法MLS形函数是“看周围若干节点脸色”的它在节点处取值既不是1也不是0是某个中间数值。所以MLS形函数不满足克罗内克δ性质这是无网格法一切边界处理麻烦的根源。下面这张表帮你看清区别。特征有限元形函数MLS无网格形函数定义域局部单元影响域内所有节点节点处取值满足克罗内克δ不满足连续性控制受单元阶次限制由权函数直接决定本质边界条件施加直接修改方程罚函数/Lagrange乘子网格依赖强依赖完全不需要网格由于MLS形函数在任意点都至少有C0甚至C1连续性它不需要像有限元那样做节点协调这是很多人愿意学它的原因。但代价就是求函数值时必须做矩阵求逆计算成本比有限元形函数高了不少。2.2 权函数与影响域形函数计算的Matlab实现MLS形函数的计算链路是给定求解点x的一组基函数和周围节点先组装矩方程再解出形函数。一维情形下基函数p(x)通常取线性基或者二次基权函数则采用三次样条权函数。我一般建议直接写一个独立的形函数函数后续改参数不需要翻主程序。function [phi, dphi] mls_shape(x, nodes, h, p_order, dmax) % x : 求解点坐标 % nodes : 所有节点坐标向量 % h : 节点平均间距 % p_order: 基函数阶次, 1线性基 [1,x], 2二次基 [1,x,x^2] % dmax : 影响域无量纲半径, 实际影响域半径 R dmax * h R dmax * h; n length(nodes); % --- 选取影响域内的节点 --- idx find(abs(nodes - x) R); m length(idx); if m p_order 1 error(影响域内节点数不足, 请调大 dmax); end % --- 组装基函数矩阵和权函数对角阵 --- P zeros(m, p_order1); W zeros(m, m); for j 1:m dx x - nodes(idx(j)); P(j,:) dx.^(0:p_order); % 基函数在局部坐标下计算 s abs(dx) / R; w 1 - 3*s^2 2*s^3; % 三次样条权函数 W(j,j) w; end % --- 计算 A 矩阵和 B 向量 --- A P * W * P; pt (x - nodes(idx)).^0; % 当前点的基函数向量 B P * W; invA inv(A); phi pt * invA * B; % 形函数行向量 % --- 导数用五点差分近似, 教学代码最省事的做法 --- eps1 1e-6; [phi_plus, ~] mls_shape(xeps1, nodes, h, p_order, dmax); [phi_minus, ~] mls_shape(x-eps1, nodes, h, p_order, dmax); dphi (phi_plus - phi_minus) / (2*eps1); end这段代码的逻辑我拆成三步讲。第一步是选节点以求解点x为圆心、半径Rdmax*h画一个圆圈进来的节点才参与计算。第二步是组装A矩阵A矩阵是一个(p_order1)维方阵它本质是加权最小二乘的法方程矩阵如果影响域内节点太少A矩阵会奇异这就是程序报错“矩阵接近奇异”的最常见原因。第三步是求形函数phi最终是一个行向量长度等于影响域内节点数对应每个节点对当前点x的“贡献权重”。别忘了权函数在s1处恰好衰减为0所以影响域边界上的节点权重是零。这也是一个隐患如果高斯点恰好落在两个节点影响域的交界上权函数可能出现数值零导致A矩阵不可逆。后面避坑章节我会专门讲。2.3 Galerkin弱式与整体刚度矩阵的组装流程形函数只是工具真正要解的是一维扩散方程-d/dx(k·du/dx) f。这里的k是材料参数f是源项。标准的伽辽金做法是在方程两边乘以测试函数然后分部积分得到弱形式。弱形式的核心变化是原来方程里有二阶导数分部积分后变成一阶导数乘积的积分加上边界项。边界上自然边界条件如力边界会直接进入边界积分项这跟有限元完全一致。对一维区间[0,L]刚度矩阵的表达式为K_IJ ∫₀ᴸ k(x)·φ₁(x)·φ_J(x)dx载荷向量为F_I ∫₀ᴸ f(x)·φ_I(x)dx再加上边界产生的修正项。矩阵/向量维度物理含义K (刚度矩阵)N×N节点之间的刚度耦合关系F (载荷向量)N×1等效节点载荷A (矩方程)(p1)×(p1)每个求解点处的加权最小二乘法方程Φ (形函数)N×1 (对单点)各节点对当前点的权重这里的K矩阵不再是有限元那种“带状”结构因为每个高斯点都跟影响域内所有节点相连节点的关联范围变大带宽由dmax决定。这也是无网格法比有限元计算量大的原因之一。把dmax从2.0调到3.5矩阵的非零元素数量会明显增加求解时间跟着涨。3. 主程序逐段拆解从解压到跑出第一条曲线3.1 解压后的文件结构与运行顺序这份资源下下来是zip包解压后你会看到几个.m脚本。我建议先不要急着在MATLAB里双击运行先把目录结构看清楚。常规组织方式是一个主程序加若干函数文件主程序负责定义模型参数、节点生成、调用组装函数、求解并画图形函数函数负责计算单点形函数还有误差分析脚本负责输出收敛阶。用命令行进入目录会比较稳避免MATLAB当前路径不对导致函数找不到。# 在任意终端里解压并查看文件结构 cd ~/Downloads unzip 一维伽辽金型无网格法MATLAB程序.zip -d efg_1d ls -l efg_1d/*.m如果你在Windows上解压直接右键“解压到当前文件夹”也行唯一需要注意的是中文文件名和中文注释的编码问题这个放到第5章讲。进入MATLAB后把当前文件夹切换到解压目录然后直接运行主程序。如果你的MATLAB版本偏低比如R2016a大概率也能跑因为程序里没用到太新的语法如果用的是2023a及以上版本反而要留意中文注释乱码问题低版本时反而没有这个毛病。这里提醒一句MATLAB在线网页版mathworks.com的MATLAB Online也能跑这类脚本但它对文件读写路径的管理比较特殊你要先把整个zip在本地解压后逐个上传比较麻烦不如本地装一个桌面版省心。3.2 刚度矩阵组装与边界条件处理核心代码主程序核心段长这样。我按工程实战习惯把“模型参数定义 - 节点生成 - 高斯积分组装 - 边界处理 - 求解画图”分成几块彼此用注释隔开方便你改参数时不会滚到乱七八糟的代码里。% 1. 模型参数 L 1.0; % 求解域长度 N 21; % 节点数 h L / (N - 1); % 节点间距 p_order 1; % 基函数阶次 dmax 2.5; % 影响域半径系数 alpha 1e8; % 罚函数系数 n_gauss 2; % 每个区间的高斯点数目 % 2. 节点与坐标 nodes linspace(0, L, N); u zeros(N, 1); % 3. 组装刚度矩阵与载荷向量 K zeros(N, N); F zeros(N, 1); for i 1 : N-1 xa nodes(i); xb nodes(i1); gauss_pts (xaxb)/2 (xb-xa)/2 * [-0.577350269189626, 0.577350269189626]; gauss_w (xb-xa)/2 * [1, 1]; for g 1 : n_gauss xg gauss_pts(g); [phi, dphi] mls_shape(xg, nodes, h, p_order, dmax); % 高斯点处组装贡献 for I 1 : N for J 1 : N K(I,J) K(I,J) dphi(I) * dphi(J) * gauss_w(g); end F(I) F(I) phi(I) * f_source(xg) * gauss_w(g); end end end % 4. 左端本质边界 u(0)0, 罚函数法 xb 0; [phi_b, ~] mls_shape(xb, nodes, h, p_order, dmax); K K alpha * (phi_b * phi_b); F F - alpha * phi_b * 0; % 强制 u0 % 5. 求解 u K \ F;代码逻辑走向很清晰。第一步生成高斯点一维线性单元的经典做法是两点高斯积分权系数和坐标都是固定的。第二步在每一个高斯点上调用mls_shape拿到所有节点对这个点的形函数值和导数值再把贡献累加进K和F。第三步处理本质边界条件这是无网格法与有限元最不同的地方因为形函数u(0)Σφ_I(0)·u_I不等于u_0对应的那个节点值必须把整个形函数向量乘以罚系数alpha加进K矩阵。罚系数取1e8是常用值取小了边界条件“不够硬”取太大又会让矩阵条件数恶化这个问题我会在第4章展开讲。3.3 边界条件的三类实现方式对比本质上无网格法处理本质边界条件有三条路子罚函数法、Lagrange乘子法、修正变分法。我实际用下来罚函数法最直观、代码最容易改也是这份程序包大概率采用的方式Lagrange乘子法精度更高但会引入额外的乘子变量矩阵变成不定矩阵求解器选择受限修正变分法理论上最美实现起来却需要对变分形式做改动不适合教学代码。方法实现难度精度表现矩阵性质适用场景罚函数法低依赖罚系数精度一般对称正定快速验证、教学Lagrange乘子法中边界精度高鞍点不定精度要求高的计算修正变分法较高稳定对称正定研究级代码如果你拿到的程序里写的是罚函数法你想换乘子法的话改动点集中在组装部分的“边界积分”上需要在边界节点处组装额外的乘子变量块。不过一维情况下意义不大因为罚函数法配合一个足够大的alpha误差通常已经小于形函数本身的插值误差再折腾未必值得。4. 影响精度的四个参数数值实验与调试策略4.1 节点数与高斯点数收敛趋势怎么看拿到程序后第一件事不是改参数而是跑一遍基础算例记录误差再系统性地扫参数。节点数N是最直接的收敛控制参数线性基MLS形函数的L2误差理论收敛阶是2阶意思是节点数翻倍误差大约降到原来的四分之一。二次基对应3阶。我习惯这样验证收敛阶% 参数扫描: 节点数从 5 到 40, 记录 L2 误差 N_list [5, 8, 12, 16, 20, 30, 40]; L2_err zeros(size(N_list)); for i 1 : length(N_list) [u_h, err] run_efg_1d(N_list(i), 1, 2.5); % 主程序封装成函数 L2_err(i) err; end loglog(N_list, L2_err, o-); % 双对数坐标下斜率就是收敛阶收敛阶的检查方法是看双对数曲线斜率线性基理论斜率是2实测在1.8到2.0之间都属于正常。如果斜率明显偏低甚至为负说明误差主导项不是插值误差而是边界罚函数项或者数值积分误差这时候你调节点数是没有用的。高斯点数方面一维问题取2个点已经足够——因为弱形式里只有一阶导数乘积被积函数是二次多项式2点高斯积分可以精确积分到3次多项式满足要求。如果你把n_gauss改到3精度几乎不变但计算时间涨了50%。4.2 影响域半径dmax那个最“玄学”的参数dmax是这份程序包里最值得你亲手去测的参数。它代表影响域半径相对于节点间距的倍数。dmax太小高斯点周围节点数不足A矩阵奇异dmax太大每个点的形函数要涉及几十个节点刚度矩阵带宽暴涨而且数值积分误差反而增大。dmax影响域内节点数L2误差趋势风险1.22~3可能直接报错A矩阵奇异1.83~5误差偏大边界节点覆盖不够2.55~7误差较优常用起点3.57~9误差略有上升矩阵带宽增大血泪经验是线性基从dmax2.5起步二次基从dmax3.0起步这是大多数算例的甜点区。不要试图把dmax调到1.5以下追求“稀疏”——对无网格法来说稀疏性一直是短板过于追求只会换来奇异矩阵警告。还有一个细节容易忽视dmax取值的密度依赖性。均匀节点分布下dmax2.5很稳但如果你在求解域里加了加密节点局部节点间距变小影响域半径也随之变小。非均匀布点时的dmax选择需要单独做一次敏感性测试别拿着均匀网格的经验硬套。4.3 基函数阶次从线性基升到二次基的代价与收益把程序里的p_order从1改成2基函数从[1,x]变成[1,x,x²]最直接的变化是每一个高斯点处的A矩阵从2×2变成3×3求逆计算量上升但收敛阶从2阶升到3阶。听起来很划算但有个隐性代价二次基要求影响域内至少4个节点如果dmax维持2.5不变有些边界附近的高斯点可能凑不齐节点。这会导致程序在边界处报错或者精度骤降。二次基的另一个问题是导数差分会更敏感。一维MLS导数如果还是用简单的中心差分近似二次基的导数误差会比线性基大因为形函数本身变化更剧烈。我的建议是只有当你确认误差瓶颈是插值逼近而不是边界处理时才考虑升阶在一维算例里线性基好的边界处理通常已经够写一篇完整分析报告了。5. 伽辽金型无网格法常见问题排查五个踩坑现场5.1 现象边界处位移跟精确解差出一大截内部倒还正常原因MLS形函数不满足克罗内克δ性质你直接把u(0)0写成了u(1)0节点索引对应的那个分量等于边界条件根本没加进去或者加错了位置。有限元里改一行代码就能施加本质边界条件的习惯在无网格法里是失效的。解决检查主程序是否有罚函数项K K alpha * (phi_b * phi_b)。如果程序里直接是K(1,1)1e15这种写法多半是用了“有限元思维”处理边界赶紧改成基于形函数向量的罚函数形式。另外边界节点坐标必须是0和L对应的精确值不是离它最近的节点坐标。5.2 现象运行时报“矩阵接近奇异或缩放错误”有时直接NaN原因高斯点处的影响域内节点数少于基函数维度。最常见的情况是dmax设太小或者是边界附近的高斯点的影响域被求解域边界切掉了一半。另一个隐秘原因是权函数在影响域边界处恰好衰减到0虽然节点在影响域内但权重为0实际参与计算的节点数少了。解决先做一次诊断在每个高斯点处打印影响域内节点数和有效权重节点数。然后加保护语句if m p_order 1, error(dmax 过小), end。一般把dmax调到2.5以上就能解决大半问题。别指望靠MATLAB的pinv函数“兜底”那会掩盖问题算出来的结果你也不敢信。5.3 现象程序跑通了但L2误差一直降不下去收敛阶几乎是平的原因误差计算方式有问题。很多人直接在节点位置上算误差范数但一维无网格法的节点数就那么二十来个节点处误差不能代表整个求解域的误差水平。正确做法是在求解域内加密采样点比如取500个均匀点在每一点重构u_h(x)Σφ_I(x)u_I再用数值积分计算L2范数。解决误差分析脚本里采样点数量要远大于节点数通常取5到10倍。具体做法% 用密集采样点计算L2相对误差, 而不是只查节点值 x_fine linspace(0, L, 500); u_exact exact_solution(x_fine); u_h zeros(size(x_fine)); for i 1 : length(x_fine) [phi, ~] mls_shape(x_fine(i), nodes, h, p_order, dmax); u_h(i) phi * u; % 重构近似解 end L2_err sqrt(sum((u_h - u_exact).^2) / sum(u_exact.^2));采样点数量加密后收敛阶会重新回到理论值附近。这个坑很多人会踩因为有限元后处理里直接在节点上比较是习惯动作但在无网格法里行不通。5.4 现象zip文件解压时报错或者解压后.m文件里的中文注释全是乱码原因两种可能。一是zip包本身被分享时用了伪加密或其他兼容性处理常规解压工具识别不了文件头。二是文件内编码是GBK而新版MATLAB默认采用UTF-8读取源文件中文注释就会乱码。解决伪加密的情况先换个工具解压WinRAR、7-Zip都能修复这类问题乱码问题如果在MATLAB 2023a及以上版本出现可以把文件用记事本打开后另存为UTF-8带BOM编码。注意程序能跑就行注释乱码不影响执行逻辑不要因为乱码就重新敲一遍整个文件——那样反而容易引入新错误。5.5 现象程序在别人电脑上跑得好好的自己这里一跑就报错原因MATLAB当前路径不对或者同名函数文件冲突。如果工作目录里有其他名字为mls_shape.m或相似名称的文件MATLAB会优先调用路径上更靠前的文件可能调用了另一个不兼容的版本。解决用which mls_shape命令看一下当前MATLAB实际调用的文件路径确认和你解压目录一致。这个是你在命令行窗口最值得养成的一个检查习惯不只在无网格法里任何多文件MATLAB工程都会遇到。6. 让这套程序不止跑通算例验证与二维扩展的实践路径6.1 验证三件套精确解对比、收敛阶曲线、误差分布图程序跑通只是最低标准。我拿到新写的无网格程序都会强制走一遍三件套验证流程。第一选一个带解析解的一维问题比如u(x)sin(πx)加上k(x)1的均匀扩散方程把计算解和精确解画在同一张图里目测是否吻合。第二计算L2相对误差并画出收敛阶曲线确认在双对数坐标下斜率与理论阶次一致这一步能发现90%的“看起来对但实际错”的隐性bug。第三画误差分布曲线e(x)u_h(x)-u_exact(x)如果误差在两端的峰值远大于中间那就是边界处理有问题如果误差呈现明显震荡那就是数值积分点数不足或者dmax过大。这三步走完代码才算真正属于你了。之后可以验证这类程序在非线性源项问题上的表现把f_source函数里加一个u的依赖项检验迭代格式是否需要更新这些改动都算在“扩展练习”范围内。6.2 从一维走向二维的扩展路径如果你接下来想碰真实的二维无网格问题现在这套一维代码就是你最好的起点。二维扩展需要改动的点很集中一是节点坐标从x变成(x,y)影响域从线段变成圆盘二是基函数向量从[1,x]变成[1,x,y]或加上xy交叉项影响域内最小节点数要求变成6个以上三是高斯积分从一维区间变成二维四边形单元但因为没有网格积分背景单元的划分自由度很大四是导数计算从对x求导变成对x和y分别求偏导可以用数值差分实现。别在这些改动完成之前就动手“猜代码”。我的习惯是每改完一个变量就回归一次一维算例确保改动没有破坏原有逻辑然后再叠加新维度。二维无网格法的调试成本远高于一维任何一处形函数计算错误都会让整个误差分析完全失真到时候你根本分不清是边界问题、积分问题还是形函数问题。从那以后我每次拿到这类教学程序包都会先跑通基准算例、加参数扫描、做收敛阶回归一气走完再考虑扩展。希望这篇笔记能帮你在无网格法的入口处少踩几个坑。本文还有配套的精品资源点击获取 一维伽辽金型无网格法MATLAB程序包：从MLS形函数到边界条件处理
简介这是一份面向无网格法与计算力学初学者的MATLAB源码资源对应一维伽辽金型无网格法的基础编程实现可用于学习从形函数构造到方程求解的完整流程。压缩包内共2个文件其中m文件为主要程序涉及移动最小二乘形函数、伽辽金弱形式离散、刚度矩阵与荷载向量组装、边界条件处理等关键模块另一个rar文件为配套的算例或辅助文件整套资源仅8KB体量轻、结构紧凑便于逐行阅读和调试。借助该代码读者可以快速掌握一维无网格问题的编程框架理解基函数、权函数、节点影响域等概念对结果的影响并以此为基础扩展到二维或非线性问题。目前已有646人学习浏览适合本科生、研究生及MATLAB编程爱好者用于数值方法课程设计、科研入门或相关课题二次开发。1. 一维伽辽金型无网格法一个MATLAB程序包能帮你少走多少弯路先说结论无网格法不是不用网格而是不用“单元”这个概念。你手里这份资源是教材第三章配套的一维伽辽金型无网格法MATLAB程序解压出来就是一套能跑的、从节点生成到误差输出全流程的脚本。它解决的核心问题很直接当你不想碰有限元里那套网格剖分、想试试只用散乱节点做数值逼近时这个程序就是你从头复现无网格法的最短路径。适合谁刚接触无网格法的研究生、打算用无网格法写课程设计或小论文的初学者以及那些已经看过公式但卡在“代码该从哪写起”的从业者。你不需要先精通MATLAB能把脚本跑通、能改参数就已经比只看公式强很多。2. 一维伽辽金型无网格法的核心机制有限元形函数与MLS形函数的本质差异2.1 为什么从一维起步MLS形函数到底在算什么无网格法的种类很多伽辽金型Element-Free Galerkin MethodEFG是其中最成熟的一支它以移动最小二乘Moving Least Squares, MLS形函数作为近似函数。你可能会问既然最终都要解一个类似有限元的线性方程组那它跟有限元到底差在哪答案是形函数。有限元形函数是在单元内定义的形函数只依赖单元内的几何形状而且形函数在节点上满足克罗内克δ性质第I个形函数在节点I处取值为1在其他节点处取值为0。这意味着有限元的未知量就是节点处的真实位移值边界条件可以直接“硬塞”进方程。MLS形函数则完全不同。它基于一组离散节点在每个求解点x处做一次加权最小二乘拟合得到该点的形函数值。换个说法MLS形函数是“看周围若干节点脸色”的它在节点处取值既不是1也不是0是某个中间数值。所以MLS形函数不满足克罗内克δ性质这是无网格法一切边界处理麻烦的根源。下面这张表帮你看清区别。特征有限元形函数MLS无网格形函数定义域局部单元影响域内所有节点节点处取值满足克罗内克δ不满足连续性控制受单元阶次限制由权函数直接决定本质边界条件施加直接修改方程罚函数/Lagrange乘子网格依赖强依赖完全不需要网格由于MLS形函数在任意点都至少有C0甚至C1连续性它不需要像有限元那样做节点协调这是很多人愿意学它的原因。但代价就是求函数值时必须做矩阵求逆计算成本比有限元形函数高了不少。2.2 权函数与影响域形函数计算的Matlab实现MLS形函数的计算链路是给定求解点x的一组基函数和周围节点先组装矩方程再解出形函数。一维情形下基函数p(x)通常取线性基或者二次基权函数则采用三次样条权函数。我一般建议直接写一个独立的形函数函数后续改参数不需要翻主程序。function [phi, dphi] mls_shape(x, nodes, h, p_order, dmax) % x : 求解点坐标 % nodes : 所有节点坐标向量 % h : 节点平均间距 % p_order: 基函数阶次, 1线性基 [1,x], 2二次基 [1,x,x^2] % dmax : 影响域无量纲半径, 实际影响域半径 R dmax * h R dmax * h; n length(nodes); % --- 选取影响域内的节点 --- idx find(abs(nodes - x) R); m length(idx); if m p_order 1 error(影响域内节点数不足, 请调大 dmax); end % --- 组装基函数矩阵和权函数对角阵 --- P zeros(m, p_order1); W zeros(m, m); for j 1:m dx x - nodes(idx(j)); P(j,:) dx.^(0:p_order); % 基函数在局部坐标下计算 s abs(dx) / R; w 1 - 3*s^2 2*s^3; % 三次样条权函数 W(j,j) w; end % --- 计算 A 矩阵和 B 向量 --- A P * W * P; pt (x - nodes(idx)).^0; % 当前点的基函数向量 B P * W; invA inv(A); phi pt * invA * B; % 形函数行向量 % --- 导数用五点差分近似, 教学代码最省事的做法 --- eps1 1e-6; [phi_plus, ~] mls_shape(xeps1, nodes, h, p_order, dmax); [phi_minus, ~] mls_shape(x-eps1, nodes, h, p_order, dmax); dphi (phi_plus - phi_minus) / (2*eps1); end这段代码的逻辑我拆成三步讲。第一步是选节点以求解点x为圆心、半径Rdmax*h画一个圆圈进来的节点才参与计算。第二步是组装A矩阵A矩阵是一个(p_order1)维方阵它本质是加权最小二乘的法方程矩阵如果影响域内节点太少A矩阵会奇异这就是程序报错“矩阵接近奇异”的最常见原因。第三步是求形函数phi最终是一个行向量长度等于影响域内节点数对应每个节点对当前点x的“贡献权重”。别忘了权函数在s1处恰好衰减为0所以影响域边界上的节点权重是零。这也是一个隐患如果高斯点恰好落在两个节点影响域的交界上权函数可能出现数值零导致A矩阵不可逆。后面避坑章节我会专门讲。2.3 Galerkin弱式与整体刚度矩阵的组装流程形函数只是工具真正要解的是一维扩散方程-d/dx(k·du/dx) f。这里的k是材料参数f是源项。标准的伽辽金做法是在方程两边乘以测试函数然后分部积分得到弱形式。弱形式的核心变化是原来方程里有二阶导数分部积分后变成一阶导数乘积的积分加上边界项。边界上自然边界条件如力边界会直接进入边界积分项这跟有限元完全一致。对一维区间[0,L]刚度矩阵的表达式为K_IJ ∫₀ᴸ k(x)·φ₁(x)·φ_J(x)dx载荷向量为F_I ∫₀ᴸ f(x)·φ_I(x)dx再加上边界产生的修正项。矩阵/向量维度物理含义K (刚度矩阵)N×N节点之间的刚度耦合关系F (载荷向量)N×1等效节点载荷A (矩方程)(p1)×(p1)每个求解点处的加权最小二乘法方程Φ (形函数)N×1 (对单点)各节点对当前点的权重这里的K矩阵不再是有限元那种“带状”结构因为每个高斯点都跟影响域内所有节点相连节点的关联范围变大带宽由dmax决定。这也是无网格法比有限元计算量大的原因之一。把dmax从2.0调到3.5矩阵的非零元素数量会明显增加求解时间跟着涨。3. 主程序逐段拆解从解压到跑出第一条曲线3.1 解压后的文件结构与运行顺序这份资源下下来是zip包解压后你会看到几个.m脚本。我建议先不要急着在MATLAB里双击运行先把目录结构看清楚。常规组织方式是一个主程序加若干函数文件主程序负责定义模型参数、节点生成、调用组装函数、求解并画图形函数函数负责计算单点形函数还有误差分析脚本负责输出收敛阶。用命令行进入目录会比较稳避免MATLAB当前路径不对导致函数找不到。# 在任意终端里解压并查看文件结构 cd ~/Downloads unzip 一维伽辽金型无网格法MATLAB程序.zip -d efg_1d ls -l efg_1d/*.m如果你在Windows上解压直接右键“解压到当前文件夹”也行唯一需要注意的是中文文件名和中文注释的编码问题这个放到第5章讲。进入MATLAB后把当前文件夹切换到解压目录然后直接运行主程序。如果你的MATLAB版本偏低比如R2016a大概率也能跑因为程序里没用到太新的语法如果用的是2023a及以上版本反而要留意中文注释乱码问题低版本时反而没有这个毛病。这里提醒一句MATLAB在线网页版mathworks.com的MATLAB Online也能跑这类脚本但它对文件读写路径的管理比较特殊你要先把整个zip在本地解压后逐个上传比较麻烦不如本地装一个桌面版省心。3.2 刚度矩阵组装与边界条件处理核心代码主程序核心段长这样。我按工程实战习惯把“模型参数定义 - 节点生成 - 高斯积分组装 - 边界处理 - 求解画图”分成几块彼此用注释隔开方便你改参数时不会滚到乱七八糟的代码里。% 1. 模型参数 L 1.0; % 求解域长度 N 21; % 节点数 h L / (N - 1); % 节点间距 p_order 1; % 基函数阶次 dmax 2.5; % 影响域半径系数 alpha 1e8; % 罚函数系数 n_gauss 2; % 每个区间的高斯点数目 % 2. 节点与坐标 nodes linspace(0, L, N); u zeros(N, 1); % 3. 组装刚度矩阵与载荷向量 K zeros(N, N); F zeros(N, 1); for i 1 : N-1 xa nodes(i); xb nodes(i1); gauss_pts (xaxb)/2 (xb-xa)/2 * [-0.577350269189626, 0.577350269189626]; gauss_w (xb-xa)/2 * [1, 1]; for g 1 : n_gauss xg gauss_pts(g); [phi, dphi] mls_shape(xg, nodes, h, p_order, dmax); % 高斯点处组装贡献 for I 1 : N for J 1 : N K(I,J) K(I,J) dphi(I) * dphi(J) * gauss_w(g); end F(I) F(I) phi(I) * f_source(xg) * gauss_w(g); end end end % 4. 左端本质边界 u(0)0, 罚函数法 xb 0; [phi_b, ~] mls_shape(xb, nodes, h, p_order, dmax); K K alpha * (phi_b * phi_b); F F - alpha * phi_b * 0; % 强制 u0 % 5. 求解 u K \ F;代码逻辑走向很清晰。第一步生成高斯点一维线性单元的经典做法是两点高斯积分权系数和坐标都是固定的。第二步在每一个高斯点上调用mls_shape拿到所有节点对这个点的形函数值和导数值再把贡献累加进K和F。第三步处理本质边界条件这是无网格法与有限元最不同的地方因为形函数u(0)Σφ_I(0)·u_I不等于u_0对应的那个节点值必须把整个形函数向量乘以罚系数alpha加进K矩阵。罚系数取1e8是常用值取小了边界条件“不够硬”取太大又会让矩阵条件数恶化这个问题我会在第4章展开讲。3.3 边界条件的三类实现方式对比本质上无网格法处理本质边界条件有三条路子罚函数法、Lagrange乘子法、修正变分法。我实际用下来罚函数法最直观、代码最容易改也是这份程序包大概率采用的方式Lagrange乘子法精度更高但会引入额外的乘子变量矩阵变成不定矩阵求解器选择受限修正变分法理论上最美实现起来却需要对变分形式做改动不适合教学代码。方法实现难度精度表现矩阵性质适用场景罚函数法低依赖罚系数精度一般对称正定快速验证、教学Lagrange乘子法中边界精度高鞍点不定精度要求高的计算修正变分法较高稳定对称正定研究级代码如果你拿到的程序里写的是罚函数法你想换乘子法的话改动点集中在组装部分的“边界积分”上需要在边界节点处组装额外的乘子变量块。不过一维情况下意义不大因为罚函数法配合一个足够大的alpha误差通常已经小于形函数本身的插值误差再折腾未必值得。4. 影响精度的四个参数数值实验与调试策略4.1 节点数与高斯点数收敛趋势怎么看拿到程序后第一件事不是改参数而是跑一遍基础算例记录误差再系统性地扫参数。节点数N是最直接的收敛控制参数线性基MLS形函数的L2误差理论收敛阶是2阶意思是节点数翻倍误差大约降到原来的四分之一。二次基对应3阶。我习惯这样验证收敛阶% 参数扫描: 节点数从 5 到 40, 记录 L2 误差 N_list [5, 8, 12, 16, 20, 30, 40]; L2_err zeros(size(N_list)); for i 1 : length(N_list) [u_h, err] run_efg_1d(N_list(i), 1, 2.5); % 主程序封装成函数 L2_err(i) err; end loglog(N_list, L2_err, o-); % 双对数坐标下斜率就是收敛阶收敛阶的检查方法是看双对数曲线斜率线性基理论斜率是2实测在1.8到2.0之间都属于正常。如果斜率明显偏低甚至为负说明误差主导项不是插值误差而是边界罚函数项或者数值积分误差这时候你调节点数是没有用的。高斯点数方面一维问题取2个点已经足够——因为弱形式里只有一阶导数乘积被积函数是二次多项式2点高斯积分可以精确积分到3次多项式满足要求。如果你把n_gauss改到3精度几乎不变但计算时间涨了50%。4.2 影响域半径dmax那个最“玄学”的参数dmax是这份程序包里最值得你亲手去测的参数。它代表影响域半径相对于节点间距的倍数。dmax太小高斯点周围节点数不足A矩阵奇异dmax太大每个点的形函数要涉及几十个节点刚度矩阵带宽暴涨而且数值积分误差反而增大。dmax影响域内节点数L2误差趋势风险1.22~3可能直接报错A矩阵奇异1.83~5误差偏大边界节点覆盖不够2.55~7误差较优常用起点3.57~9误差略有上升矩阵带宽增大血泪经验是线性基从dmax2.5起步二次基从dmax3.0起步这是大多数算例的甜点区。不要试图把dmax调到1.5以下追求“稀疏”——对无网格法来说稀疏性一直是短板过于追求只会换来奇异矩阵警告。还有一个细节容易忽视dmax取值的密度依赖性。均匀节点分布下dmax2.5很稳但如果你在求解域里加了加密节点局部节点间距变小影响域半径也随之变小。非均匀布点时的dmax选择需要单独做一次敏感性测试别拿着均匀网格的经验硬套。4.3 基函数阶次从线性基升到二次基的代价与收益把程序里的p_order从1改成2基函数从[1,x]变成[1,x,x²]最直接的变化是每一个高斯点处的A矩阵从2×2变成3×3求逆计算量上升但收敛阶从2阶升到3阶。听起来很划算但有个隐性代价二次基要求影响域内至少4个节点如果dmax维持2.5不变有些边界附近的高斯点可能凑不齐节点。这会导致程序在边界处报错或者精度骤降。二次基的另一个问题是导数差分会更敏感。一维MLS导数如果还是用简单的中心差分近似二次基的导数误差会比线性基大因为形函数本身变化更剧烈。我的建议是只有当你确认误差瓶颈是插值逼近而不是边界处理时才考虑升阶在一维算例里线性基好的边界处理通常已经够写一篇完整分析报告了。5. 伽辽金型无网格法常见问题排查五个踩坑现场5.1 现象边界处位移跟精确解差出一大截内部倒还正常原因MLS形函数不满足克罗内克δ性质你直接把u(0)0写成了u(1)0节点索引对应的那个分量等于边界条件根本没加进去或者加错了位置。有限元里改一行代码就能施加本质边界条件的习惯在无网格法里是失效的。解决检查主程序是否有罚函数项K K alpha * (phi_b * phi_b)。如果程序里直接是K(1,1)1e15这种写法多半是用了“有限元思维”处理边界赶紧改成基于形函数向量的罚函数形式。另外边界节点坐标必须是0和L对应的精确值不是离它最近的节点坐标。5.2 现象运行时报“矩阵接近奇异或缩放错误”有时直接NaN原因高斯点处的影响域内节点数少于基函数维度。最常见的情况是dmax设太小或者是边界附近的高斯点的影响域被求解域边界切掉了一半。另一个隐秘原因是权函数在影响域边界处恰好衰减到0虽然节点在影响域内但权重为0实际参与计算的节点数少了。解决先做一次诊断在每个高斯点处打印影响域内节点数和有效权重节点数。然后加保护语句if m p_order 1, error(dmax 过小), end。一般把dmax调到2.5以上就能解决大半问题。别指望靠MATLAB的pinv函数“兜底”那会掩盖问题算出来的结果你也不敢信。5.3 现象程序跑通了但L2误差一直降不下去收敛阶几乎是平的原因误差计算方式有问题。很多人直接在节点位置上算误差范数但一维无网格法的节点数就那么二十来个节点处误差不能代表整个求解域的误差水平。正确做法是在求解域内加密采样点比如取500个均匀点在每一点重构u_h(x)Σφ_I(x)u_I再用数值积分计算L2范数。解决误差分析脚本里采样点数量要远大于节点数通常取5到10倍。具体做法% 用密集采样点计算L2相对误差, 而不是只查节点值 x_fine linspace(0, L, 500); u_exact exact_solution(x_fine); u_h zeros(size(x_fine)); for i 1 : length(x_fine) [phi, ~] mls_shape(x_fine(i), nodes, h, p_order, dmax); u_h(i) phi * u; % 重构近似解 end L2_err sqrt(sum((u_h - u_exact).^2) / sum(u_exact.^2));采样点数量加密后收敛阶会重新回到理论值附近。这个坑很多人会踩因为有限元后处理里直接在节点上比较是习惯动作但在无网格法里行不通。5.4 现象zip文件解压时报错或者解压后.m文件里的中文注释全是乱码原因两种可能。一是zip包本身被分享时用了伪加密或其他兼容性处理常规解压工具识别不了文件头。二是文件内编码是GBK而新版MATLAB默认采用UTF-8读取源文件中文注释就会乱码。解决伪加密的情况先换个工具解压WinRAR、7-Zip都能修复这类问题乱码问题如果在MATLAB 2023a及以上版本出现可以把文件用记事本打开后另存为UTF-8带BOM编码。注意程序能跑就行注释乱码不影响执行逻辑不要因为乱码就重新敲一遍整个文件——那样反而容易引入新错误。5.5 现象程序在别人电脑上跑得好好的自己这里一跑就报错原因MATLAB当前路径不对或者同名函数文件冲突。如果工作目录里有其他名字为mls_shape.m或相似名称的文件MATLAB会优先调用路径上更靠前的文件可能调用了另一个不兼容的版本。解决用which mls_shape命令看一下当前MATLAB实际调用的文件路径确认和你解压目录一致。这个是你在命令行窗口最值得养成的一个检查习惯不只在无网格法里任何多文件MATLAB工程都会遇到。6. 让这套程序不止跑通算例验证与二维扩展的实践路径6.1 验证三件套精确解对比、收敛阶曲线、误差分布图程序跑通只是最低标准。我拿到新写的无网格程序都会强制走一遍三件套验证流程。第一选一个带解析解的一维问题比如u(x)sin(πx)加上k(x)1的均匀扩散方程把计算解和精确解画在同一张图里目测是否吻合。第二计算L2相对误差并画出收敛阶曲线确认在双对数坐标下斜率与理论阶次一致这一步能发现90%的“看起来对但实际错”的隐性bug。第三画误差分布曲线e(x)u_h(x)-u_exact(x)如果误差在两端的峰值远大于中间那就是边界处理有问题如果误差呈现明显震荡那就是数值积分点数不足或者dmax过大。这三步走完代码才算真正属于你了。之后可以验证这类程序在非线性源项问题上的表现把f_source函数里加一个u的依赖项检验迭代格式是否需要更新这些改动都算在“扩展练习”范围内。6.2 从一维走向二维的扩展路径如果你接下来想碰真实的二维无网格问题现在这套一维代码就是你最好的起点。二维扩展需要改动的点很集中一是节点坐标从x变成(x,y)影响域从线段变成圆盘二是基函数向量从[1,x]变成[1,x,y]或加上xy交叉项影响域内最小节点数要求变成6个以上三是高斯积分从一维区间变成二维四边形单元但因为没有网格积分背景单元的划分自由度很大四是导数计算从对x求导变成对x和y分别求偏导可以用数值差分实现。别在这些改动完成之前就动手“猜代码”。我的习惯是每改完一个变量就回归一次一维算例确保改动没有破坏原有逻辑然后再叠加新维度。二维无网格法的调试成本远高于一维任何一处形函数计算错误都会让整个误差分析完全失真到时候你根本分不清是边界问题、积分问题还是形函数问题。从那以后我每次拿到这类教学程序包都会先跑通基准算例、加参数扫描、做收敛阶回归一气走完再考虑扩展。希望这篇笔记能帮你在无网格法的入口处少踩几个坑。本文还有配套的精品资源点击获取