Three Dimensional Finite Element Method
eigenvaluelanczos
-19

(8) N(λ)がmodeより小さい場合はa1をλで置き換え、そうでない場合はa2をλで置き換えます。そして(5)へジャンプします。つまり、上の図からN(λ)がmodeより小さいとき、固有値はλとa2の範囲にあることが解ります。しかし、N(λ)がmodeより大きいか同じの場合、固有値はa1とλの範囲にあることになります。
(9) λをmodeの固有値とします。

以上、modeを1からnまでを繰り返すと全ての固有値が計算できます。例題で用いた[A]をLanczos法で [T]行列に変換し、そしてBisection法を適用すると、以下の固有値が得られます。上図の階段状になっている位置と一致していますよね。

値はこれまでと全く同じです。上で説明した計算手順を使うと、固有値の小さい順に計算されます。 プログラムLANCZOS_PRINCIPLE_BASIC.FORの中の SUBROUTINE BSECTIONを探し、Bisection法で固有値を計算する部分を見て下さい。そして、計算手順が上の説明と同じになっていることを確認して下さい。[T]マトリクスを生成せず、αiとβiのみで固有値を計算するBisection法のプログラム LANCZOS_PRINCIPLE_BASIC2.FORも作成しましたので見て下さい。
ところで、計算をスタートさせるために必要な{u}1ですが、 プログラムでは{u}1T={1/√(2) 0 1/√(2)}を使っています。しかし、他の値を つかってもかまいません。例えば、{u}1T={1 0 0}。|{u}1}|=1であれば 上と同じ固有値が計算されます。

Bisection法を数値的に扱う時にちょっとした心配ごとがあります。それは、Pi(λ)の計算値がオーバーフローかアンダーフローする可能性があるということです。これは、[T]を構成しているαiとβi値に依存します。Pi(λ)を完全に展開し記述するには、大変ですが、しかし、以下のような展開式になっていそうです。

つまり、Pi(λ)の値としては、αをi回掛け算することになります。例えば、αが10‐kだとすると、Pi(λ)=O(10‐ik)となりアンダーフローが心配されます。逆にαが10kだとすると、Pi(λ)=O(10ik)となりこんどはオーバーフローが発生することになります。 これを回避するには、αiとβiの値を事前にある数値(γ)で割っておいて、出来るだけα'iとβ'iを1に近づけておくことが必要になります。固有値λ'が計算されたら、真のλは、λ=λ'γで得られます。例えば、αiとβiの値を固有値の範囲の上限値bで割ります。すると、固有値の範囲は、[‐1, 1]になります。そして、固有値λ'が計算されたら、真の値はλ=λ'bということになります。しかし、この状況下でもアンダーフローの心配があります。この対処には、各々の固有値λ'を計算する際に最適なγ値を使うとよいでしょう。しかしこのような対処でもアンダー/オーバーフローの心配があります。あまり賢い対処方法ではありませんが、求めたい固有値の数mを少なくすると良いでしょう。例えばm=40とか。

次のセクションでは、Bisection法にNewton-Raphson法を活用し、計算速度を上げる方法を紹介します。

BACK NEXT
Menu LU Decompo Stiff 3D Solid 3D Fluid Eigen&Lanczos Sound Eigen Solid Eigen Solid Axisym