数値計算的には、繰り返し計算が必要になりますので、下式の右辺にある初期ベクトル{q0}iを何らかの方法で決めます。そして、以下の連立方程式を{q1}iについて計算します。
{q0}iに上で計算された{q1}iを代入し、再度、上式を計算します。
これを数回繰り返すと、{q0}iと{q1}iはほぼ同じになり、
{q}iが得られたことになります。
初期ベクトルですが、{q0}i={1,1,・・・,1}Tを使うと全ての固有ベクトル計算の初期値に使えるので便利です。
ここで新たに以下の行列を定義します。
そして、n回の計算で得られた{q1}iを{qn}iとします。 すると、{qn}iは{q0}iに対し以下の式で書くことができ、 丁度Krylov列の逆数になっています。
上の連立方程式を解くときですが、逆行列は計算しません。その代わりにLU分解法を使います。
そして、その連立方程式を解く際、det[[T]−λi[I]]は非常にゼロ近いので、
[[T]−λi[I]]の対角要素がLU分解の途中でゼロになる可能性があります。
そのときは、εの値を少し大きくして再度計算をトライしてみて下さい。実際の数値計算では、
ε=0でも問題なく計算できます。有限桁の計算には必ず誤差があるため、意図的にゼロにしない限り
対角要素がゼロになることはありません。
面倒なのは、重複固有値の場合です。私も対処法に悩みました。このあたりの代数的な対処法は、
文献等を調査してみて下さい。
ここでは本方法とBisection法とで併用する側の面からの対処法を考えてみましょう。
後で紹介する弾性解析の固有値解析で起こる可能性がある問題を取り上げます。
断面が2m(Y座標軸)×2m(Z座標軸)で長さが100m(X座標軸)のビームを取り上げます。
六面体1次要素を使いY座標軸方向に2分割、Z座標軸方向にも2分割、X座標軸方向は100分割で離散化したとします。
拘束条件はfree-freeです。
すると1次モードとしてY軸とZ軸方向に曲がる2つの固有値が存在します。これらの固有値は全く同じになります。
つまりこれら2つの固有値は重複固有値またはdet[A]の重根ということになります。これをどう回避するかが問題になります。
この場合の最も有効的な方法として、Y軸方向の分割数を2、Z軸方向の分割を3にすると固有値が変化します。
こうすることによりYとZ軸方向の固有値を互いに変えることができます。要素分割数は、固有値に敏感ですからね。
六面体2次要素で分割した場合も同じ対処方法でOKです。
話は固有ベクトル計算に戻りますが、ここで紹介した原点移動付逆ベキ乗法(Shifted Inverse Power Method)の収束速度ですが、 3回程度の繰り返しで所定の精度に到達します。しかし、アンダーフローが発生すると、繰り返し数は約15回に上昇します。
BACK | NEXT |
---|
Menu | LU Decompo | Stiff | 3D Solid | 3D Fluid | Eigen&Lanczos | Sound Eigen | Solid Eigen | Solid Axisym |