8 ケーリー・ハミルトンの公式の応用

本節では、行列乗を計算するもう一つの方法を紹介する。 これは、wikipedia のケーリー・ハミルトンの項 ([2]) に 書かれているもので、実は私もこれまで知らなかった。

$n$ 次正方行列に対しては、ケーリー・ハミルトンの公式が成り立つことが 知られている。すなわち、 $A\in M_n(\mbox{\boldmath$K$})$ の固有多項式を $f_A(x)=\vert xE-A\vert$ ($n$ 次式) とすると、

  $$f_A(A)=O$$ (30)

となる。

これにより、$A^n$ が $E,A,A^2,\ldots,A^{n-1}$ の一次結合、 すなわち $A$ の $(n-1)$ 次以下の多項式として表され、 さらに任意の $A^m$ ($m\geq 1$) も同様に $(n-1)$ 次以下の $A$ の多項式で 表されることになる。よって、そのような項の無限和である $e^{At}$ も $A$ の $(n-1)$ 次以下の多項式で表される、 ということが想像されるが、それを考えてみる。

まず、「無限次の多項式」のような解析関数 (ベキ級数) に対する剰余計算 に関する定理を紹介する。

定理 8.1

関数 $g(x)$ が

$$g(x) = \sum_{m=0}^\infty \alpha_m x^m$$

で定義され $\vert x\vert<r$ で収束する解析関数であるとし、 $p(x)$ を $n$ 次多項式

$$p(x)=x^n+p_{n-1}x^{n-1}+\cdots+p_0$$

とする。このとき、解析関数 $q(x)$ (商) と $(n-1)$ 次以下の 多項式 $r(x)$ (余り) が存在して、 $g(x) = q(x)p(x)+r(x)$ が成り立つ。

なお、$\alpha_m$, $p_j$ は一般に複素数とする。

証明

$p(x)$ を ($C$の範囲で) 1 次式の積に因数分解する。

$$p(x)=(x-c_1)(x-c_2)\cdots (x-c_n)$$

このとき、

