22.4 使用稀疏矩阵的实际应用示例¶
稀疏矩阵的一个常见应用是在有限元模型的求解中。有限元模型能够对没有解析解的偏微分方程进行数值求解，这通常是由于求解域的几何形状非常复杂。
为了说明这一应用，我们考虑边值拉普拉斯方程。该方程组可以模拟标量势场，例如热势或电势。给定一个区域 Omega，其边界为 dOmega。在 dOmega 上的所有点，边界条件都是已知的，我们希望计算 Omega 内部的势。边界条件可以指定势（狄利克雷边界条件）、势在边界上的法向导数（诺伊曼边界条件），或势与其导数的加权和（柯西边界条件）。
在热模型中，我们想要计算 Omega 内部的温度，并已知边界温度（狄利克雷条件）或热通量（通过除以边界处的热导率，可由热通量计算出诺伊曼条件）。类似地，在电模型中，我们想要计算 Omega 内部的电压，并已知边界电压（狄利克雷条件）或电流（通过除以电导率后的诺伊曼条件）。在电模型中，大部分边界通常是电隔离的——这对应于电流为零的诺伊曼边界条件。
最简单的有限元模型会将 Omega 划分为单纯形（二维为三角形，三维为四面体）。我们以 EIDORS 项目中的三维示例为例：一个装有液体的圆柱形储罐，内部有一个小的非导电球11。该模型旨在反映电阻抗断层扫描的应用，其中对这样的储罐施加电流模式，以对内部电导率分布进行成像。为了描述有限元几何，我们有一个顶点矩阵 nodes 和一个单纯形矩阵 elems。
以下示例创建了一个简单的矩形二维导电介质，在相对两侧分别施加 10 V 和 20 V 电压（狄利克雷边界条件）。所有其他边均为电隔离。
node_y = [1;1.2;1.5;1.8;2]*ones(1,11);
node_x = ones(5,1)*[1,1.05,1.1,1.2, ...
1.3,1.5,1.7,1.8,1.9,1.95,2];
nodes = [node_x(:), node_y(:)];
[h,w] = size (node_x);
elems = [];
for idx = 1:w-1
widx = (idx-1)*h;
elems = [elems; ...
widx+[(1:h-1);(2:h);h+(1:h-1)]'; ...
widx+[(2:h);h+(2:h);h+(1:h-1)]' ];
endfor
E = size (elems,1); # 单纯形数量
N = size (nodes,1); # 顶点数量
D = size (elems,2); # 维度+1
这将创建一个 N×2 的矩阵 nodes 和一个 E×3 的矩阵 elems，其中的值定义了有限元三角形：
nodes(1:7,:)'
1.00 1.00 1.00 1.00 1.00 1.05 1.05 ...
1.00 1.20 1.50 1.80 2.00 1.00 1.20 ...
elems(1:7,:)'
1 2 3 4 2 3 4 ...
2 3 4 5 7 8 9 ...
6 7 8 9 6 7 8 ...
使用一阶有限元法，我们将 Omega 内部的电导率分布近似为每个单纯形上的常数（由向量 conductivity 表示）。基于有限元几何，我们首先计算每个单纯形的系统（或刚度）矩阵（在单元系统矩阵 SE 的对角线上表示为 3×3 的元素）。基于 SE 和一个 N×DE 的连接矩阵 C（表示单纯形与顶点之间的连接关系），计算全局系统矩阵 S。
## 单元电导率
conductivity = [1*ones(1,16), ...
2*ones(1,48), 1*ones(1,16)];
## 连接矩阵
C = sparse ((1:D*E), reshape (elems', ...
D*E, 1), 1, D*E, N);
## 计算系统矩阵
Siidx = floor ([0:D*E-1]'/D) * D * ...
ones(1,D) + ones(D*E,1)*(1:D) ;
Sjidx = [1:D*E]'*ones (1,D);
Sdata = zeros (D*E,D);
dfact = factorial (D-1);
for j = 1:E
a = inv ([ones(D,1), ...
nodes(elems(j,:), :)]);
const = conductivity(j) * 2 / ...
dfact / abs (det (a));
Sdata(D*(j-1)+(1:D),:) = const * ...
a(2:D,:)' * a(2:D,:);
endfor
## 单元系统矩阵
SE = sparse(Siidx,Sjidx,Sdata);
## 全局系统矩阵
S = C'* SE *C;
系统矩阵 S 在欧姆定律 S * V = I 中起到电导率的作用。基于狄利克雷和诺伊曼边界条件，我们可以求解每个顶点的电压 V。
## 狄利克雷边界条件
D_nodes = [1:5, 51:55];
D_value = [10*ones(1,5), 20*ones(1,5)];
V = zeros (N,1);
V(D_nodes) = D_value;
idx = 1:N; # 不含狄利克雷边界条件的顶点
idx(D_nodes) = [];
## 诺伊曼边界条件。注意：
## N_value 必须根据边界长度和单元电导率进行归一化
N_nodes = [];
N_value = [];
Q = zeros (N,1);
Q(N_nodes) = N_value;
V(idx) = S(idx,idx) \ ( Q(idx) - ...
S(idx,D_nodes) * V(D_nodes));
最后，为了显示解，我们在 z 轴上展示每个单纯形顶点的求解电压值。参见图 22.6。
elemx = elems(:,[1,2,3,1])'; xelems = reshape (nodes(elemx, 1), 4, E); yelems = reshape (nodes(elemx, 2), 4, E); velems = reshape (V(elemx), 4, E); plot3 (xelems,yelems,velems,"k"); print "grid.eps";
图 22.6: 示例有限元模型，显示三角形单元。每个顶点的高度对应于该点的求解值。
