22.4使用稀疏矩阵的真实生活示例¶

稀疏矩阵的一个常见应用是在有限元模型的求解中。有限元模型允许数值求解并没有闭合形式解的偏微分方程,这通常是因为域的形状很复杂。

为了促进这一应用,我们考虑了边值拉普拉斯方程。这个系统可以模拟标量势场,例如热或电势。给定一个具有边界dOmega的中等Omega。在dOmega上的所有点上,边界条件都是已知的,我们希望计算电势inOmega。边界条件可以指定势(Dirichlet边界条件)、其在边界上的法向导数(Neumann边界条件)或势与导数的加权和(Cauchy边界条件)。

在热模型中,我们想计算Omegaa中的温度,并知道边界温度(狄利克雷条件)或热通量(从中我们可以通过除以边界处的热导率来计算Neumann条件)。类似地,在电学模型中,我们想计算Omegaa中的电压,并知道边界电压(Dirichlet)或这里的(通过电导率潜水后的Neumann条件)。在电气模型中,边界的大部分被电气隔离是很常见的;这是这里的等于零的Neumann边界条件。

最简单的有限元模型将把10兆分为简单的(2D中的三角形,3D中的金字塔)。我们以EIDORS项目中一个带有小型非导电球的圆柱形充液罐为三维示例¹¹该模型旨在反映电阻抗断层扫描的应用,其中将这里的模式应用于这样的储罐,以对内部电导率分布进行成像。为了描述有限元几何,我们有一个垂直矩阵节点和单纯形elems.

以下示例创建了一个简单的矩形二维电导介质,在相对两侧施加10V和20V(狄利克雷边界条件)。所有其他边缘都进行了电气化处理。

node_y=[1;1.2;1.5;1.8;2]*个一(1,11);node_x=ones(5,1)*[1,1.05,11,1.2,…1.3,1.5,1.7,1.8,1.9,1.95,2];nodes=[node_x(:),node_y(:)];[h,w]=大小(node_x);elems=[];对于idx=1:w-1 widx=(idx-1)*h;elems=[elems;…widx+[(1:h-1);(2:h);h+(1:h-2)]'。。。widx+[(2:h);h+(2:h,h+(1:h-1)]'];endfor E=大小(elems,1);#单形数N=大小(节点,1);#顶点数D=大小(elems,2);#尺寸+1

这将创建一个N乘2的矩阵节点和E-by-3矩阵elems具有定义有限元三角形的值:

节点(1:7,:)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,:)'12 3 4 2 3 4。。。2    3    4    5    7    8    9 ...    6    7    8    9    6    7    8 ...

使用一阶有限元法,我们近似了每个单纯形上Omegaas常数的电导率分布(从向量表示电导率).基于有限元几何,我们首先计算每个单纯形的系统(或刚度)矩阵(在按元系统矩阵的对角线上表示为3乘3个元素)东南方). 基于东南方和N-by-DE连接矩阵C,表示单形和顶点之间的连接,全局连接矩阵S已计算。

##元素电导率=[1*个(1,16),…2*个(1,48),1*个(116];##连通性矩阵C=稀疏((1:D*E),整形(elems',…D*E,1),1,D*E、N);##计算系统矩阵Siidx=楼层([0:D*E-1]'/D)*D*。。。一个(1,D)+一个(D*E,1)*(1:D);Sjidx=[1:D*E]'*ones(1,D);Sdata=零(D*E,D);dfact=因子(D-1);对于j=1:E a=inv([ones(D,1),…nodes(elems(j,:),:)]);const=电导率(j)*2/。。。dfact/abs(det(a));Sdata(D*(j-1)+(1:D),:)=常量*。。。a(2:D,:)'*a(2:D,:);endfor##元素系统矩阵SE=稀疏(Siidx,Sjidx,Sdata);##全局系统矩阵S=C'*SE*C;

系统矩阵的作用类似于电导率S在欧姆定律中S*V=I基于Dirichlet和Neumann边界条件,我们能够求解每个顶点处的电压五、.

##狄利克雷边界条件D_nodes=[1:5,51:55];D_value=[10*ones(1,5),20*ones[1,5)];V=零(N,1);V(D_nodes)=D_value;idx=1:N;#没有Dirichlet的顶点#边界条件idx(D_nodes)=[];##Neumann边界条件。注意,##N_值必须通过##边界长度和元素电导率N_nodes=[]进行归一化;N_值=[];Q=零(N,1);Q(N_nodes)=N_value;V(idx)=S(idx,idx)\(Q(idx。。。S(idx;

最后,为了显示解,我们在z轴上显示每个单纯形顶点的每个求解电压值。详见图22.6.

elemx=elems(:,[1,2,1])';xelems=整形(节点(elemx,1),4,E);yelems=整形(节点(elemx,2),4,E);velems=整形(V(elemx),4,E);情节3(xelems,yelems,velems,“k”);打印“grid.eps”;

图22.6:示例有限元模型显示三角形单元。每个顶点的高度对应于解的值。