$$q_1(x) = \left\{\begin{array}{ll} \displaystyle \frac{g(x)-g(... ...x\neq c_1$\ のとき})\\ [2ex] g'(c_1) & (\mbox{$x=c_1$\ のとき}) \end{array}\right.$$

とすると、$x=c_1$ は $q_1(x)$ の除去可能特異点で、 $q_1(x)$ も $\vert x\vert<r$ で正則となり、よって $\vert x\vert<r$ で解析的となる。 そして、

  $$g(x) = q_1(x)(x-c_1)+g(c_1)$$ (31)

と書ける。 同様に、

$$q_2(x) = \left\{\begin{array}{ll} \displaystyle \frac{q_1(x)-... ...neq c_2$\ のとき})\\ [2ex] q_1'(c_2) & (\mbox{$x=c_2$\ のとき}) \end{array}\right.$$

とすれば、$q_2(x)$ も $\vert x\vert<r$ で解析的で、

$$q_1(x)=q_2(x)(x-c_2)+q_1(c_2)$$

と書ける。これを繰り返せば、最後に

$$q_{n-1}(x)=q_n(x)(x-c_n)+q_{n-1}(c_n)$$

となり、これらを (31) に代入すれば、

$$\begin{eqnarray*}r(x) &=& g(c_1)+q_1(c_2)(x-c_1)+q_2(c_3)(x-c_1)(x-c_2)+\cdots \\ &&\mbox{}+q_{n-1}(c_n)(x-c_1)\cdots(x-c_{n-1}) \end{eqnarray*}$$

と $q(x)=q_n(x)$ に対して $g(x) = q(x)p(x)+r(x)$ が成り立つ。

この定理 8.1 により、

  $$e^{\alpha x} = q_{A,\alpha}(x)f_A(x)+r_{A,\alpha}(x)$$ (32)

となるような解析関数 $q_{A,\alpha}(x)$ と $(n-1)$ 次以下の 多項式 $r_{A,\alpha}(x)$ が存在することになる。

今、$A$ の固有値を $\lambda_1,\ldots,\lambda_n$ とすると、 $f_A(\lambda_j)=0$ なので、(32) に代入すると

  $$e^{\alpha\lambda_j}=r_{A,\alpha}(\lambda_j) \hspace{1zw}(1\leq j\leq n)$$ (33)

が成り立つことになる。 $r_{A,\alpha}(x)$ は高々 $(n-1)$ 次式であるからその係数は $n$ 個以下で、 よってもしすべての固有値 $\lambda_j$ が異なれば、 (33) からその係数を決定でき ( $e^{\alpha\lambda_j}$ で表され)、 $r_{A,\alpha}(x)$ が求まる。 すなわち、未定係数 $\beta_j$ を用いて

$$r_{A,\alpha}(x)=\beta_0+\beta_1 x+\cdots+\beta_{n-1}x^{n-1}$$

とすれば、(33) より $\beta_0,\ldots,\beta_{n-1}$ に 対する連立方程式

$$\left[\begin{array}{cccc}1&\lambda_1&\cdots&\lambda_1^{n-1}\\ \... ...{array}{c}e^{\alpha\lambda_1}\\ \vdots\\ e^{\alpha\lambda_n}\end{array}\right]$$

が得られ、これを $\beta_0,\ldots,\beta_{n-1}$ について解けばよい。 ちなみに、この方程式の係数行列はファンデルモンド行列と 呼ばれるものになっていて、 $\lambda_1,\ldots,\lambda_n$ が すべて異なれば、確かに逆行列が存在し、 $\beta_0,\ldots,\beta_{n-1}$ が 求まる。

次に固有値に重解がある場合を考える。 固有値 $x=\lambda_j$ が $f_A(x)=0$ の $k$ 重解である場合、 (33) の方程式のうち $k$ 本が 1 本の同じ方程式に なってしまうが、 この場合 $f_A(x)$ は $(x-\lambda_j)^k$ で割り切れるので、 $f_A$ の $(k-1)$ 階までの導関数に対して

$$f_A^{(i)}(\lambda_j)=0\hspace{1zw}(0\leq i\leq k-1)$$

が成り立つ。 (32) を $i$ 階微分すると ( $1\leq i\leq k-1$)、

$$\alpha^ie^{\alpha x} =\sum_{\ell=0}^i\left(\begin{array}{c} \!\!... ...{array}\right)q_{A,\alpha}^{(i-\ell)}(x)f_A^{(\ell)}(x) +r_{A,\alpha}^{(i)}(x)$$

となるので、$x=\lambda_j$ とすると $f_A^{(\ell)}(\lambda_j)=0$ より

  $$\alpha^ie^{\alpha\lambda_j}=r_{A,\alpha}^{(i)}(\lambda_j) \hspace{1zw}(1\leq i\leq k-1)$$ (34)

となり、これにより足りない $(k-1)$ 本の方程式が補えることになる。

以上により、固有値に重解がある場合も含めて $r_{A,\alpha}(x)$ が求められ ることになり、 そして (32) に $x=A$ を代入すると、 (30) により

  $$e^{\alpha A} = r_{A,\alpha}(A) = \beta_0E + \beta_1 A+\cdots+\beta_{n-1}A^{n-1}$$ (35)

のように、行列乗が $A$ の $(n-1)$ 次以下の多項式で表されることになる。

前節で計算した $e^{At}$ を、この方法で計算してみる。 まずは、

$$A=\left[\begin{array}{cc}0&1\\ -4&0\end{array}\right]$$

の場合。この場合、$f_A(x)=x^2+4$, $\lambda=\pm 2i$ だったので、 $r_{A,\alpha}(x)=\beta_0+\beta_1 x$ とすると、 (33) は

$$e^{-2i\alpha}=r_{A,\alpha}(-2i) = \beta_0-2i\beta_1, \hspace{1zw} e^{2i\alpha}=r_{A,\alpha}(2i) = \beta_0+2i\beta_1$$

となるので、

$$\beta_0 = \frac{e^{-2i\alpha}+e^{2i\alpha}}{2}=\cos 2\alpha, \hspace{1zw} \beta_1 = \frac{e^{2i\alpha}-e^{-2i\alpha}}{4i}=\frac{\sin 2\alpha}{2}$$

となり、よって

$$r_{A,\alpha}(x) = \cos 2\alpha + \frac{x}{2}\,\sin 2\alpha$$

と求まるから、(35) より

$$e^{\alpha A} = r_{A,\alpha}(A) = E\cos 2\alpha + \frac{A}{2}\sin... ...{cc}\cos2\alpha &\sin2\alpha/2\\ -2\sin 2\alpha &\cos2\alpha\end{array}\right]$$

となって、実質的に (29) と同じものが 求まっていることがわかる。しかもこちらの方が計算は楽そうである。

次に、

$$A = \left[\begin{array}{cccc}0&0&1&0\\ 0&0&0&1\\ 1&2&0&2\\ 1&-1&3&0\end{array}\right]$$

