最近のエンジニアリング分野では、市販のソフトウェアを利用することが定着したゆえか、はたまた技術計算プログラムの新規開発があまりない時代背景ゆえか、数値解析のことを詳細に記した新刊本がほとんど出ていないようだ。書店の店頭に並ぶその手の本といえば、昔執筆された書籍が店頭に並んでいるだけと感じる。工学問題の中で非常に重要な問題の一つである固有値解析の算法(アルゴリズム)を勉強しようとすれば、まずこれらの参考書を手にすることになるだろうが、これらの本で紹介されている固有値算法は、たいていが比較的小規模サイズのマトリックスを対象とするものである。エンジニアが遭遇する大規模サイズのマトリックスを対象とする実用的固有値算法を説く最近の参考書は皆無と言っていいのではないだろうか。

1960～80年代には、当時活躍されていた数値解析の大家戸川隼人先生が、非常に有用な数値解析の出版物を何冊か出されていて、若輩のエンジニアであった筆者などは、大いに勉強させられたものである。また、当時の工学部応用力学系の研究室では、アプリケーションプログラムの開発が多くなされており、数値解析部門が活況を呈していた時代でもあった。

もし今、若い研究者もしくはエンジニアの人たちの中で、固有値解析のプログラムを開発しようという奇特な人がいたら、ちょっと困ることになるのではないかと想像する。と言うのは、例えば、それが有限要素法プログラムであるならば(大規模サイズのマトリックスを対象とするという意味で)、公開されたよほど豊富な数学ライブラリー(MathLib)を利用できる環境を持つ人でもない限り、固有値算法のプログラムは自分で作成する必要があるからである。一方、一般ユーザーが利用できるMathLibで用意されている固有値算法というのは、有限要素法プログラムなどのアプリケーションプログラムでは全面的には利用し難いものである。その理由は後述する。

したがって、現役プログラマーがアプリケーションプログラムの開発に固有値算法を必要とする場合、もはや絶版となっている昔の参考書を専門系図書館で探すことになるのだろうか。

筆者は、若い頃から何度か固有値解析に接触してきた。構造解析分野での座屈解析に始まり、いくつかの振動解析で固有値解析プログラムを開発してきた。また、通信分野での固有値解析も経験したことがある。そこで、以下、若い数値解析プログラマーあるいは、少しは固有値解析の性格を知りたいという若きエンジニアへの遺言のつもりで、実務レベルの視点で見た固有値解析の雑談をしてみたいと思う。実用的固有値算法の勉強への入門ガイドとしてお読みいただければ幸いに思う。ただ、固有値解析の理論を展開したり、話の都合上、固有値算法の要点を触れることはあっても詳細を記述するわけでないことだけは予め断っておく。それらについては、それぞれ適当な数学書をお読みいただきたい。

 

工学問題での固有値解析と言えば、振動解析(モード解析)が真っ先に頭に浮かぶことに異論はないであろう。工学が扱う数ある振動現象の中でも、弾性体の振動が一番の伝統的問題であることも間違いないであろう。弾性力学とは全く離れた異分野の振動解析においても、“質量マトリックス”などの用語が慣習的に使用されている例をみても、この分野での歴史を物語っているかと思う。
固有値解析を誘導する振動解析の方程式といえば、工学部生が習う微分方程式の中でも一番有名な時間変数に関する二階常微分方程式(1自由度の場合)であろう。

式中のUは言わずと知れた変位振幅ベクトルを表し、 Üが加速度を意味することなどはもちろんである。

もし式(1)右辺に強制外力項が存在する場合の応答解析であれば、われわれが大学教養課程で習う微分方程式論での一般解、特解の重ねあわせで応答解が求まることになる。しかし、右辺がゼロである場合のモード解析はそう簡単ではない。

式(1)のタイプの固有値解析は、ときおりMCK型の固有値解析と呼ばれるが、このタイプの固有値、固有ベクトルを求める場合、非常に苦労するはずである。この問題を厄介にしているのは、減衰マトリックスCの存在である。Cの存在を厳密に考慮した場合、固有値解析の解は複素数になってしまう。もちろん、Cの存在がなくても、剛性マトリックスKあるいは質量マトリックスMにそもそも複素数を含む場合のモード解析もあり、この場合も固有値、固有ベクトルが複素数になる。どちらにしても、解が複素数になる固有値解析では、固有値算法の選択に制約を受けることになる。

もし、全ての複素固有値を求めるというのであれば、少しは救われるだろうか。QR法という算法がこの解決に役に立つことになるからである。ところが、モード解析に出てくる各マトリックスのサイズNが非常に大きくなる場合は、計算機のメモリと計算時間に課題が生じてくることになる。というのも、QR法は全て(N個)の固有値を求める算法だからである。

一方、求める固有値の数がNに比べて非常に少ない場合、例えば最小(最大)固有値からp個(p≦N)を求める 場合 – 通常の工学問題のほとんどがこのタイプである– この問題は超難問となる。おそらく数学者の助けなしには解決しないだろうと思う。以前、筆者もこの問題にチャレンジしかかったことがある。一度は相談相手を探す段階までいったことがある経験から言うと、日本の応用数学者で、この問題の相談に乗ってくれる人は数少ないと想像する。幸か不幸か、筆者の場合、そのプロジェクトは中断となり、難を避けることができた(笑)。

