本エッセイ第52話では、“マトリックス三重積”という話をした。そこでの話題は、 C^tBCという形の三つのマトリックス乗算の話だった。今回、さらに中央のBが逆マトリックス(インバース)である場合、すなわち、

という形のマトリックスの演算式を考えてみる。この形式のマトリックス積は有限要素法の世界では、ときおりお目にかかる。ここでは、これを俎上に載せてみたい。

数ある要素剛性マトリックスの定式化タイプの中には、その誘導過程の中で、所属要素内でのみで通用する仮の節点を要素内部に設けるものや、本エッセイ第1話で紹介したHCT要素のように原要素をさらに細分割したために、要素内部に新たにできる節点、というようなものがある。あるいは非適合要素の場合のような付加モードという直接、既存節点とは関係しない未知自由度が追加されることもある。これら、新たに生成された節点に属する自由度や付加モードは、最終的には、要素剛性マトリックスが全体剛性マトリックスへアセンブリされる前に消去される運命にある。この一連の作業を“静的縮退”という。

静的縮退は、未知変数(通常は、要素変位自由度)を目的にあった二つのグループに分けることから始める。ここでは、アセンブリ時に残す自由度群と縮退対象の自由度群の二つとなる。
次の式は静的縮退の説明するものである。ただし、説明に都合いいように、マトリックス、ベクトルをわざと分割した変数記述をしている。静的縮退のない通常の要素剛性マトリックスの場合では、 k₁₁u₁＝P₁ となるものだ。

式(2)のマトリックス式を上下2式に分解して、下式の方からU₂を消去し、それをもう一つの式に代入すると、次の式にまとめられることが容易にわかる。

ここで、 k₂₁=K₁₂^t であることを考慮して、式(3)を書き換えると、左辺括弧内第2項に式(1)と同形式の次のマトリックス3重積が出てくることがわかる。

さて、式(4)のマトリックス計算をプログラミングするには、マトリックスサイズが小さいものを対象とする間は、数式が示す通り素直にプログラミングすればよいかと思う。この形式が現れる場合、大抵K₂₂ のサイズが相対的に K₁₁のサイズより小さいうえ、要素剛性マトリックス自身さほど大きなマトリックスでないため、式(4)中央のマトリックスのインバース計算に気遣いすることもないと考えるからである。それで、この計算、まずK₂₂のインバースを求めておいてから、後ろからK₂₁を掛ける算段となるのが素直な手法である。

数値計算で注意しないといけないのは、マトリックスサイズが大きくなってきた場合である。マトリックスのインバース計算はその演算量がマトリックスサイズNの三乗に比例するという大変不利な数値計算であることがよく知られている。ところが、k₂₂^-1k₂₁の計算のところでは、実はもう一つの選択肢があるのである。

K₂₂の次元を(N×N)とすると、K₂₁ の次元は(N×M :Mは任意)となる。 K₂₁を要素数Nの縦ベクトルがM本並んでいると思えば、 k₂₂^-1k₂₁は連立方程式K₂₂X=K₂₁の縦ベクトル の求解にほかならないのだ。したがって、M本の右辺ベクトルの連立方程式を一挙に解いてしまえば、 k₂₂^-1k₂₁を求めたことになるわけである。連立方程式の解法での定番的手法、三角分解法では同マトリックスのインバース計算に比べて、演算量が数分の一で済むことが分かっている。

メモリの点でもインバース計算が不利なことが存在する。下の式は左側が元のマトリックスを右側がそのインバースを表している。

この例のように、三重対角マトリックスという最高レベルのスパースなバンドマトリックスであっても、そのインバースはフルマトリックスとなってしまうように、いかに元のマトリックスがスパース性を持っていても、そのインバースはもはやスパース性とは限らないのである。したがってマトリックスの記憶において、そのスパース性を利用したメモリ節約ができないことになる。これが、連立方程式の解式とみなして、係数マトリックスを三角分解する場合では、元マトリックスのバンド領域のみの値保持だけで済むのである。

さて、今までは、あくまでも要素剛性マトリックスのようにマトリックスサイズが小さめの世界での話であった。これが、全体剛性マトリックスレベルの大次元のマトリックスを扱う世界になると、さらに厄介な問題が待ち構えている。

振動解析では、大次元の係数マトリックスを保有する固有値方程式が現れる。固有値の数値解析は、連立方程式の解式に比べて格段に時間の掛かるものである。そこで、なるべく全体の自由度を減らす手段を考えようとして、ギヨンの縮退法といったものが採用されたりしている。この縮退は、先の要素剛性マトリックス段階での縮退ではなく、全体剛性マトリックス段階での縮退であることが、大きな違いである。ただし、この場合でも、式(3)のような形式の方程式が現れることは、全く同様である。

ところで、同じように式(3)形式で数式表現される要素剛性マトリックスと全体剛性マトリックスの両者の違いに、後者が式(3)左辺の括弧内全体のインバースを求める-実際には連立方程式を解く宿命があるということです。

そして、この時に大きな問題が立ちはだかるのである。一般的に、工学問題での大型係数マトリックスを持つ連立方程式の解式では、係数マトリックスのスパース性を利用した解法が採られることが常識のようになっている。式(3)を見つめると、左辺第1項はたしかに、スパース性を保持しているが、第2項の加算により、もはやスパース性は失われてしまうのである。したがって、式(3)形式の連立方程式を解くというアルゴリズムを迂闊にプログラムに採用できないことになってしまう。

 

せっかく数学的にかっこよく、式の展開を図ってきても、いざ数値解析の段になると、その数式が使えないという状況に陥ることは、別に珍しいことではなく、数値解析の世界ではよく経験することである。

2012年5月記