2.4節全体剛性の組み立て

　2.4節は「全体剛性マトリックス」である。Dマトリックス・Bマトリックスの作成と来たのでこの順番は妥当である。妥当なのは順番だけな気がするが。
　実際の全体剛性マトリックスの作成方法は次のようなステップを踏む。

全体剛性マトリックスKを0クリアする。
各要素剛性マトリックスの各要素を対応する節点自由度に足しこむ。

　要するにこれだけのことなのだが、図入りで説明しているところを見ると、こんなものでも理解できない連中がごっそりいるのかと思ってしまう。三好本でも同じことをしていたからイメージできない連中が多いということなのだろう。
　この部分だけなら単に面倒なだけで、何をするかさえ理解できていれば、難しいことをしているわけではない。
　この後に境界条件の組み込みがるのだが、実は全体剛性マトリックスの作成では、各要素剛性マトリックスの値を全体剛性に組み込むのと同時に境界条件を考慮することが出来るが、初心者にこのようなことを要求するのは酷なので、とにかく初めに全体剛性マトリックスを作成し、この後に境界条件の処理を行うようにしているのは仕方がないところである。

　要素剛性マトリックス$K_e$は角要素一つについての剛性(強度)を定義する。この要素剛性マトリックスでも適当な荷重と変位を作用させると次の式が成り立つ。

\mathbf f_e=\mathbf k_e\mathbf u_e

　各項は行列とベクトルなので、ボールド体にしている。二つのベクトルはそれぞれコネクティビティ順に節点番号が並んでいる。今回の場合、構造解析だから各節点番号と自由度の順に並べることになる。本書の要素1番の場合、X,Y軸方向の変位をu,vとした場合、$u_1,v_1,u_2,v_2,u_5,v_5$の順に並ぶことになる。添え字は節点番号である。今回は関係ないが、逆に変位成分を優先することもできる。しかし、全体剛性マトリックスの容量を削減するために、通常は節点番号優先で並べる。
　各要素剛性を別個に解いても正しい解は得られない。まず、各要素剛性方程式の同じ変位成分はどこも同じ値になる。節点1番は要素1と要素4で共有されているが、この節点変位は両方の要素で同じになる。というか同じにならないとおかしい。
　次に荷重である。荷重については節点変位と異なり、合計が0になる。外力をかけるとその外力をキャンセルするように荷重が各要素の当該節点にかかることになる。先の例でいうと、要素1と要素4は向きは置くとして外力が無ければ同じ値の荷重がかかる。要素4を半分にしたらどうなるだろうか。荷重の大きさの比は他所1と要素4(分割済み)で2:1:1になるだろう。実際には荷重の向きを考慮する必要があるから要素1と要素4の荷重の向きは逆になる。
　これが全体剛性マトリックスを作る際に要素剛性マトリックスの各成分を足しこむ理由である。
　最初に0クリアする理由は単純で、各要素剛性マトリックスの値を足しこむためで、最初は0とする必要があるためである。なんかの値が入っていたら、その部分の剛性がおかしくなるからである。

　全体剛性マトリックスの作成手順と詳細説明を端的に行うと以上のようになる。

　この後プログラミングについて説明をしているが、剰余演算子を用いるのはどうかなとは思う。そのようなものを使わない場合、自由度のループを追加する必要が生じるから、一長一短ではある。参考までに自作のVBA版では各節点において、全自由度の位置を力業でいれているから、この辺りは好みと効率の問題になる。ただし自由度の数が変わる場合この手は使えないから、コネクティビティの数と自由度でそれぞれループを回した方がいいだろう。

2.5節境界条件