ところで、減衰マトリックス c≠0 の場合でも、複素固有値を避けられる場合もある。次の関係式がある場合である。

この減衰式はレイリー減衰と呼ばれている。α、βは比例定数である。レイリー減衰は、物理的要請から創出されたものでなく、計算上の都合から生まれたものではないだろうか。たしかに、この関係式が成り立つなら、MCK型の固有値解析はMK型の固有値解析に還元されてしまう。もっともレイリー減衰は、応答解析では多用されるが、モード解析においては、そもそも c=0 で始めることがほとんどのようだ。

ここからは、 Mü+KU=0 であるMK型の固有値解析に限定して話を進めることにする。また、弾性体の振動解析に限らない固有値解析を鮮明にさせるためにも、記号も改めて、固有値解析の教科書内でよく表現される次の一般固有値解析の方程式に書き換える。

蛇足とはなるが、マトリックスAは弾性体振動解析の場合の剛性マトリックスKに相当し、マトリックスBが質量マトリックスMに相当することは言うまでもない。また、ここでの固有値算法を紹介するに当たって、言っておかなければならない前提がある。

• まず、マトリックスA、Bは対称マトリックスであること。これによって、固有値、固有ベクトルが実数であることが保証される。
• 固有値算法を論じるときは、マトリックスA、Bの正定値性、非特異性のことも確認しておく必要もあるのだが、ここでは、固有値解析の理論を述べているのではなく、あくまでも実用的算法を話題にしているので、それらの小難しい話題は避けることにする。
  もっとも、きちっと拘束された弾性体(構造物)の振動をイメージする通常のモード解析であれば、こういった条件は自然にクリアーしているので、気にすることはない。
• マトリックスA、BのサイズNが大きくて、しかもリクエストする固有値の数pとの関係は、あくまでもp«Nの場合である。例えば、N=100,000、p=10といった、有限要素法解析で出現するような固有値解析である。

いよいよ固有値算法を話す段になったが、その前にMathLibを利用することを少し考えてみる。結論から言ってしまえば、マトリックスA、Bとも引数渡しにするような完全にMathLibに頼る手段はお勧めできない。固有値算法の一部でMathLib利用することは、大いに結構なことだと思うが、これについては、その都度言及するつもり。

一般にMathLib内で用意されている固有値算法は式(3)のような一般固有値解析ではなく、次の形式の標準固有値解析であることが多い。

これは、一般固有値解析の方程式は標準固有値解析のそれに変換することができるからである。早い話、式(3)両辺にBの逆マトリックスを前から掛ければ、原理的には標準固有値解析の形式になってしまう。しかし、工学問題のアプリケーションプログラムでは、そんなバカなことはしない。そんなことをすれば係数マトリックスの対称性を崩すことになったり、Bの逆マトリックスをとるにしても、そもそもBが特異マトリックスであることもあるわけだからだ。それよりもなによりも、プリ処理で標準固有値解析へ変換する算法は計算効率を考慮した手段で実施することが求められるものだ。後述するランチョス法などはその典型的な優秀算法です。

次にMathLib内で用意された固有値算法に、一般固有値解析用のものが用意されていたとする。この場合、アプリケーション側で、巨大サイズになるA、B二つのマトリックスを用意する必要となる。いかに圧縮化した保存手段を採っていても、有限要素法などで登場するA、Bマトリックスはそれでも巨大サイズとなる。巨大サイズのマトリックスを二つもメモリ展開することは、メモリ面で非常に問題となる。

一方、どんな固有値算法であろうと、二つのマトリックスのうち、Aの方は、それを対象に原理的には逆マトリックスを求める必要があり(実際は連立方程式の解式を実行)、マトリックスAはそのままの形で保存することが要求されるが、マトリックスBについてはそんな必要はない。固有値算法の中では、中間(作業)ベクトルとの積で必要とされるだけなので、何も巨大サイズのBのままで保存する必要はないのである。この点は、少し分かり難いと思うので、有限要素法での用語を使って説明する。

今、仮に作業ベクトルをWとする。するとBWは全体質量マトリックスBとWの積を意味する。ところで、マトリックスの線形性を利用すると、何も巨大マトリックスの積計算をしなくても、サイズの小さい要素質量マトリックスの段階で、Wとの積計算をしておき、その結果を全要素について集計する手段でもいいわけである。このアルゴリズムでは、マトリックスBを保存する必要がなくなり、メモリ面で有利になり、かつ計算効率においても、乗算回数の圧倒的減少が図られて非常に有利となる。こんなことができるのも、アプリケーションプログラムのスキーム内で固有値算法を組み込むゆえ出来ることであって、単純なMathLibのコール利用では出来ないワザなのである。

そんなこんなで、もしあなたがプロ級のプログラマーを自認される方なら、固有値解析を含むアプリケーションプログラムの開発する場合、(部分的なMathLibコールは利用するにしても)全面的なMathLib利用は放棄して、アプリケーション内に固有値算法を採り入れたプログラム開発をお勧めする次第でる。

さて、ここからいよいよ本論に入っていこうとしたのだが、前置きが思いの外長くなってしまった。そこで、続きは次回へとさせていただく。悪しからず。

2011年10月記