固有値解析¶

一般化固有値問題¶

連続体の自由振動解析を行う場合、空間的離散化を行い、図2.3.1に示すような集中質点による多自由度系でモデル化される。減衰のない自由振動問題の場合、支配方程式（運動方程式）は以下のとおりである。

\[\begin{equation} M \ddot{u} + K u = 0 \label{eq:2.3.1} \end{equation}\]

ただし、$u$ は一般化変位ベクトル、$M$ は質量マトリックス、$K$は剛性マトリックスである。ところで、固有角振動数を$\omega$とし、$a$、$b$を同時には0でない任意定数、$x$をベクトルとして、関数

\[\begin{equation} u(t) = (a \sin \omega t + b \cos \omega t ) x \label{eq:2.3.2} \end{equation}\]

を定義する。ここで、この式の2階の微分は、

\[\begin{equation} \ddot{u}(t) = -\omega^2 (a \sin \omega t + b \cos \omega t) x \label{eq:2.3.3} \end{equation}\]

となる。式$\eqref{eq:2.3.2}$および式$\eqref{eq:2.3.3}$を式$\eqref{eq:2.3.1}$に代入すれば、

\[\begin{equation} M \ddot{u} + K u = (a \sin \omega t + b \cos \omega t) (K-\omega^2 M) x = 0 \label{eq:2.3.4} \end{equation}\]

となる。

非自明な振動を考えると、$a \sin \omega t + b \cos \omega t$は恒等的には0ではないため、

\[ (K-\omega^2M)x=0 \]

を得る。したがって、$\lambda=\omega^2$とおけば、

\[\begin{equation} K x = \lambda M x \label{eq:2.3.5} \end{equation}\]

を得る。

係数$\lambda$を固有値、ベクトル$x$を固有ベクトルと呼び、式$\eqref{eq:2.3.5}$で表される問題を一般化固有値問題と呼ぶ。

固有値$\lambda=\omega^2$から固有角振動数$\omega$が得られ、対応する固有ベクトル$x$は振動モードを表す。

減衰のない自由振動の多自由度系の例

図 2.3.1　減衰のない自由振動の多自由度系の例

行列の性質と仮定¶

前節で得た一般化固有値問題 $Kx=\lambda Mx$ に対し、本マニュアルでは対象となる行列の対称性を仮定する。複素行列の場合にはエルミート行列、実行列の場合には対称行列に相当する。

行列$K$の$ij$成分を$k_{ij}$とした時、エルミート性は

\[\begin{equation} k_{ij} = \bar{k}_{ji} \label{eq:2.3.6} \end{equation}\]

と表される。ここで、$\bar{k}_{ji}$は$k_{ji}$の複素共役である。実行列の場合には、この関係は$k_{ij}=k_{ji}$となる。

また、実対称行列$H$が正定値であるとは、任意の非零ベクトル$x$に対して

\[\begin{equation} x^{t} H x > 0 \label{eq:2.3.7} \end{equation}\]

が成り立つことをいう。この場合、$H$の固有値はすべて正である。

構造固有値問題では、質量マトリックス$M$は通常正定値として扱われる。一方、剛性マトリックス$K$は拘束条件によっては半正定値となり、剛体モードに対応する零固有値を持つ場合がある。

シフト付逆反復法¶

有限要素法による構造解析では、実用上、全ての固有値は必要とせず、高々数個の低次の固有値で十分な場合が多い。ところで、HEC-MWでは大規模な問題を扱うことを想定しており、行列はサイズが大きく非常に疎（零要素が多い）である。したがって、この事を念頭に低次のモードの固有値を効率よく求めることが重要である。

シフト量を$\sigma$とした時、$-\sigma$が固有値と一致せず、$K+\sigma M$が正則であれば、式$\eqref{eq:2.3.5}$は次式のように変形できる。

\[\begin{equation} (K + \sigma M)^{-1} M x = \frac{1}{\lambda+\sigma} x \label{eq:2.3.8} \end{equation}\]

この変換では固有ベクトル$x$は変化せず、固有値$\lambda$は$1/(\lambda+\sigma)$に写される。

したがって、$\lambda$が$-\sigma$に近いほど、変換後の固有値の絶対値は大きくなる。構造固有値問題では$\lambda \geq 0$、$\sigma \geq 0$であるから、最も低次の固有値が絶対値最大の固有値に写る。この性質を利用し、絶対値の大きい固有値から収束しやすい反復法を式$\eqref{eq:2.3.8}$に適用することで、低次の固有値から順に効率よく求めることができる。

この手法をシフト付逆反復法と呼ぶ。

FrontISTRでは、拘束条件のある解析では$\sigma = 0$とし、式$\eqref{eq:2.3.8}$は$K^{-1} M x = \frac{1}{\lambda} x$、すなわちシフトなしの逆反復となる。拘束条件のないフリーフリー解析では$K$が剛体モードに対応する零固有値を持ち特異となるため、$\sigma$に正の値を与えて$K+\sigma M$を正則化する。このときの$\sigma$は!EIGENのSIGMAで指定する。

ランチョス法¶

採用理由（Jacobi 法との比較）¶

古典的な方法ではJacobi法がよく知られている。

この方法は、行列サイズが小さく密行列である時、有効である。しかしながら、HEC-MWで扱う行列は大規模で疎であるため、この方法は採用せずランチョス(Lanczos)反復解法を採用している。

アルゴリズムと特徴¶

1950年代にC. Lanczosにより提案されたこの手法は、行列を三重対角化する計算算法であり、下記のような特徴を有している。