　まず「x成分とy成分が交互に並んでいる」は「各節点ごとに各自由度(座標成分)が並ぶ」でいいだろう。今回の連では自由度が2だから強ち間違いとは言えないが、自由度数がこれ以上になると正確さを欠く記述は避けたほうがいい。
　続くP68は変数U,Fの処置であるが、ここからわけがわからない説明が出てくる。当プログラムもこの辺は踏襲しているがここには設定した値がそのまま入れている。なお、これらの配列は最初に0にセットして計算時に変な値を出さないようにするためである。伊達や酔狂で取り敢えず0にしているわけではない。
　通常変位を0に設定する自由度を拘束自由度というが、0でない変位を入れる項を強制変位などと言っている。この本では書いていないが、それ以外の自由度は未知変位と言って、計算で求める量ということになる。
　同じことは荷重にも言えるが、荷重の場合殆どの項が0、つまり値を設定しないから特に呼び方は無かったりする。ただ、自由度を拘束した場合、外荷重をかけても無視されることになる。その代わり、荷重項には反力が入ることになる。この本では何故かこのことが省かれているが、反力も未知だから未知変位諸共求めることが出来る。。
　　この辺りは全体剛性マトリックスの計算との兼ね合いもあって少し難しいのだが、要するに、実際のプログラミングでは既知項を一つのベクトルにまとめたあと、消去法などを用いて方程式を解くことになる。
　FEMでは連立一次方程式を解く際、後述するように拘束変位や強制変位の処理を行う必要がある。このため、全剛性方程式を作ってそれでおしまいとはならないのである。この処理を行うには、拘束変位項を覚えておかなければならない。そうしないと計算そのものに支障をきたす。本書のUmはそのために用意されている。
　ここでは反力を含めた計算方法を説明する。本書にならって3×3の方程式を解くことを考える。このとき2番目の変位$u_2$が拘束されているとすると、式は次のようになる。

\left[ \matrix{
k_{11}&k_{12}&k_{13}\\
k_{21}&k_{22}&k_{23}\\
k_{31}&k_{32}&k_{33}\\
}\right]\left[ \matrix{
u_1\\u_2\\u_3
}\right]=\left[ \matrix{
f_1\\R_2\\f_3
}\right]

　2行目の荷重項が$R_2$となってるが、先ほど説明したとおり、この項には反力が入るため区別している。
　この式を展開すると次のようになる。

k_{11}u_1+k_{12}u_2+k_{13}u_3=f_1 \\
k_{21}u_1+k_{22}u_2+k_{23}u_3=R_2 \\
k_{31}u_1+k_{32}u_2+k_{33}u_3=f_3 \\

$u_2$のかかる項は基地だから右辺に移すことが出来る。

k_{11}u_1+k_{13}u_3=f_1-k_{12}u_2 \\
k_{21}u_1+k_{23}u_3=R_2-k_{22}u_2 \\
k_{31}u_1+k_{33}u_3=f_3-k_{32}u_2 \\

　次に$R_2$は未知だから左辺に移す。

k_{11}u_1+k_{13}u_3=f_1-k_{12}u_2 \\
k_{21}u_1+k_{23}u_3-R_2=-k_{22}u_2 \\
k_{31}u_1+k_{33}u_3=f_3-k_{32}u_2 \\

　これを行列表示すると次のようになる。

\left[ \matrix{
k_{11}&0&k_{13}\\
k_{21}&-1&k_{23}\\
k_{31}&0&k_{33}\\
}\right]\left[ \matrix{
u_1\\R_2\\u_3
}\right]=\left[ \matrix{
f_1-k_{12}u_2\\-k_{22}u_2\\f_3-k_{32}u_2
}\right]

これを一般の行列式で表すと次のようになる。

\left[ \matrix{
\mathbf K_{11} & \mathbf K_{12} \\
\mathbf K_{21} & \mathbf K_{22} \\
} \right]
\left[ \matrix{
\mathbf u \\
\bar{\mathbf u}
} \right]=
\left[ \matrix{
\mathbf F \\ \mathbf R
} \right]

　ここで$\mathbf K_{ij}$は自由度の拘束の有無で分けられた部分行列である。$\mathbf u$は未知変位、$\bar{\mathbf u}$は既知変位とする。また、$\mathbf F,\mathbf R$はそれぞれ外荷重と反力である。
　この式を展開すると次のようになる。

