22.4 疎行列を使った実際の例
疎行列の一般的な用途は、有限要素モデルの解法です。有限要素モデルでは、通常、領域の形状が複雑なために閉じた形式の解を持たない偏微分方程式を数値的に解くことができます。
このアプリケーションを動機付けるために、境界値ラプラス方程式を検討します。このシステムは、熱や電位などのスカラー電位場をモデル化できます。境界 dOmega を持つ媒体 Omega が与えられます。dOmega 上のすべてのポイントで境界条件が既知であり、Omega の電位を計算します。境界条件は、電位 (ディリクレ境界条件)、境界を横切るその法線導関数 (ノイマン境界条件)、または電位とその導関数の加重和 (コーシー境界条件) を指定できます。
熱モデルでは、オメガで温度を計算し、境界温度 (ディリクレ条件) または熱流束 (境界での熱伝導率で割ることでノイマン条件を計算できます) を知りたいと思います。同様に、電気モデルでは、オメガで電圧を計算し、境界電圧 (ディリクレ) または電流 (電気伝導率で割った後のノイマン条件) を知りたいと思います。電気モデルでは、境界の大部分が電気的に分離されているのが一般的です。これは、電流がゼロに等しいノイマン境界条件です。
最も単純な有限要素モデルでは、オメガを単体（2D では三角形、3D ではピラミッド）に分割します。3D の例として、EIDORS プロジェクト11の小さな非導電性ボールを備えた円筒形の液体充填タンクを取り上げます。このモデルは、内部の導電率分布を画像化するためにこのようなタンクに電流パターンを適用する電気インピーダンス トモグラフィーのアプリケーションを反映するように設計されています。FEM ジオメトリを記述するために、頂点nodesと単体のマトリックスを使用しますelems。
次の例では、反対側に 10 V と 20 V が課せられた (ディリクレ境界条件) 単純な長方形の 2D 導電性媒体を作成します。他のすべてのエッジは電気的に分離されています。
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); # No. of simplices
N = size (nodes,1); # No. of vertices
D = size (elems,2); # dimensions+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 ...
一次 FEM を使用して、オメガの電気伝導率分布を各単体（ベクトル で表されるconductivity）上の定数として近似します。有限要素ジオメトリに基づいて、まず各単体（要素単位のシステム マトリックス の対角線上の 3 行 3 列の要素として表される）のシステム（または剛性）マトリックスを計算します。 および単体と頂点間の接続を表す N 行 DE 列の接続マトリックスSEに基づいて、グローバル接続マトリックスが計算されます。 SECS
## Element conductivity
conductivity = [1*ones(1,16), ...
2*ones(1,48), 1*ones(1,16)];
## Connectivity matrix
C = sparse ((1:D*E), reshape (elems', ...
D*E, 1), 1, D*E, N);
## Calculate system matrix
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
## Element-wise system matrix
SE = sparse(Siidx,Sjidx,Sdata);
## Global system matrix
S = C'* SE *C;
システム行列は S オームの法則における 導電率のような働きをしますS * V = I。ディリクレ境界条件とノイマン境界条件に基づいて、各頂点の電圧を解くことができますV。
## Dirichlet boundary conditions
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; # vertices without Dirichlet
# boundary condns
idx(D_nodes) = [];
## Neumann boundary conditions. Note that
## N_value must be normalized by the
## boundary length and element conductivity
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";
Figure 22.6: 三角形要素を示す有限要素モデルの例。各頂点の高さは解の値に対応します。
Footnotes
(11)
EIDORS - Electrical Impedance Tomography and Diffuse optical Tomography Reconstruction Software http://eidors3d.sourceforge.net