の場合を計算してみる。この場合は、 $f_A(x)=x^4-6x^2-8x-3$ であり、 固有値は $-1$ (3 重根) と $3$ であったので、

$$r_{A,\alpha}(x)=\beta_0+\beta_1 x+\beta_2 x^2 + \beta_3 x^3$$

とすると、 (33), (34) により

$$\left\{\begin{array}{ll} e^{3\alpha} & =r(3) = \beta_0+3\beta_1+... ...a_3,\\ \alpha^2 e^{-\alpha} & = r''(-1) = 2\beta_2-6\beta_3\end{array}\right.$$

となるので、連立方程式

$$\left[\begin{array}{cccc}1&3&9&27\\ 1&-1&1&-1\\ 0&1&-2&3\\ 0&0&2&... ...ha}\\ e^{-\alpha}\\ \alpha e^{-\alpha}\\ \alpha^2e^{-\alpha}\end{array}\right]$$

を解けばよい。 それはクラメルの公式でも、通常の消去法でも計算できるが、 ここでは $\alpha$ を含まない数だけの計算にするために、 この係数行列 $B$ の逆行列を消去法で計算する。 これも一応分数を出さないように計算してみる。

$$\begin{eqnarray*}\lefteqn{[B\ \vert\ E] \ =\ \left[\begin{array}{cccc\vert cc... ...0\\ 0&0&64&0 &3&-3&-12&8\\ 0&0&0&64 &1&-1&-4&-8\end{array}\right]\end{eqnarray*}$$

となるので、よって、

$$\left[\begin{array}{c}\beta_0\\ \beta_1\\ \beta_2\\ \beta_3\end{array}\right] = \frac{1}{64} \left[\begin{array}{cccc}1&63&60&24\\ 3&-3&52&40\\ ... ...}\\ e^{-\alpha}\\ \alpha e^{-\alpha}\\ \alpha^2e^{-\alpha}\end{array}\right]$$ (36)

$$= \frac{1}{64} \left[\begin{array}{c}e^{3\alpha}+(63+60\alpha+24\al... ...2)e^{-\alpha}\\ e^{3\alpha}-(1+4\alpha+8\alpha^2)e^{-\alpha}\end{array}\right]$$

と各 $\beta_j$ が求まる。この $\beta_j$ に対して

$$e^{\alpha A} = \beta_0E + \beta_1 A+\beta_2 A^2 + \beta_3 A^3$$

と表されることになるので、あとは $A^2$, $A^3$ を計算すればよい。

$$\begin{eqnarray*}A^2 &=& \left[\begin{array}{cccc}0&0&1&0\\ 0&0&0&1\\ 1&2&0&2\... ...cccc}2&-2&7&2\\ 3&6&1&5\\ 9&12&8&12\\ 6&-3&18&8\end{array}\right]\end{eqnarray*}$$

となる。さらに最終形として、$\alpha$ によらない定数行列 $P_j$ により

$$e^{\alpha A} = \beta_0E + \beta_1 A+\beta_2 A^2 + \beta_3 A^3 = P_1e^{3\alpha}+P_2e^{-\alpha}+P_3\alpha e^{-\alpha}+P_4\alpha^2 e^{-\alpha}$$

の形にするなら、(36) より

  $$\left\{\begin{array}{ll} P_1 &= \displaystyle \frac{1}{64}(E+3A... ...ac{1}{64}(24E+40A+8A^2-8A^3) \ =\frac{1}{8}(3E+5A+A^2-A^3) \end{array}\right.$$ (37)

を計算すればよい。

$$\begin{eqnarray*}\lefteqn{64P_1 = E+3A+3A^2+A^3} \\ &=&\hspace{-.8em} \left... ...ccc}2&4&-2&0\\ -2&-4&2&0\\ -2&-4&2&0\\ 2&4&-2&0\end{array}\right]\end{eqnarray*}$$

となる。

こちらの計算では、固有ベクトル、広義固有ベクトルの計算はないが、 $\beta_j$ を求める連立方程式の計算と $A^{n-1}$ の計算が必要で、 さらに $P_j$ による表現を求める場合は (37) の行列の 定数倍の和の計算が必要になる。 理論展開にはジョルダン標準形の方が良いかもしれないが、 具体的な $e^{At}$ の計算をするには、 こちらのケーリー・ハミルトンの公式を利用する方が最終目的にも少し近く、 計算もやや易しい気がする。

なお、$A^k$ を経由せずに直接 $P_j$ が計算できればもっと楽になるが、 そのような計算法があるかどうかまではわからなかった。