\mathbf K_{11} \mathbf u + \mathbf K_{12} \bar{\mathbf u} = \mathbf F \\
\mathbf K_{21} \mathbf u + \mathbf K_{22} \bar{\mathbf u} = \mathbf R

　$\mathbf K_{i2}\bar{\mathbf u}$は計算できるから、右辺に移動する。

\mathbf K_{11} \mathbf u = \mathbf F - \mathbf K_{12} \bar{\mathbf u}\\
\mathbf K_{21} \mathbf u = \mathbf R - \mathbf K_{22} \bar{\mathbf u}

　一方$\mathbf R$は未知項なので左辺に移動する。

\mathbf K_{11} \mathbf u = \mathbf F - \mathbf K_{12} \bar{\mathbf u}\\
\mathbf K_{21} \mathbf u - \mathbf R = - \mathbf K_{22} \bar{\mathbf u}

　これを再度行列表示すると次のようになる。

\left[ \matrix{
\mathbf K_{11} & \mathbf 0 \\
\mathbf K_{21} & -\mathbf I \\
} \right]
\left[ \matrix{
\mathbf u \\
\mathbf R
} \right]=
\left[ \matrix{
\mathbf F - \mathbf K_{12} \bar{\mathbf u}\\ 
-\mathbf K_{22} \bar{\mathbf u}
} \right]

　ここで$\mathbf 0$は0行列で、$\mathbf I$は単位行列である。以下、この式をそのまま解けば未知変位と反力をそのまま求めることが出来る。

　と、こんな感じで書いてほしかったよ。意図が分からないのだが、本文にあるような変位項をそのままにしているという訳の分からないことをするより、反力を同時に求めたほうがいい。実際には後半の説明にある$\mathbf K_{11}$のみ保持すれば未知変位項を求めることが出来る。また全体剛性マトリックスの対称性を崩すこともなくなる。反力は後述するように等価節点力から求めることが出来る。これをやると初心者向けではなくなるから、この点についてはコラムで触れる方がいいだろう。
　コラムでスパースマトリックスについて説明しているが、そのことを生かすとメモリ容量を抑えることが出来るくらいは書くべきだっただろう。初めのうちはフルマトリックスを作り、それで問題を解いても構わないが、そこで止まられても困る。それゆえにこの程度のことを書いておかないと、あとで苦労することは目に見えている。
　このようなことを実現するためには、各行について起点と終点の位置を保持するための変数が必要になるのと、全体剛性マトリックスを作成する際、拘束をかける自由度のある行は作成しないため、プログラムが複雑になることである。だからプログラムに慣れない最初のうちはフルマトリックスで問題を解き、それから拘束条件のある自由度を抜く方法を考えるといい。それから、各行について必要な領域を保持してメモリ容量を減らすといいだろう。この状態をバンドマトリックスという。これでも0の項が残るが、それをさらに減らすにはどうしたらいいか考えてみるのもいい。構造問題では全体剛性マトリックスは対称になるから上半分だけ記憶しておけば計算時間はざっと半分に減ることになる。バンドマトリックスから上半分だけの計算に変えるだけなら実は思ったほど手間はかからない。ピボットから下の行の要素がなくなるが、対称マトリックスだることから、ピボットと同じ行の要素を用いることで計算が出来ることになる。
　できるはずのことを書かないで、奇怪な説明を行うべきではないし、現実的にどうするかの入り口位は書いておく方がよかっただろう。

　前にも書いたが、このプログラムでは、データーを直接コードにしている。このような方法をハードコードというが、今回はその是非は置いておく。メッシュデーターがハードコードされているから、境界条件もそうなっている。この部分のコードだけを示す。

    ' 全体変位ベクトル
    U(1) = 0#: U(2) = 0#: U(3) = 0#: U(4) = 0#: U(5) = 0#: U(6) = 0#
    U(7) = 0#: U(8) = 0#: U(9) = 0#: U(10) = 0#: U(11) = 0#: U(12) = 0#

    ' 全体荷重ベクトル
    F(1) = 0#: F(2) = 0#: F(3) = 0#: F(4) = 0#: F(5) = 0#: F(6) = 0#
    F(7) = 0#: F(8) = -100#: F(9) = 0#: F(10) = 0#: F(11) = 0#: F(12) = 0#

    ' 拘束目印ベクトル
    Um(1) = True: Um(2) = True: Um(3) = False: Um(4) = False
    Um(5) = False: Um(6) = False: Um(7) = False: Um(8) = False
    Um(9) = False: Um(10) = False: Um(11) = True: Um(12) = False

