固有値

正方行列 A のシンボリックな固有値、または正方行列 A のシンボリックな固有値と固有ベクトルは、それぞれコマンド E = eig(A) および [V,E] = eig(A) で計算されます。

可変精度の行列の場合は、それぞれコマンド E = eig(vpa(A)) および [V,E] = eig(vpa(A)) で計算されます。

A の固有値は、A の特性多項式 det(A-x*I) の零点で、charpoly(A) で計算されます。

最初の例として、前の節の行列 H を使います。

H = sym([8/9 1/2 1/3; 1/2 1/3 1/4; 1/3 1/4 1/5])

H =
[ 8/9, 1/2, 1/3]
[ 1/2, 1/3, 1/4]
[ 1/3, 1/4, 1/5]

これは特異行列なので、その固有値の 1 つは 0 でなければなりません。次のステートメントを見てみましょう。

[T,E] = eig(H)

行列 T と E が作成されます。T の列は H の固有ベクトルで、E の対角要素は H の固有値です。

T =
[ 3/10, 218/285 - (4*12589^(1/2))/285, (4*12589^(1/2))/285 + 218/285]
[ -6/5,     292/285 - 12589^(1/2)/285,     12589^(1/2)/285 + 292/285]
[    1,                             1,                             1]
 
E =
[ 0,                       0,                       0]
[ 0, 32/45 - 12589^(1/2)/180,                       0]
[ 0,                       0, 12589^(1/2)/180 + 32/45]

T や E を 10 進数に変換すると、固有ベクトル T と固有値 E の行列の構造が理解しやすくなる可能性があります。そこで、次のコマンドを使います。コマンド

Td = double(T)
Ed = double(E)

は以下を返します。

Td =
    0.3000   -0.8098    2.3397
   -1.2000    0.6309    1.4182
    1.0000    1.0000    1.0000

Ed =
         0         0         0
         0    0.0878         0
         0         0    1.3344

最初の固有値はゼロです。対応する固有ベクトルは Td の 1 列目で、前の節で求めたヌル空間の基底と同じです。他の 2 つの固有値は、factor(charpoly(H, x)) の 2 次因子である $x^{2} - \frac{64}{45} x + \frac{253}{2160}$ に 2 次の公式を適用した結果です。

syms x
g = factor(charpoly(H, x))/x
solve(g(3))

g =
[ 1/(2160*x), 1, (2160*x^2 - 3072*x + 253)/x]
ans =
 32/45 - 12589^(1/2)/180
 12589^(1/2)/180 + 32/45

特性多項式が 4 次以下の有理数多項式の積として表されるときのみ、固有値に対する閉じた形のシンボリック式が得られます。古典的な数値解析のテスト行列である Rosser 行列を使って、この条件を示します。次のステートメントを見てみましょう。

R = sym(rosser)

は、次を生成します。

R =
[  611,  196, -192,  407,   -8,  -52,  -49,   29]
[  196,  899,  113, -192,  -71,  -43,   -8,  -44]
[ -192,  113,  899,  196,   61,   49,    8,   52]
[  407, -192,  196,  611,    8,   44,   59,  -23]
[   -8,  -71,   61,    8,  411, -599,  208,  208]
[  -52,  -43,   49,   44, -599,  411,  208,  208]
[  -49,   -8,    8,   59,  208,  208,   99, -911]
[   29,  -44,   52,  -23,  208,  208, -911,   99]

コマンド

p = charpoly(R, x);
factor(p)

結果は以下のようになります。

ans =
 
[ x, x - 1020, x^2 - 1040500, x^2 - 1020*x + 100, x - 1000, x - 1000]

この 8 次の特性多項式は、きれいに 2 つの線形な項と 3 つの 2 次の項の積で表されます。この結果から、固有値のうちの 4 つは、0、1020 および 1000 における重根であることがすぐにわかります。他の 4 つの根は、残りの 2 次式から得られます。次を使用します。

eig(R)

すべての固有値が得られます。

ans =
                  0
               1000
               1000
               1020
 510 - 100*26^(1/2)
 100*26^(1/2) + 510
    -10*10405^(1/2)
     10*10405^(1/2)

Rosser 行列は典型的な例ではありません。8 行 8 列の行列が、このような簡単な形に因数分解できる特性多項式をもつことは稀です。次のコマンドを使って行列 R の 2 隅の要素を 29 から 30 に変更し、

S = R;
S(1,8) = 30;
S(8,1) = 30;

さらに次のコマンドを実行すると、

p = charpoly(S, x)

次が得られます。

p =
x^8 - 4040*x^7 + 5079941*x^6 + 82706090*x^5...
 - 5327831918568*x^4 + 4287832912719760*x^3...
 - 1082699388411166000*x^2 + 51264008540948000*x...
 + 40250968213600000

また、factor(p) が p 自体であることがわかります。つまり、この特性多項式は有理数上で因数分解できません。

この修正した Rosser 行列に対して、

F = eig(S)

は以下を返します。

F = 
 -1020.053214255892
  -0.17053529728769
 0.2180398054830161
  999.9469178604428
  1000.120698293384
  1019.524355263202
  1019.993550129163
  1020.420188201505

これらの値は、オリジナルの Rosser 行列の固有値に非常に近い値です。

シンボリック行列の固有値を求めようとすることもできますが、閉じた形の解が得られることはほとんどありません。ギブンス回転行列は、次の初等行列の行列指数として生成されます。

$A = [\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix}] .$

Symbolic Math Toolbox™ コマンド

syms t
A = sym([0 1; -1 0]);
G = expm(t*A)

は以下を返します。

G =
[           exp(-t*1i)/2 + exp(t*1i)/2,
    (exp(-t*1i)*1i)/2 - (exp(t*1i)*1i)/2]
[ - (exp(-t*1i)*1i)/2 + (exp(t*1i)*1i)/2,
            exp(-t*1i)/2 + exp(t*1i)/2]

この式は simplify を使用して単純化することができます。

G = simplify(G)

G = 
[  cos(t), sin(t)]
[ -sin(t), cos(t)]

次に、コマンド

g = eig(G)

は、以下の結果を出力します。

g =
 cos(t) - sin(t)*1i
 cos(t) + sin(t)*1i

g は指数を用いて書き直すことができます。

g = rewrite(g, 'exp')

g = 
 exp(-t*1i)
  exp(t*1i)