応用数理実験
固有値問題 (1) 固有値問題に対する冪乗法

桂田 祐史

2005年6月28日

今回は固有値問題を取り扱うための冪乗法とその周辺の種々の手法を実験し てもらう。

実験に用いる行列は小さなもので良い。自分で選ぶのが面倒なら、次の行列 を使ってよい¹。

$$A= \left( \begin{array}{ccc} 1& 1& 1 \\ 1& 2& 2 \\ 1& 2& 3 \end{array}\right).$$

• 選んだ行列のサイズが小さい場合は、固有多項式を計算して固有値に関する 情報を調べよ²。 サイズが小さくない場合は、 Mathematica か (後で例をあげる)、 MATLAB (Octave) の関数 eig() を用いて固有値を調べよ。
• 冪乗法により絶対値最大の固有値とそれに属する固有ベクトルを求めよ。収 束の速さを調べよ。 固有値を反復ベクトルの成分の比として計算した場合と、 Rayleigh 商で計算 した場合とで、収束の速さに違いが出るか？
• 逆反復法により (連立 $1$ 次方程式を解くのには LU 分解を使え)、絶対値 最小の固有値とそれに属する固有ベクトルを求めよ。
• シフト法の実験をせよ。すなわち絶対値最大、最小の固有値以外の固有値を 見い出せ。(例えば、上に掲げた行列 $A$ に対しては $0.6$ に近い固有値を 探してみよ。)
• (余裕がある人向け) 絶対値最大の固有値に属する固有ベクトルを求めてから、 その成分を含まない初期ベクトルを用いて冪乗法の反復を実行することにより、 絶対値が二番目に大きな固有値の固有ベクトルを求めるアルゴリズムがある。 適当な文献で調べるか、自分で考えて実験してみよ。
• (余裕がある人向け) 絶対値最大の固有値が二つ以上ある場合、 冪乗法で何が起こるか観察せよ。

Mathmatica の使用例

oyabun%	math
Mathematica 4.0 for Solaris
Copyright 1988-1999 Wolfram Research, Inc.
 -- Motif graphics initialized -- 

In[1]:= A={{1,1,1},{1,2,2},{1,2,3}}

Out[1]= {{1, 1, 1}, {1, 2, 2}, {1, 2, 3}}

In[2]:= Eigenvalues[A]

Out[2]= (省略 --- 興味があればやってみて下さい)

In[3]:= N[%,40]
                                                          -46
Out[3]= {5.048917339522305313522214407023369723596 + 0. 10    I, 

                                                       -42
>    0.3079785283699041303721851029979308598027 + 0. 10    I, 

                                                       -42
>    0.6431041321077905561056004899786994166009 + 0. 10    I}

In[4]:= na=N[A,50]

Out[4]= {{1.0000000000000000000000000000000000000000000000000, 
>     1.0000000000000000000000000000000000000000000000000, 
>     1.0000000000000000000000000000000000000000000000000}, 
>    {1.0000000000000000000000000000000000000000000000000, 
>     2.0000000000000000000000000000000000000000000000000, 
>     2.0000000000000000000000000000000000000000000000000}, 
>    {1.0000000000000000000000000000000000000000000000000, 
>     2.0000000000000000000000000000000000000000000000000, 
>     3.0000000000000000000000000000000000000000000000000}}

In[5]:= Eigenvalues[na]

Out[5]= {5.0489173395223053135222144070233697235963877860565, 
>    0.64310413210779055610560048997869941660087281326538, 
>    0.30797852836990413037218510299793085980273940067811}

In[6]:= Eigenvectors[na]

Out[6]= {{-0.32798527760568176779603202500454990640870532112133, 
>     -0.59100904850610352545794571799937980844713881955014, 
>     -0.73697622909957824233808630700517009796156650157119}, 
>    {0.73697622909957824233808630700517009796156650157119, 
>     0.32798527760568176779603202500454990640870532112133, 
>     -0.59100904850610352545794571799937980844713881955014}, 
>    {-0.59100904850610352545794571799937980844713881955014, 
>     0.73697622909957824233808630700517009796156650157119, 
>     -0.32798527760568176779603202500454990640870532112133}}

In[7]:= NSolve[Det[x IdentityMatrix[3]-A]==0, x, WorkingPrecision->30]}

Out[7]= {{x -> 0.30797852836990413037218510300},
         {x -> 0.64310413210779055610560048998},
         {x -> 5.0489173395223053135222144070}}