連立一次方程式における既知量と未知量の取扱
ここで紹介しているFEM プログラムでは，剛性マトリックスの大きさを変えずに境界条件を導入し，連立方程式を解いている．この方法を用いる場合，連立一次方程式において，未知数と既知数の入れ替えが必要であるため，その考え方を以下に示す．
簡単な事例として，以下の連立一次方程式を考える．
\begin{align}
& k_{11} x_1 + k_{12} x_2 + k_{13} x_3 = f_1 \\
& k_{21} x_2 + k_{22} x_2 + k_{23} x_3 = f_2 \\
& k_{31} x_3 + k_{32} x_2 + k_{33} x_3 = f_3 \\
\end{align}
ここで，未知数が $x_1, x_3, f_2$，既知数が $f_1, f_3, x_2$だとすれば，行列演算を行うためには，以下のように書き換える必要がある．
\begin{align}
& k_{11} x_1 + 0 + k_{13} x_3 = f_1 - k_{12} x_2 \\
& k_{21} x_2 - f_2 + k_{23} x_3 = 0 - k_{22} x_2 \\
& k_{31} x_3 + 0 + k_{33} x_3 = f_3 - k_{32} x_2 \\
\end{align}
上記関係を，拡張して考えると，以下の通りとなる．
元の剛性方程式
\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & k_{1i} & \ldots & k_{1j} & \ldots & k_{1n} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{i1} & \ldots & k_{ii} & \ldots & k_{ij} & \ldots & k_{in} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{j1} & \ldots & k_{ji} & \ldots & k_{jj} & \ldots & k_{jn} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{n1} & \ldots & k_{ni} & \ldots & k_{nj} & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
\delta_1 \\
\vdots \\
\delta_i \\
\vdots \\
\delta_j \\
\vdots \\
\delta_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
f_i \\
\vdots \\
f_j \\
\vdots \\
f_n
\end{Bmatrix}
\end{equation}
境界条件導入後の剛性方程式
$\delta_i$ および $\delta_j$ が既知とすれば，剛性マトリックスの $k_{ii}$ および $k_{jj}$ を $1$ とし，それ以外の $i$ 列および $j$ 列の要素をゼロとする処理を行う．また剛性マトリックスにおける $i$ 列および $j$ 列の効果を右辺に移項する必要がある．
\begin{equation}
\begin{bmatrix}
k_{11} & \ldots & 0 & \ldots & 0 & \ldots & k_{1n} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{i1} & \ldots & 1 & \ldots & 0 & \ldots & k_{in} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{j1} & \ldots & 0 & \ldots & 1 & \ldots & k_{jn} \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
k_{n1} & \ldots & 0 & \ldots & 0 & \ldots & k_{nn}
\end{bmatrix}
\begin{Bmatrix}
\delta_1 \\
\vdots \\
-f_i \\
\vdots \\
-f_j \\
\vdots \\
\delta_n
\end{Bmatrix}
=\begin{Bmatrix}
f_1 \\
\vdots \\
0 \\
\vdots \\
0 \\
\vdots \\
f_n
\end{Bmatrix}
-\delta_i
\begin{Bmatrix}
k_{1i} \\
\vdots \\
k_{ii} \\
\vdots \\
k_{ji} \\
\vdots \\
k_{ni}
\end{Bmatrix}
-\delta_j
\begin{Bmatrix}
k_{1j} \\
\vdots \\
k_{ij} \\
\vdots \\
k_{jj} \\
\vdots \\
k_{nj}
\end{Bmatrix}
\end{equation}
この操作により，剛性マトリックスは非対称となりますが，数千自由度の連立一次方程式であれば，numpy の $x = np.linalg.solve(A, b)$ で大きなストレスなく解くことができる．
実際のプログラムでの処理
実際のプログラム（構造解析）での処理を以下に示す．連立方程式を解いた後，既知変位に相当する自由度の項に，既知変位を入れ直すことを忘れないように．
# treatment of boundary conditions
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
fp[iz]=0.0
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
fp=fp-rdis[j,i]*gk[:,iz]
gk[:,iz]=0.0
gk[iz,iz]=1.0
# solution of simultaneous linear equations
disg = np.linalg.solve(gk, fp)
# recovery of restricted displacements
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
disg[iz]=rdis[j,i]
| npoin | ：総節点数 |
| nfree | ：１節点自由度 |
| mpfix | ：各節点毎の拘束条件を格納した配列 |
| fp | ：節点外力ベクトル |
| rdis | ：変位拘束する節点の変位量 |
| gk | ：全体剛性マトリックス |
| disg | ：全体剛性方程式の解（変位） |