　この部分の何がまずいかは前にも書いた。よって改良案だけ書く。要するにデーター設定用のサブルーチンを書き、当該コードをデーターを設定するサブルーチンの呼び出しに置き換えるだけである。まず値を設定するルーチンを書く。

`拘束条件の設定
sub set_SPC(byval node as integer, byval dof as integer, optional byval value as double = 0)
    dim i as integer
    i=(node-1)*2+dof
    U(i)=value
    Um(i)=true
end sub

'荷重条件の設定
sub set_pload(byval node as integer, byval fx as double, byval fy as double)
    F(node*2-1)=fx
    F(node*2)=fy
end sub

　これらのサブルーチンはそれぞれ拘束条件と荷重条件を与える。set_SPCには節点番号と自由度番号を与える。変位は任意の値を与えることが出来るが、省略することもでき、その時は0になる。前回書くのを忘れたが、パラメーターの定義で、byvalは値のコピーを渡す。こうしておけば、サブルーチン内で変数を書き換えても呼び出し元の変数が書き換わることはない。またoptioinalはパラメーターを省略するためのものである。通常拘束条件は定義通りに変位を0にするためこの設定でよい。
　なお、SPCとはSpecific Point Constraintの略で、NastranやLS-DYNAの入力データーで使われている。日本語に訳すと特定の点の拘束といった意味になる。
　set_ploadも基本的にset_SPCと同じだが、X,Y軸方向の荷重を与える点が異なる。
　実際の設定部分は次のようになる。

    dim i as integer
    for i=1 to NODES
        U(i)=0: Um(i)=false: F(i)=0
    next i

    set_spc 1,1
    set_spc 1,2
    set_spc 6,1

    setpload 4,0,-100

　あのウザいだけの部分はこの程度の工夫でこれだけ見やすくなる。

2.6節「連立方程式ソルバー」

　この部分について言うと「はじめての有限要素法プログラム」に挑戦するなら、これでもいいんじゃない、という所である。ただし最初のうちはだが。ここで紹介している方法を墨守していると笑われるので、もう少し実際的な方法について説明する。

先の連立一次方程式を例にとる。ただし$u_2$も未知としている。

\left[ \matrix{
k_{11}&k_{12}&k_{13}\\
k_{21}&k_{22}&k_{23}\\
k_{31}&k_{32}&k_{33}\\
}\right]\left[ \matrix{
u_1\\u_2\\u_3
}\right]=\left[ \matrix{
f_1\\f_2\\f_3
}\right]

　ガウスの消去法というのは対角項$k_{ii}$に着目して、その下の行のi列目の項を0にする操作を繰り返すことを言う。この$k_{ii}$の項のことをピボットという。各列は、このピボットとピボットの直下の対応する項の比をかけることになる。なお、ピボットが0だと方程式を解くことが出来ないので、この行の下にある列のうちピボット行が0でないものを探して交換するか、エラーを報告してプログラムを終了させるなどの対策を取る必要がある。
　次に0にした行については、各列の項と右辺の既知量に対しても0にした列と同様の操作を行う。これを最終行まで切り返すのである。この操作を前進消去という。
　次に最後の行は対角項だけだから、右辺の既知量をこの対角項で割れば最終行の未知項が求まる。であれば、最終列についてのみ代入操作をすれば、最終列を消去、つまり0にできる。この調子で最終列から2列目迄の項を消去していくと、最終的に一行の道量を求めることで全未知量の計算を終えることが出来る。この操作を後退代入という。
　文章で書くとこの様なものである。

　初めに1行目から処理を開始する。この場合、$k_{21},k_{31}$を0にするわけである。この書籍での説明では次のような方法を取っていた。なお、右上の添え字は演算回数を表す

\begin{eqnarray*}
&\frac{k_{21}}{k_{21}}-\frac{k_{11}}{k_{11}} &= 0 \\
&\frac{k_{22}}{k_{21}}-\frac{k_{12}}{k_{11}} &= k^1_{22} \\
&\frac{k_{23}}{k_{21}}-\frac{k_{13}}{k_{11}} &= k^1_{23} \\
&\frac{u_2}{k_{21}}-\frac{u_1}{k_{11}} &= u^1_2\\
\end{eqnarray*}

この式を次のように変形する。

\begin{eqnarray*}
&k_{21}-\frac{k_{11}}{k_{11}}{k_{21}} &= 0 \\
&k_{22}-\frac{k_{12}}{k_{11}}{k_{21}} &= k^1_{22} \\
&k_{23}-\frac{k_{13}}{k_{11}}{k_{21}} &= k^1_{23} \\
&f_2-\frac{f_1}{k_{11}}{k_{21}} &= f^1_2\\
\end{eqnarray*}

　右辺の記号に変化はないが、ここにはこの式を用いた計算結果が入る。
　3行目も同じようになる。

\begin{eqnarray*}
&k_{31}-\frac{k_{11}}{k_{11}}{k_{31}} &= 0 \\
&k_{32}-\frac{k_{12}}{k_{11}}{k_{31}} &= k^1_{32} \\
&k_{33}-\frac{k_{13}}{k_{11}}{k_{31}} &= k^1_{33} \\
&f_3-\frac{f_1}{k_{11}}{k_{21}} &= f^1_3\\
\end{eqnarray*}

これを行列表示すると次のようになる。

\left[ \matrix{
k_{11}&k_{12}&k_{13} \\
k_{21}-\frac{k_{11}}{k_{11}}{k_{21}} & k_{22}-\frac{k_{12}}{k_{11}}{k_{21}} & k_{23}-\frac{k_{13}}{k_{11}}{k_{21}} \\
k_{31}-\frac{k_{11}}{k_{11}}{k_{31}} & k_{32}-\frac{k_{12}}{k_{11}}{k_{31}} & k_{33}-\frac{k_{13}}{k_{11}}{k_{31}}
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f_1 \\ f_2-\frac{f_1}{k_{11}}{k_{21}} \\ f_3-\frac{f_1}{k_{11}}{k_{21}}
}\right]

記号を置き換える。

\left[ \matrix{
k_{11}&k_{12}&k_{13} \\
0 & k^1_{22} & k_{23}^1 \\
0 & k^1_{32} & k_{33}^1
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f_1 \\ f_2^1 \\ f_3^1
}\right]

　次に$k^1_{32}$を0にする。これも先と同じ処理を行えばいい。

\begin{eqnarray*}
&k^1_{32}-\frac{k^1_{22}}{k^1_{22}}{k^1_{32}} &= 0 \\
&k^1_{33}-\frac{k^1_{23}}{k^1_{22}}{k^1_{32}} &= k^2_{33} \\
&f^1_3-\frac{f^1_2}{k^1_{22}}{k^1_{32}} &= f^2_3
\end{eqnarray*}

　これらの式をマトリックスに反映すると次のようになる。

\left[ \matrix{
k_{11}&k_{12}&k_{13} \\
0 & k^1_{22} & k_{23}^1 \\
0 & 0 & k_{33}^2
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f_1 \\ f_2^1 \\ f_3^2
}\right]

　$k_{ii}$をピボットした時のj行の一般項表記は次のようになる。添え字kは列番号である。

k_{jk}-\frac{k_{ik}}{k_{ii}}k_{ji}

　つまり、ピボット行より下は上の式に沿って行列内の値を更新する。$k=i$とすればこの項は0になる。これはとりもなおさず、ピボットから下の要素を0にすることを意味している。
　そして、既知項も以下の式で更新する。

f_j-\frac{f_i}{k_{ii}}k_{ji}

　これを見ると分かるが、j行の各項に対して、ピボット行の対応する列要素に$k_{ji}/k_{ii}$を掛けた値を引けばよい。これを各列に対して行う。
　以上が前進代入と言われる作業である。効率化を考えるならば、どうせ0にするのだから対角項から右だけ計算すればよい。さらに対称行列の場合、左下は全部省略できる。この場合、ピボットから下の列は常に対称になることが確認されているので、そのまま計算が可能である。

　次に後退代入と呼ばれる作業を行う。最後の行に着目すると、

k^2_{33}u_3=f^2_3

となっているので、$u_3$を簡単に求めることが出来る。

u_3=f^2_3/k^2_{33}

　マトリックスに反映すると次のようになる。

\left[ \matrix{
k_{11}&k_{12}&k_{13} \\
0 & k^1_{22} & k_{23}^1 \\
0 & 0 & 1
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f_1 \\ f_2^1 \\ \frac{f_3^2}{k_{33}^2}
}\right]

　これを一度展開する。

\begin{eqnarray*}
k_{11}u_1+k_{12}u_2+k_{13}u_3&=&f_1 \\
k^1_{22}u_2+k_{23}^1u_3&=&f^1_2 \\
u_3&=&\frac{f_3^2}{k_{33}^2}
\end{eqnarray*}

　$u_3$が求まっているので、これを上の2式に代入する。

\begin{eqnarray*}
k_{11}u_1+k_{12}u_2+k_{13}\frac{f_3^2}{k_{33}^2}&=&f_1 \\
k^1_{22}u_2+k_{23}^1\frac{f_3^2}{k_{33}^2}&=&f^1_2 \\
u_3&=&\frac{f_3^2}{k_{33}^2}
\end{eqnarray*}

　更に既知となった項を右辺に移す。

\begin{eqnarray*}
k_{11}u_1+k_{12}u_2&=&f_1-k_{13}\frac{f_3^2}{k_{33}^2} \\
k^1_{22}u_2&=&f^1_2-k_{23}^1\frac{f_3^2}{k_{33}^2} \\
u_3&=&\frac{f_3^2}{k_{33}^2}
\end{eqnarray*}

　記号を置き換える。

\begin{eqnarray*}
k_{11}u_1+k_{12}u_2&=&f^3_1 \\
k^1_{22}u_2&=&f^3_2 \\
u_3&=&\frac{f_3^2}{k_{33}^2}
\end{eqnarray*}

　再度マトリックス表示をすると次のようになる。

\left[ \matrix{
k_{11}&k_{12}& 0 \\
0 & k^1_{22} & 0 \\
0 & 0 & 1
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f^3_1 \\ f^3_2 \\  \frac{f_3^2}{k_{33}^2}
}\right]

　そうすると第2列も第3列と同じように$u_2$を求めることが出来る。

u_2=\frac{f^3_2}{k^1_{22}}

　これもマトリックスに反映する。

\left[ \matrix{
k_{11}&k_{12}& 0 \\
0 & 1 & 0 \\
0 & 0 & 1
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
f^3_1 \\ \frac{f^3_2}{k^1_{22}} \\  \frac{f_3^2}{k_{33}^2}
}\right]

　残るは1行目でこれも現在はこのようになる。

k_{11}u_1+k_{12}u_2=f^3_1

　$u_2$を代入し、$u_1$について解く。

u_1=\frac{f^3_1-k_{12}u_2}{k_{11}}

　最終的なマトリックスは次のようになる。

\left[ \matrix{
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
}\right]
\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
\frac{f^3_1-k_{12}u_2}{k_{11}} \\ \frac{f^3_2}{k^1_{22}} \\  \frac{f_3^2}{k_{33}^2}
}\right]

　右辺第一項は単位行列に他ならないから省略できるので、最終的に次のようになる。

\left[ \matrix{
u_1\\u_2\\u_3
}\right]=
\left[ \matrix{
\frac{f^3_1-k_{12}u_2}{k_{11}} \\ \frac{f^3_2}{k^1_{22}} \\  \frac{f_3^2}{k_{33}^2}
}\right]

　これで後退代入が出来たことになる。
　添え字表記をするとこうなる。

　対角項のみとなった場合、以下の式を計算して変位を求める。

u_i=\frac{f_i}{k_{ii}}

　それ以外の項は以下の式を計算した結果を格納する。

u_j=u_j-k_{ji}u_i

　実際のプログラムでは対角項より右は計算に使用するのは、説明の通りだが、いちいち0クリアするのは時間を余計に消費するだけだったりするので、計算に用いるだけで、消去しない。

　最後に、元のプログラムはこれ。数学が苦手な人にはこの手はいいかもしれない。最もこの手が使えるのは最初のうちだけだが。このプログラムを見ると本文にある通りピボットが非常に小さな値を取った時の省略をしていないが、エラーを報告してプログラムから抜けとけよと言いたくなる。

Sub solve()
    Dim r As Integer, c As Integer      ' マトリックスの行と列のカウンター
    Dim rr As Integer, cc As Integer    ' マトリックスの行と列のカウンター
    Dim pivot As Double                 ' マトリックスの対角成分
    Dim p As Double                     ' 計算に使用するマトリックスの成分
    
    For r = 1 To DOF_TOTAL              ' 前進消去
        pivot = Kc(r, r)                ' 対角成分をpivotに代入
        
        For c = r To DOF_TOTAL          ' r行目の処理
            Kc(r, c) = Kc(r, c) / pivot
        Next c
        F(r) = F(r) / pivot
        
        For rr = r + 1 To DOF_TOTAL     ' r+1行目以下の処理
            p = Kc(rr, r)
            For cc = r To DOF_TOTAL
                Kc(rr, cc) = Kc(rr, cc) - p * Kc(r, cc)
            Next cc
            F(rr) = F(rr) - p * F(r)
       Next rr
    Next r

    For r = DOF_TOTAL To 1 Step -1      ' 後退代入
        U(r) = F(r)
        For c = r + 1 To DOF_TOTAL
            U(r) = U(r) - Kc(r, c) * U(c)
        Next c
    Next r
End Sub

　自作品からの抜粋。デフォルトパラメーターにピボットの許容値(1e-10)を入れており、ピボットがこれ未満になると計算が停止する。尚、エラーメッセージにMsgBoxを使っているのはマクロ上で動かすことを意図していたためである。上と同じく、一度コピーを取ってから計算を行うようにしている。なお、ピボットの下の行の処理で$k_ji$が0に近いと処理を飛ばすようにしている。元の式を見ればわかるが、ここが0だと何も引きようがなくなるので処理を省略することが出来る。そのほかはほぼ元の式そのままである。

Sub solve(Optional ByVal tol As Double = 0.0000000001)
    Dim pivot As Double, kij As Double
    Dim i As Integer, j As Integer, k As Integer

    copyK = cp_mat2(sysK)
    '前進消去
    For i = 1 To ntnode * 2
        pivot = copyK(i, i)
        If Abs(pivot) < tol Then
            MsgBox "エラー:" & Str(i) & "行のピボットが許容値未満になりました"
            End
        End If
        For j = i + 1 To ntnode * 2
            kij = copyK(j, i)
            If Abs(kij) >= tol Then
                For k = 1 To ntnode * 2
                    copyK(j, k) = copyK(j, k) - copyK(i, k) / pivot * kij
                Next k
                rhs(j) = rhs(j) - rhs(i) / pivot * kij
            End If
        Next j
    Next i
    '後退代入
    For i = ntnode * 2 To 1 Step -1
        pivot = copyK(i, i)
        rhs(i) = rhs(i) / pivot
        For j = 1 To i - 1
            rhs(j) = rhs(j) - rhs(i) * copyK(j, i)
        Next j
    Next i
End Sub

訂正 20190908 文面修正