pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Maximaで行列の固有値・固有ベクトルを求め、対角化する

数学の具体的な計算にMaximaを使って、数学もMaximaも同時に学んでしまいましょう。今回はMaximaを使って行列の計算をしてみたいと思います。行列の固有値固有ベクトルを求め、対角化する方法をみてみます。

Maximaのごく基本的な使い方については以下の記事を参照してください:

pianofisica.hatenablog.com

また、Maximaでの行列の扱いの基本については

pianofisica.hatenablog.com

にまとめていますので、参考にしてみてください。



固有値固有ベクトル

行列  A固有値固有ベクトルとは、方程式

 Au=\lambda u

をみたす数  \lambda とベクトル  u のことをいいます。

非自明な固有ベクトル  u が存在するためには

 (A-\lambda\,I)u=0

としたとき、行列

 A-\lambda\,I

の逆元が存在してはいけません。なぜなら、もし逆元が存在すれば、それを両辺に作用させて  u=0 が導かれてしまい、  u が非自明であることに反してしまうからです。

以上から、固有値に対して以下の固有方程式

 \det(A-\lambda\,I)=0

が要請されます。この固有方程式を  \lambda について解くことにより、行列  A固有値が求まります。

一般に、 n 次正方行列に対しては固有方程式は  \lambda n 次方程式になるので、複素数の範囲では(縮重度を含めて) n 個の固有値があります。

行列の対角化

以下では固有値は互いに全て異なるものとします。

固有値  \lambda_i \ (i=1,\dots,n) に対する固有ベクトル  u_i から作られる  n 次正方行列

 \displaystyle{P=\left(\begin{array}{cccc} u_1 & u_2 & \dots & u_n \end{array}\right)}

を導入すると、この行列  P は正則であることがわかります。このことは、互いに異なる固有値固有ベクトル  u_i が線型独立であることから結論できます。

一番簡単な場合として、2次正方行列の場合を考えてみましょう。このとき一般に2つの固有値固有ベクトルがあります:

 Au_1=\lambda_1 u_1 \qquad Au_2=\lambda_2 u_2

ここで2つの固有値が互いに異なるとき、  u_1 u_2 は線型独立です。

これを背理法で示してみましょう。すなわち、2つの固有値が互いに異なるときに  u_1 u_2 が線型独立でない、と仮定してみます。すると、ある数  c が存在して

 u_2=c u_1

と書けることになります。仮定から、  u_1 u_2 は行列  A固有ベクトルでしたから

 \lambda_2u_2 =A u_2=A(cu_1)=cAu_1=c\lambda_1u_1=\lambda_1(cu_1)=\lambda_1u_2

です。以上から

 (\lambda_2-\lambda_1)u_2=0

となります。いま、固有ベクトル  u_2 は非自明ですから、この式が成り立つためには

 \lambda_1=\lambda_2

が導かれます。ところがこれは2つの固有値が互いに異なるという前提と矛盾します。よって  u_1 u_2 は線型独立であることがわかります。

 

同様の考察によって一般の場合についても、互いに異なる固有値固有ベクトル  u_i が線型独立であることがわかります。したがって線型独立なベクトルから構成された行列  P は、その行列式がゼロではあり得ず、その逆行列  P^{-1} が存在すること、つまり正則であることがわかります。

さて、 u_i が行列  A固有ベクトルであることから

 \displaystyle{AP=\left(\begin{array}{cccc} \lambda_1u_1 & \lambda_2u_2 & \dots &\lambda_n u_n \end{array}\right)}

より、これに左から行列  P^{-1} を作用させて

 \displaystyle{\Lambda=P^{-1}AP=\left(\begin{array}{cccc} \lambda_1 & 0 & \dots & 0 \\  0 & \lambda_2 & \dots & 0 \\ \vdots & & & \vdots \\  0 & 0 & \dots & \lambda_n \end{array}\right)}

という対角行列を作ることができます。このような対角行列を作ることのご利益ですが、行列  A のベキ乗が

 A^n=(P\Lambda P^{-1})^n=P\Lambda^nP^{-1}

によって簡単に計算できるということがあげられます。



具体例:2次正方行列の場合

固有値固有ベクトルを具体的にみるために行列

 \displaystyle{A=\left(\begin{array}{cc} 1 & 2 \\ 3 & 2 \end{array}\right)}

を考えてみます。行列をMaximaに入力するには

A:matrix([1,2],[3,2]);

とします。

固有方程式  \det(A-\lambda,I)=0行列式

 \displaystyle{\det\left(\begin{array}{c} 1-\lambda & 2 \\ 3 & 2-\lambda \end{array}\right)}

D:determinant(A-lambda*ident(2));

をゼロとして lambda について方程式を解きます:

L:solve(D=0, lambda);

により、固有値

 \displaystyle{\lambda_1=4 \qquad \lambda_2=-1}

と求まります。これらの固有値に対する固有ベクトルを求めていきます。

まず固有値  \lambda_1=4 に対する固有ベクトル  u_1={}^t(x_1,y_1)

 \displaystyle{Au_1=-4u_1}

