STEP５：線形代数方程式の解法

大型マトリクスを有する線形弾性方程式を解く際には、採用する解法によって計算時間に大きな影響を与えます。そのため、扱うマトリクスの性質にしたがって、それに適した解法で解くことが望ましいです。材料マトリクス（\(\boldsymbol{D}\,\)マトリクス）が線形の場合には、先に述べたように全体剛性マトリクスが対称型となるため、対角成分と上三角行列のみをストアするのが合理的となります。そのうえで伝統的にはGaussの消去法を用います。また、このとき可変バンド消去法が利用可能であり、それに応じた計算アルゴリズムとすることで、追加的に計算コストを下げられます。こうした解法は直接法といわれるものです。
さらに、コンピュータの演算処理速度の向上とともに、昨今では、直接法の10分の1の計算時間を実現するスパースなどの反復法によるソルバも普及しています。
ここでは、伝統的な直接法について整理しておきましょう。

Gaussの消去法

線形弾性方程式 \begin{equation} \boldsymbol{K}\vec{u}=\vec{F} \label{79} \end{equation} を解く際、Gaussの消去法は下三角行列を0にするような値\(\,\boldsymbol{L^{-1}_i}\,\)を両辺に左から掛けていきます。すなわち全体剛性マトリクスが対角行列になるとき、式\eqref{79}には \begin{equation} \boldsymbol{L^{-1}_n} \cdots \boldsymbol{L^{-1}_2}\boldsymbol{L^{-1}_1}\boldsymbol{K}\vec{u}=\boldsymbol{L^{-1}_n} \cdots \boldsymbol{L^{-1}_2}\boldsymbol{L^{-1}_1}\vec{F} \end{equation} という操作を行っていることに相当します。ただし\(\,\boldsymbol{L^{-1}_i}\,\)は、 \begin{eqnarray} L^{-1}_i= \begin{bmatrix} 1 & & & & \\ 0 & 1 & & & \\ 0 & l_{j,i} & & & \\ \vdots & & & \ddots& \\ 0 & -l_{n,i} & 0 & 0 & 1 \end{bmatrix} & \quad & l_{j,i}=\frac{K^{\prime}_{j,i}}{K^{\prime}_{i,i}} \end{eqnarray} で表される下三角マトリクスです。ここで、 \begin{equation} \boldsymbol{L^{-1}_n} \cdots \boldsymbol{L^{-1}_2}\boldsymbol{L^{-1}_1}\boldsymbol{K}=\boldsymbol{U} \end{equation} とおくと全体剛性マトリクスは、 \begin{equation} \boldsymbol{K}=\boldsymbol{L_1}\boldsymbol{L_2}\cdots\boldsymbol{L_n}\boldsymbol{U}=\boldsymbol{L}\boldsymbol{U} \end{equation} と表されます。それでは、具体的に\(\,\boldsymbol{L}\,\)と\(\,\boldsymbol{U}\,\)を求めてみましょう。まず\(\,\boldsymbol{L}\,\)について、 \begin{equation} \boldsymbol{K}= \begin{bmatrix} 1 & & & & & \\ L_{21} & 1 & & & & \\ \vdots & & \ddots & & & \\ L_{i1} & \cdots & L_{ij} & 1 & & \\ \vdots & & & & \ddots & \\ L_{n1} & L_{n2} & \cdots & & \cdots & 1 \end{bmatrix} \begin{bmatrix} U_{11} & U_{12} & \cdots & U_{1j} & \cdots & U_{1n} \\ & U_{22} & & U_{2j} & & U_{2n} \\ & & \ddots & \vdots & & \vdots \\ & & & U_{jj} & & \\ & & & & \ddots & \\ & & & & & U_{nn} \end{bmatrix} \notag \end{equation} ですから、その成分\(\,K_{ij} \, \)は \begin{eqnarray} K_{ij} &=& L_{i1}U_{1j}+L_{i2}U_{2j} + \cdots +L_{i\,j-1}U_{j-1 \, j}+L_{ij}U_{jj} &=& \sum^{j-1}_{m=1}L_{im}U_{mj}+L_{ij}U_{jj} \end{eqnarray} と表され、これにより\(\,L_{ij} \, \)を求めると、 \begin{eqnarray} L_{ij}=\frac{ \left( K_{ij}-\displaystyle \sum^{j-1}_{m=1}L_{im}U_{mj} \right)}{U_{jj}} & \quad & (i=2～n \, , \quad j=1～i-1) \end{eqnarray} のように帰納的に定まります。次に\(\,\boldsymbol{U}\,\)については、 \begin{equation} \boldsymbol{K}= \begin{bmatrix} 1 & & & & & \\ L_{21} & 1 & & & & \\ \vdots & & \ddots & & & \\ L_{i1} & \cdots & L_{jj-1 \, j} & 1 & & \\ \vdots & & & & \ddots & \\ L_{n1} & L_{n2} & \cdots & & \cdots & 1 \end{bmatrix} \begin{bmatrix} U_{11} & U_{12} & \cdots & U_{1j} & \cdots & U_{1n} \\ & U_{22} & & U_{2j} & & U_{2n} \\ & & \ddots & U_{ji} & & \vdots \\ & & & U_{ii} & & \\ & & & & \ddots & \vdots \\ & & & & & U_{nn} \end{bmatrix} \end{equation} と考えると、その成分\(\,K_{ij} \, \)は、 \begin{eqnarray} K_{ij}=L_{j1}U_{1i}+L_{j2}U_{2i} + \cdots +L_{jj-1}U_{j-1 \, i}+U_{ji}=\sum^{j-1}_{m=1}L_{jm}U_{mi}+U_{ji} \end{eqnarray} と表せるので、\(\,U_{ji} \,\)を求めると、 \begin{eqnarray} U_{ji} = K_{ji} - \sum^{j-1}_{m=1}L_{jm}U_{mi} & \quad & (i=1～n & \, , \quad & j=1～i) \end{eqnarray} のように帰納的に定まります。ここで、ともに下三角マトリクスである場合、その積は再び下三角マトリクスとなりますから、\(\boldsymbol{L}\)は下三角マトリクスとなっています。以上のことから式\eqref{79}は、 \begin{equation} \boldsymbol{LU}\vec{u}=\vec{F} \end{equation} ゆえに、 \begin{equation} \boldsymbol{U}\vec{u}=\boldsymbol{L^{-1}}\vec{F} \end{equation} となってますので、ガウスの消去法では、まず前進消去によって\(\vec{y}=\boldsymbol{L^{-1}}\vec{F}\)を求めた後に、後退代入 \begin{equation} \boldsymbol{U}\vec{u} = \vec{y} \end{equation} することで\(\, \vec{u}\,\)を求めています。

可変バンド消去法

全体系における剛性マトリクスは前節(2)で述べたように式\eqref{77}のように、①構造的に定まった位置に0成分を含む、②対称行列となる、ことがわかっています。この性質を利用して対角成分と上三角成分だけを記憶領域にとれば充分であり、計算時間も短縮できるため、STRUCT（インハウスプログラム）ではこれらの成分を次のように１次元配列でストアしています。 \begin{equation} \begin{bmatrix} G(1) & G(3) & & & G(13) & G(19) \\ & G(2) & G(5) & G(8) & G(12) & G(18) \\ & & G(4) & G(7) & G(11) & G(17) \\ & & & G(6) & G(10) & G(16) \\ & & & & G(9) & G(15) \\ & & & & & G(14) \end{bmatrix} \end{equation} これをLU分解しても形状は変わらずに、 \begin{equation} \begin{bmatrix} U(1) & U(3) & & & U(13) & U(19) \\ & U(2) & U(5) & U(8) & U(12) & U(18) \\ & & U(4) & U(7) & U(11) & U(17) \\ & & & U(6) & U(10) & U(16) \\ & & & & U(9) & U(15) \\ & & & & & U(14) \end{bmatrix} \end{equation} と置くことができます。

Step-by-step一覧へ戻る