MATLAB等几何分析(IGA)程序拆解：从节点向量到刚度矩阵的完整链路与避坑指南

简介：这份资源是面向初学者的MATLAB等几何分析（IGA）程序包nliga_v1.3，用于在MATLAB环境中学习并实践将CAD参数化几何直接用于数值模拟的方法。包内共240个文件，以194个.m脚本与函数文件为主体，辅以12个.cc与7个.c底层NURBS基函数实现、11个.msh网格数据、6个.mexw64编译接口及少量头文件与说明文档，压缩包约12.34MB，结构上覆盖基础函数库、模型构建、求解器与后处理等模块。已有1415人学习下载，说明其在入门群体中具有一定参考价值。读者可借助其中的NURBS曲线曲面构造、B样条基函数计算与示例案例，理解IGA基本原理并动手运行简单分析，逐步掌握MATLAB编程与结构、流体等工程问题数值模拟的衔接方式。

1. 拿到一份 MATLAB 版 IGA 程序，先别急着跑

前阵子帮一个做结构优化的朋友看代码，他丢过来一个压缩包，说“这是别人用 MATLAB 写的等几何分析程序，我跑不通”。我打开一看，目录里全是.m文件，没有 README，没有示例脚本，函数名还都是buildK、assemble、refine这种通用词。这种“裸源码”在工程圈太常见了——作者自己跑得通，换台机器就玄学报错。等几何分析（IGA）本身不复杂：用 NURBS 基函数代替传统有限元的拉格朗日插值，让几何模型和求解空间共用一套样条表示，省掉网格划分这一步。这份 MATLAB 程序就是把 IGA 的完整链路——控制点、节点向量、基函数求值、刚度矩阵组装、边界条件施加、线性方程组求解——用脚本串起来。适合两类人：学过有限元但没碰过样条基的学生，以及想把 IGA 嵌进自己优化流程的工程师。下面按“能干什么、怎么跑、怎么改”拆开讲。

2. 程序骨架拆解：从节点向量到刚度矩阵的完整链路

2.1 目录与函数分工

拿到源码先别逐行读，先看文件命名和调用关系。这类 IGA 程序通常分四层：几何层（控制点坐标、权重、节点向量）、基函数层（B 样条/NURBS 基函数及其导数）、组装层（单元循环、高斯积分、刚度矩阵）、求解层（边界条件、线性求解、后处理）。我见过的典型文件清单长这样：

这个分层不是作者拍脑袋定的，是 IGA 的数学结构决定的：基函数只认节点向量和参数坐标，组装只认基函数值和雅可比，求解只认全局矩阵。你把basisFuns.m单独拎出来测，输入一个参数坐标和节点向量，输出基函数值，对不上就是这里的问题，不用往下查。常见做法是每个函数开头写一段注释说明输入输出维度，但这份源码大概率没有——我一般会先补一个test_basis.m做单元测试，拿一个已知的二次 B 样条节点向量[0 0 0 1 2 2 2]验证。

2.2 节点向量与基函数：IGA 的地基

IGA 和传统 FEM 最大的区别就在这一层。传统 FEM 的形函数是固定的多项式，IGA 的基函数由节点向量决定，改变节点向量就改变了基函数的连续性和支撑域。看这段典型实现：

function [N, dN] = basisFuns(i, u, p, U) % 计算第 i 个 p 次基函数在参数 u 处的值及一阶导数 % i: 基函数索引（1-based） % u: 参数坐标 % p: 次数 % U: 节点向量 N = zeros(p+1, 1); dN = zeros(p+1, 1); left = zeros(p+1, 1); right = zeros(p+1, 1); N(1) = 1.0; for j = 1:p left(j+1) = u - U(i+1-j); right(j+1) = U(i+j) - u; saved = 0.0; for r = 0:j-1 temp = N(r+1) / (right(r+2) + left(j-r+1)); N(r+1) = saved + right(r+2) * temp; saved = left(j-r+1) * temp; end N(j+1) = saved; end % 导数计算略，按同样递推结构展开 end

这段是 Cox-de Boor 递推的原地版本，left和right存的是参数点到左右节点的距离，saved用来滚动保存中间量。参数说明里最容易翻车的是索引：MATLAB 从 1 开始，而文献里的公式从 0 开始，U(i+1-j)这种写法就是在做偏移。如果你改次数p，节点向量长度必须满足length(U) = n + p + 1（n 是控制点数），否则递推直接越界。我一般会在main.m里加一句断言：

assert(length(U) == size(P,1) + p + 1, '节点向量长度与控制点数不匹配');

这句断言能挡掉八成“跑不通”的问题——很多人改了几何参数忘了同步节点向量。

2.3 组装层：单元循环与高斯积分

