このサイトで紹介している2次元のFEMプログラムの掛算と割算にカウンターを付けて、計算が終了するまでに 掛算数と割算数を数えて見ました。すると、連立方程式の計算は全体のほぼ80%でした。つまり計算時間を 短縮するには連立方程式をもっと効率的に解く方法を考えなければなりません。FEMの計算では、同じ剛性 マトリクス[K]を用い様々な荷重で変位や応力を計算することが有ります。ここでは、そのような計算に 適した方法であるLU Decomposition法を紹介します。この方法は流体の計算でも使うことができます。

■LU Decomposition法の原理■
通常私たちは連立方程式を[A]{x}={B}で表します。ここに{x}は未知数です。 もし[A]を[L][D][L]^Tに分解します。すると、連立方程式の計算を2つのステップに分解できます。 三角行列[L]とは、その対角要素は1で、対角要素より上の要素にはゼロが入っているマトリクスです。

とりあえず、何らかの方法で[L]が計算出来たとします。[L]の計算方法については後で紹介しますのでご安心下さい。 したがって、計算しなければならない連立方程式は[L][D][L]^T{x}={B}になります。

まず以下のように定義します。

すると[L][D][L]^T{x}={B} は、いかのように書き表すことができますね。

すると{x'}は代入法で計算できます。なぜかと言うと、x'_n=B_nになります。L_nnは1であることに 注意して下さい。

次に以下のように定義します。

すると、[D][L]^T{x}={x'}は[D]{x"}={x'}になりますから、{x"}は[D]^-1{x'}で間単に計算できます。 [D]^-1{x'}は[D]の各対角要素の逆数で計算できます。つまり、1/D_iです。

そして最後に[L]^T{x}={x"}を代入法で計算すると、目的の{x}が計算されます。

こんどは正方マトリクス[A]を用いて、上で説明した手順をFortranで書いてみましょう。