u1:matrix([x1],[y1])$
B1:A.u1$
solve([B1[1][1]=rhs(L[1])*u1[1][1],B1[2][1]=rhs(L[1])*u1[2][1]],[x1,y1]);

より、ベクトルの成分ついて

 \displaystyle{x_1=\frac{2}{3}y_1}

という条件が課されます。よって、例えば

 \displaystyle{x_1=2\qquad y_1=3 }

ととることができます。同様にして固有値  \lambda_2=-1 に対する固有ベクトル  u_2={}^t(x_2,y_2)

u2:matrix([x2],[y2])$
B2:A.u2$
X:solve([B2[1][1]=rhs(L[2])*u2[1][1],B2[2][1]=rhs(L[2])*u2[2][1]],[x2,y2]);

より

 \displaystyle{x_2=1 \qquad y_2=-1 }

ととれます。

 

さて、行列  A を対角化するために行列  P

 \displaystyle{P=\left(\begin{array}{cc} u_1 & u_2 \end{array}\right)=\left(\begin{array}{cc} 2 & 1 \\ 3 & -1 \end{array}\right)}

によって導入します:

P:matrix([2,1],[3,-1]);

すると

invert(P).A.P;

より

 \displaystyle{\Lambda=P^{-1}AP=\left(\begin{array}{cc} 4 & 0 \\ 0 & -1 \end{array}\right)}

と計算されます。ただし

invert(P);

より

 \displaystyle{P^{-1}=\left(\begin{array}{cc} \frac15 & \frac15 \\ \frac35 & -\frac25 \end{array}\right)}

です。

ここで  A^n を考えると

 \displaystyle{A^n=P\Lambda^n P^{-1}=\frac15\left(\begin{array}{cc} 2 & 1 \\ 3 & -1 \end{array}\right)\left(\begin{array}{cc} 4^n & 0 \\ 0 & (-1)^n \end{array}\right)\left(\begin{array}{cc} 1 & 1 \\ 3 & -2 \end{array}\right)=\frac15\left(\begin{array}{cc} 2\,4^{n}+3(-1)^n & 2\,4^{n}-2(-1)^n \\ 3\,4^{n}-3(-1)^n & 3\,4^n+2(-1)^n \end{array}\right)}

と計算できます。実際

Apow[n]:=1/5*matrix([2*4^n+3*(-1)^n, 2*4^n-2*(-1)^n],[3*4^n-3*(-1)^n, 3*4^n+2*(-1)^n])$

から

Apow[0];
Apow[1];
Apow[10];
A^^10;

などによって確認することができます。

行列の対角化(Maximaの組み込み)

固有値固有ベクトルは原理的には以上のようにして求められるのですが、これらを求める処理がMaximaにはすでに組み込まれています:まず固有値ですが

eigenvalues(A);

により、その値と縮重度を知ることができます。さらに

eigenvectors(A);

では、固有値とその縮重度に加え、それらに対応する固有ベクトルを出力します。

よってこれまでの計算は、例えば、上で考えた  A とはまた違った行列

 \displaystyle{B=\left(\begin{array}{cc} 1 & 2 \\ 2 & 1 \end{array}\right)}

に対して

B:matrix([1,2],[2,1])$
E:eigenvectors(B)$
v1:matrix([E[2][1][1][1]],[E[2][1][1][2]])$
v2:matrix([E[2][2][1][1]],[E[2][2][1][2]])$
Q:matrix([v1[1][1],v2[1][1]],[v1[2][1],v2[2][1]]);
Lamb:ratsimp(invert(Q).B.Q);

などとして計算できます。

入力する行列を変えることで、それに応じた対角行列が出力されます。

エラーが出る場合は、固有値に縮重度があって線型独立な固有ベクトルがとれない場合になります。

縮重度が2で線型独立な固有ベクトルがない場合は例えば

C1:matrix([1,2],[-2,5])$
eigenvectors(C1);

などです。詳しくは以下の記事で解説しています:

pianofisica.hatenablog.com

 


3次正方行列への拡張

これまで  2\times 2 行列についてみてきましたが、  3\times 3 行列に拡張するには

B:matrix([1,2,1],[2,1,1],[0,2,1])$
E:eigenvectors(B)$
v1:matrix([E[2][1][1][1]],[E[2][1][1][2]],[E[2][1][1][3]])$
v2:matrix([E[2][2][1][1]],[E[2][2][1][2]],[E[2][2][1][3]])$
v3:matrix([E[2][3][1][1]],[E[2][3][1][2]],[E[2][3][1][3]])$
Q:expand(ratsimp(matrix(
[v1[1][1],v2[1][1],v3[1][1]],
[v1[2][1],v2[2][1],v3[2][1]],
[v1[3][1],v2[3][1],v3[3][1]])));
Lamb:expand(ratsimp(invert(Q).B.Q));

などとして同様に計算できます。


本記事で学んだ内容の応用として、次の記事もあわせて読んでみてください:

pianofisica.hatenablog.com



キーワードMaxima、行列、固有値、対角化