组装层的逻辑和标准 FEM 几乎一样，区别只在形函数来源。核心循环长这样：

function K = assembleStiffness(P, W, U, V, p, q, ngp) % P: 控制点坐标 (n x 2) % W: 权重 (n x 1) % U, V: 两个参数方向的节点向量 % p, q: 两个方向的次数 % ngp: 每个方向的高斯点数 n = size(P, 1); K = sparse(n, n); [gp, gw] = gaussQuad(ngp); % 高斯点与权重 % 遍历参数域单元（简化：单单元情形） for a = 1:ngp for b = 1:ngp u = 0.5*(1+gp(a)); % 映射到 [0,1] v = 0.5*(1+gp(b)); [R, dRdu, dRdv] = nurbsBasis(u, v, U, V, p, q, W); J = dRdu' * P; % 雅可比矩阵 detJ = det(J); dRdx = J \ [dRdu dRdv]'; % 物理域导数 K = K + (dRdx' * dRdx) * detJ * gw(a) * gw(b); end end end

几个关键点：sparse预分配是必须的，控制点数上千后稠密矩阵直接爆内存；J \ [dRdu dRdv]'是解雅可比线性系统求物理域导数，比显式求逆稳；detJ必须为正，出现负值说明单元参数方向反了或者控制点顺序错了。高斯点数ngp的经验值是ceil((p+q+2)/2)，二次基函数取 3 点足够，取多了浪费，取少了积分不精确——这个坑我在板壳单元上踩过，刚度矩阵不对称，查了半天才发现是积分欠采样。

2.4 边界条件与求解：本质边界条件的特殊处理

IGA 最反直觉的地方在这里：NURBS 基函数不满足 Kronecker delta 性质，控制点上的值不等于物理场在该点的值。所以传统 FEM 里“把对应自由度置零”的做法不能直接用。常见处理有三种：罚函数法、Nitsche 法、控制点插值法。这份源码大概率用的是最粗暴的罚函数：

function [K, F] = applyBC(K, F, fixedDofs, penalty) % fixedDofs: 被约束的自由度索引 % penalty: 罚因子，通常取 1e8 ~ 1e12 for i = 1:length(fixedDofs) d = fixedDofs(i); K(d, d) = K(d, d) + penalty; F(d) = F(d) + penalty * 0; % 齐次边界 end end

罚因子取太小约束不紧，取太大条件数爆炸。我一般取1e10量级，然后检查残差norm(K*x - F)/norm(F)，超过1e-6就说明罚因子不合适或者基函数有问题。如果你要做非齐次边界（比如给定位移），F(d)那行要改成penalty * prescribedValue，但注意这个值不是控制点值，得先做插值换算——这是二次开发时最容易写错的一行。

3. 跑通第一遍：环境、参数与验证

3.1 环境准备与最小运行路径

MATLAB 版本不用太新，R2016b 以上都行，因为代码里大概率没用到字符串数组或隐式扩展这些新特性。不需要额外工具箱，sparse、det、\都是基础功能。跑之前先确认三件事：当前目录切到源码根目录（否则函数找不到）、main.m里的几何参数和节点向量自洽、输出目录有写权限。最小运行路径就是直接敲：

>> main

如果报“未定义函数或变量”，八成是路径问题；如果报“索引超出数组范围”，去查节点向量长度断言；如果跑完但位移全零，查边界条件索引是不是 1-based 对上了。我习惯在main.m末尾加一段自检：

residual = norm(K*x - F) / norm(F); fprintf('相对残差: %.3e\n', residual); assert(residual < 1e-6, '求解残差过大，检查罚因子或基函数');

这段自检能帮你区分“跑完了”和“跑对了”——很多人看到不报错就以为成了，结果位移场是错的。

3.2 用已知解验证：悬臂梁算例

验证 IGA 程序最稳的办法是拿一个解析解已知的算例。二维悬臂梁端部受集中力，Timoshenko 梁理论有解析解，但 IGA 程序多半是平面应力或平面应变，得用更简单的：单向拉伸。左端固定，右端均布载荷，解析解是线性位移场。跑完提取右端控制点位移，和理论值比：

% 理论值：u = F*L/(E*A) u_theory = F_total * L / (E * thickness * height); u_num = max(abs(x(1:2:end))); % 取 x 方向最大位移 err = abs(u_num - u_theory) / u_theory; fprintf('相对误差: %.2f%%\n', err*100);

误差在 1% 以内算正常，超过 5% 就要查：基函数次数够不够、高斯积分够不够、边界条件施加方式对不对。注意 NURBS 的权重要全设为 1 才退化成 B 样条，如果源码里权重是随机给的，先改成全 1 再验证。

3.3 网格细化：h-refinement 与 p-refinement

IGA 的细化比 FEM 灵活，因为可以只改节点向量不改几何。h-refinement 是插入节点，p-refinement 是升阶，k-refinement 是两者结合。源码里如果有refine.m，大概率只实现了 h-refinement 的节点插入：

function [U_new, P_new, W_new] = hRefine(U, P, W, p, numInsert) % 在参数域中点插入 numInsert 个节点 for k = 1:numInsert u_new = 0.5; % 简化：固定中点 [U, P, W] = knotInsert(U, P, W, p, u_new); end U_new = U; P_new = P; W_new = W; end

节点插入的公式（Boehm 算法）不难，但索引偏移容易错。验证方法是插入前后几何不变：在参数域取 100 个点，用插入前后的基函数分别算物理坐标，最大偏差应该小于1e-12。这个测试我每次都跑，因为节点插入写错的话，几何会“悄悄变形”，后面所有结果都不可信。

4. 避坑与排查：五个真实翻车现场

4.1 现象：刚度矩阵不对称

原因：高斯积分点数不足，或者基函数导数在单元边界上求值不连续。IGA 的基函数在单元内部是多项式，高斯积分应该精确，但如果ngp取少了，积分误差会破坏对称性。解决：把ngp提到p+q+2，再检查norm(K - K')/norm(K)，应该小于1e-14。

4.2 现象：位移场出现振荡

原因：罚因子过大导致条件数爆炸，或者 NURBS 权重出现负值。解决：先把罚因子降到1e8试，再检查权重是否全正。NURBS 权重必须为正，这是数学要求，源码里如果从文件读权重，读进来先assert(all(W > 0))。

4.3 现象：细化后结果反而变差

原因：节点插入时控制点更新公式用错，几何发生了微小变形。解决：跑 3.3 节的几何不变性测试，偏差超过1e-10就说明插入算法有问题。我见过一个版本把U(i+1-j)写成了U(i-j)，几何偏差 0.3%，肉眼看不出来，但应力集中区域误差翻倍。

4.4 现象：detJ出现负值

原因：控制点顺序反了，或者参数方向映射搞反。解决：检查P的排列是不是按参数方向递增，二维问题里P通常按u快v慢排列。如果几何是从 CAD 导出的，控制点顺序可能和程序约定不一致，需要重排。

4.5 现象：换台机器结果不一样

原因：MATLAB 版本差异导致sparse求解器默认选项不同，或者浮点累加顺序不同。解决：在main.m开头固定随机种子（如果有随机），线性求解显式指定x = K \ F而不是用迭代求解器。IGA 的刚度矩阵是对称正定的，直接求解最稳。

5. 二次开发：把 IGA 嵌进自己的流程

5.1 改造方向：从单 patch 到多 patch

这份源码大概率只支持单 patch（一个 NURBS 面）。实际工程里几何往往是多 patch 拼的，改造的关键是 patch 间的连续性约束。常见做法是在公共边界上做控制点合并或者用 mortaring 方法。我一般先做最简单的：公共边界上的控制点坐标强制相等，用拉格朗日乘子施加。代码上就是在assembleStiffness之后加一个约束矩阵C，解鞍点系统：

% 多 patch 连续性约束：C*x = 0 K_aug = [K, C'; C, zeros(size(C,1))]; F_aug = [F; zeros(size(C,1),1)]; x_aug = K_aug \ F_aug; x = x_aug(1:size(K,1));

这个改法的代价是系统不再正定，得用\的通用求解。如果 patch 数量多，鞍点系统会变大，可以考虑静凝聚消去内部自由度。

5.2 改造方向：接外部几何数据

从 CAD 导出 NURBS 数据通常是 IGES 格式，MATLAB 没有原生读取函数。我一般让上游导出成简单的文本格式：控制点坐标、权重、节点向量各一段，用dlmread或readmatrix读。读进来第一件事是单位统一——CAD 常用毫米，程序里如果按米算，刚度差1e6倍，结果全错。我习惯在读入后加一句：

scale = 1e-3; % 毫米转米 P = P * scale;

然后跑 3.2 节的拉伸验证，确认量纲对了再往下做。

5.3 验证习惯：每次改动都跑回归

二次开发最怕“改一处坏三处”。我的习惯是维护一个小型回归集：拉伸算例验量纲、悬臂梁算例验精度、细化算例验几何不变性。每次改完basisFuns或assembleStiffness，三个算例全跑一遍，残差和误差都打印出来对比。这套习惯是被坑出来的——有次改罚因子施加方式，拉伸算例过了，悬臂梁误差从 1% 跳到 8%，查了一下午才发现是边界条件索引偏移了一位。从那以后我每次动边界条件相关代码，都强制走一遍三个算例的回归，不跳步。希望这份拆解能帮你少走几个弯路，把这份 IGA 程序真正跑起来、改下去。

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