数値解析 第14回 (1) 手計算と同じ方法
手計算のアルゴリズム
$n$ 次正方行列 $A$ の固有値・固有ベクトル $(\lambda_j,\vvv_j)$ ( $j = 1,2,\cdots,n$ ) を手計算で
求めるアルゴリズムは次の通りでした:
Algorithm 1
- 固有多項式 $\varphi_A(x)=\det(xE-A)$ を求める
- 固有方程式 $\varphi_A(x)=0$ の $n$ 個の解 $\lambda_1$, $\lambda_2$, $\cdots$, $\lambda_n$ を求める
- 各 $\lambda_j$ に対し $(\lambda_j E-A)\vvv=\ooo$ の非零解 $\vvv_j$ を求める
これをコンピュータで実現しましょう。
1° 固有多項式 $\varphi_A(x)$ の計算
多項式の計算なので数式処理ソフトを使おう、と思うかもしれませんが、
汎用性とのトレードオフで、行列サイズが大きくなると結構な計算時間が掛かるのが普通です。
次のアルゴリズムを使うと数式処理ソフトが無くても大丈夫です:
Algorithm 2 ( Frame 法 ) ( Faddeev-Leverrier 法 とも言う )
入力:$n$ 次正方行列 $A$
出力:$A$ の固有多項式 $\varphi_A(x)=x^n+c_1\,x^{n-1}+\cdots+c_n$ の係数
$X=E$
for $k=1$ to $n$
$X=AX$
$c_k=-$tr$(X)\,/\,k$
$X=X+c_kE$
ただし tr は行列の対角成分和を表す。
Ex.3 2次行列 $\dps{A=\mat{cc}{a & b \\ c & d}}$ の場合
実行例を折りたたむ
- $X$ の初期値 $\dps{=\mat{cc}{1 & 0 \\ 0 & 1}}$
- $k=1$ :
-
新 $\dps{X=AX=\mat{cc}{a & b \\ c & d}}$,
tr$(X)=a+d$,
$c_1=-(a+d)\,/\,1=-(a+d)$,
新 $\dps{X=X+c_1E=\mat{cc}{-d & b \\ c & -a}}$
- $k=2$ :
-
新 $\dps{X=AX=\mat{cc}{a & b \\ c & d}\mat{cc}{-d & b \\ c & -a}=\mat{cc}{-ad+bc & 0 \\ 0 & -ad+bc}}$,
tr$(X)=2(-ad+bc)$,
$c_2=-2(-ad+bc)\, /\, 2=ad-bc$
これで
$\varphi_A(x)=x^2-(a+d)x+(ad-bc)$
が求まりました。
Ex.4 5次行列
$\dps{
A=\mat{rrrrr}{\
-4 & 6 & -9 & -6 & -8\ \\[-0.2em]
4 & -9 & 5 & 4 & 4\ \\[-0.2em]
4 & -2 & 8 & 3 & -1\ \\[-0.2em]
1 & 0 & 9 & -1 & 2\ \\[-0.2em]
-2 & -7 & 1 & 1 & 6\
}
}$ の場合
実行例を折りたたむ
k = 1
X: -4 6 -9 -6 -8
4 -9 5 4 4
4 -2 8 3 -1
1 0 9 -1 2
-2 -7 1 1 6
c[1] = 0
k = 2
X: 14 -4 -68 19 5
-36 67 -1 -45 -41
13 33 44 -12 -48
27 -26 56 24 -7
-27 7 6 -8 25
c[2] = -87
k = 3
X: 13 -301 269 204 704
97 253 -230 137 -127
-8 203 -452 -111 -241
-37 333 -499 -42 -544
276 197 192 154 -150
c[3] = 126
k = 4
X: -1888 -1923 1936 -731 2047
599 -1480 288 -20 482
-89 464 -2761 -248 -466
656 1587 -1782 -571 -969
654 -333 1399 -470 -1448
c[4] = 2037
k = 5
X: -5369 0 0 0 0
0 -5369 0 0 0
0 0 -5369 0 0
0 0 0 -5369 0
0 0 0 0 -5369
c[5] = 5369
(
Alg.2 の証明は例えば Graduate Texts in Mathematics 138, pp.53-54 にあります。)
2° 固有方程式 $\varphi_A(x)=0$ の解の計算
これはニュートン法でできます。
第3回 Rem.12 で述べたように
複素数解もニュートン法で得られますが、
C言語なら complex.h を include するか、自分で複素数ライブラリを作るかしなければなりません。
( python なら complex 型が標準で使えます。)
3° $(\lambda_j E-A)\vvv=\ooo$ の非零解の計算
$(\lambda_j E-A)$ が逆行列を持ちませんので、
この講義で紹介したガウスの消去法などはそのままでは使えません。
ひとつの手としては
- $i=$ ランダムな番号 ( $1 \leqq i \leqq n$ )
- $\vvv$ の第 $i$ 成分 $=1$
と仮定して、$\vvv$ の他の $n-1$ 個の成分に関する連立一次方程式を解く、
ということが考えられます。( 固有値 $\lambda_j$ が重解である場合にはさらに工夫が必要です。)