反復収束解法であり、行列を疎のまま計算を進めることができる。
算法は行列、ベクトル積が中心となっており並列化に適している。
有限要素メッシュに伴う幾何学的領域分割法に適している。
求める固有値の個数やモード範囲を限定して効率よい計算を行える。

ランチョス法は、初期ベクトルからスタートして順次直交ベクトルを作成し、Krylov部分空間の基底を求める計算を行うものである。

この手法では、有限精度演算による丸め誤差の影響でベクトルの直交性が失われる場合がある。 FrontISTRの実装では、この影響を抑えるため、既に得られたLanczos基底に対して再直交化を行っている。

幾何学的意味（Krylov 部分空間）¶

式$\eqref{eq:2.3.8}$を次のように変数変換することにより

\[ A = (K + \sigma M)^{-1} M \]

\[\begin{equation} \frac{1}{\lambda+\sigma}= \zeta \label{eq:2.3.9} \end{equation}\]

問題を書き直すと

\[\begin{equation} A x = \zeta x \label{eq:2.3.10} \end{equation}\]

を得る。

適当な非零ベクトル$q_0$に対して、

\[ q_0,\ Aq_0,\ A^2q_0,\ldots,A^{m-1}q_0 \]

が張る空間

\[ \mathcal{K}_m(A,q_0) = \operatorname{span} \{q_0,Aq_0,A^2q_0,\ldots,A^{m-1}q_0\} \]

をKrylov部分空間と呼ぶ。

ランチョス法では、このKrylov部分空間の基底を逐次構成する。

FrontISTRでは、質量マトリックス$M$に関する内積

\[ \langle x,y\rangle_M = x^T M y \]

を用いて基底を正規直交化する。以下の図に示す内積$\langle x,y\rangle$およびノルム$\|x\|$は、 FrontISTRの計算ではそれぞれこの$M$-内積と、それに対応する$M$-ノルム

\[ \|x\|_M=\sqrt{x^T M x} \]

として解釈する。

適当なベクトル$q_0$に対して行列$A$による一次変換を行う（図2.3.2参照）。

$行列$A$による$q_0$の一次変換$

図 2.3.2　行列$A$による$q_0$の一次変換

変換されたベクトルは、元のベクトルとつくる空間の中で直交化される。すなわち、図2.3.3のようなグラム・シュミットの直交化を行う。そうして得られたベクトルを$r_1$として、それを正規化して$q_1$を得る。

$q_0\)に直交なベクトル\(q_1$

図 2.3.3　$q_0$に直交なベクトル$q_1$

同様な算法により$q_1$から$q_2$を得る。このとき$q_2$は$q_1$、$q_0$の両方に直交している（図2.3.4）。

$q_1\)と$q_0$に直交なベクトル\(q_2$

図 2.3.4　$q_1$と$q_0$に直交なベクトル$q_2$

このように、Lanczos法ではKrylov部分空間の正規直交基底を逐次構成する。理論上、対象となる固有値問題の対称性を利用することで、この反復は直近の基底ベクトルを用いる三項漸化式として表すことができる。

一方、FrontISTRの実装では、有限精度演算による直交性の喪失を抑えるため、既に得られたLanczos基底に対して$M$-内積による再直交化を行っている。

三重対角化¶

FrontISTRのLanczos反復では、前節で述べた$M$-内積に対して基底ベクトルを正規直交化するため、

\[ q_i^T M q_j = \delta_{ij} \]

が成り立つ。

式$\eqref{eq:2.3.10}$の行列$A$を用いると、理論上のLanczos反復は

\[\begin{equation} A q_i = \beta_i q_{i-1} + \alpha_i q_i + \beta_{i+1} q_{i+1} \label{eq:2.3.11} \end{equation}\]

と表せる。

まず、$\alpha_i$を

\[ \alpha_i = q_i^T M A q_i \]

で定め、

\[ \tilde{r}_{i+1} = Aq_i - \beta_i q_{i-1} - \alpha_i q_i \]

とする。

FrontISTRの実装では、有限精度演算による直交性の喪失を抑えるため、 $\tilde{r}_{i+1}$を既に得られたLanczos基底に対して$M$-内積で再直交化する。再直交化後の残差を$r_{i+1}$とすると、

\[\begin{equation} \beta_{i+1} = \sqrt{r_{i+1}^T M r_{i+1}}, \qquad q_{i+1} = \frac{r_{i+1}}{\beta_{i+1}} \label{eq:2.3.12} \end{equation}\]

となる。

Lanczos反復によって得られた$m$本の基底を

\[ Q_m=[q_0,q_1,\ldots,q_{m-1}] \]

とすると、有限回のLanczos反復では

\[\begin{equation} A Q_m = Q_m T_m + \beta_m q_m e_m^T \label{eq:2.3.13} \end{equation}\]

の関係が成り立つ。

ここで、$e_m$は第$m$成分のみが1である$m$次元の単位ベクトルであり、

\[\begin{equation} T_m= \begin{pmatrix} \alpha_0 & \beta_1 & & &\\ \beta_1 & \alpha_1 & \beta_2 & &\\ & \ddots & \ddots & \ddots &\\ & & \beta_{m-2} & \alpha_{m-2}& \beta_{m-1}\\ & & & \beta_{m-1} & \alpha_{m-1} \end{pmatrix} \label{eq:2.3.14} \end{equation}\]

は対称三重対角行列である。

すなわち、三重対角行列$T_m$について固有値計算を行うことにより、元の大規模な固有値問題の固有値を近似できる。