掛ける順序が絡む場所だけ壊れる〜行列の指数関数と連立微分方程式

$\frac{dx}{dt}=ax$ の解は $x=e^{at}{x}_0$ です。未知数が並んで $\frac{d\mathbf{x}}{dt}=A\mathbf{x}$ になっても、$\mathbf{x}=e^{At}{\mathbf{x}}_0$ と書けます。

記号の見た目はそのまま。ただし $e^{A+B}=e^A{e}^B$ は成り立ちません^[4]。行列が掛ける順序を選ぶためです。

どこまでがスカラーと同じで、どこから違うのか。$e^{At}$ を実際に計算しながら見ていきます。

定義

スカラーの級数をそのまま持ちこみます^[1]。

$$e^A=\sum _{n=0}^{\infty }\frac{A^n}{n!}=I+A+\frac{A^2}{2!}+\frac{A^3}{3!}+\cdots$$

$A^0=I$ と約束します。無限和なので、まず収束を確かめる必要がある。

行列のノルムには $\Vert XY\Vert \le \Vert X\Vert \Vert Y\Vert$ という性質があります。作用素ノルムをとればこれが成り立つ。

$$\Vert \frac{A^n}{n!}\Vert \le \frac{\Vert A{\Vert }^n}{n!}$$

右辺を足すと $e^{\Vert A\Vert }$ で、有限の値です。したがって級数は絶対収束し、どんな正方行列でも $e^A$ が定まります^[2]。

$t$ を入れた $e^{At}$ についても同じで、すべての $t$ で収束します。

微分すると $A$ が出る

項別に微分します。

$$\frac{d}{dt}{e}^{At}=\sum _{n=1}^{\infty }\frac{n{t}^{n-1}{A}^n}{n!}=\sum _{n=1}^{\infty }\frac{t^{n-1}{A}^n}{(n-1)!}=A{e}^{At}$$

$A$ を前にくくり出しました。後ろにくくっても同じで、$e^{At}A$ とも書けます。$A$ と $A$ のべきは可換だからです。

これで $\mathbf{x}=e^{At}{\mathbf{x}}_0$ が ${\mathbf{x}}^{\prime }=A\mathbf{x}$ を満たすことが分かります。$t=0$ で $e^O=I$ なので、初期値も合う^[2]。

解がこれ 1 つであること

$\mathbf{x}(t)$ を任意の解とし、$\mathbf{y}(t)=e^{-At}\mathbf{x}(t)$ を考えます。

$${\mathbf{y}}^{\prime }=-A{e}^{-At}\mathbf{x}+e^{-At}{\mathbf{x}}^{\prime }=-A{e}^{-At}\mathbf{x}+e^{-At}A\mathbf{x}$$

$e^{-At}$ と $A$ は可換なので、2 項が打ち消して ${\mathbf{y}}^{\prime }=0$。$\mathbf{y}$ は定数ベクトルです。

$\mathbf{y}(0)=\mathbf{x}(0)={\mathbf{x}}_0$ なので $\mathbf{y}\equiv {\mathbf{x}}_0$。両辺に $e^{At}$ を掛けて $\mathbf{x}=e^{At}{\mathbf{x}}_0$ になります。

$(e^{At}{)}^{-1}=e^{-At}$ を使いました。$A$ と $-A$ は可換なので、この 1 本だけは安全に使える等式です。

$e^{A+B}=e^A{e}^B$ が崩れる例

次の 2 つを用意します。

$$A=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),B=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)$$

どちらも 2 乗すると零行列なので、級数は 2 項で止まります。

$$e^A{e}^B=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)\left(\begin{array}{cc}1& 0\\ 1& 1\end{array}\right)=\left(\begin{array}{cc}2& 1\\ 1& 1\end{array}\right)$$

一方 $A+B$ は 2 乗すると単位行列になります。偶数乗が $I$、奇数乗が $A+B$ なので、級数が双曲線関数に分かれる。

$$e^{A+B}=(\cosh 1)I+(\sinh 1)(A+B)=\left(\begin{array}{cc}1.543& 1.175\\ 1.175& 1.543\end{array}\right)$$

まったく別の行列です。$AB\ne BA$ であることが、そのまま食い違いになって出ました^[4]。

可換なら成り立ちます。証明では級数の積を展開して二項定理を使いますが、その並べ替えに $AB=BA$ が要ります。

計算法 1：対角化する

$A=PD{P}^{-1}$ と書ければ、$A^n=P{D}^n{P}^{-1}$ なので級数がそのまま移ります^[4]。

$$e^{At}=P{e}^{Dt}{P}^{-1}$$

$e^{Dt}$ は対角成分を指数関数にするだけ。固有値を並べれば終わりです。

例：対称行列で計算する

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

固有値は $3$ と $1$、固有ベクトルは $(1,1)$ と $(1,-1)$ です。

$$e^{At}=\frac{1}{2}\left(\begin{array}{cc}{e}^{3t}+e^t& e^{3t}-e^t\\ e^{3t}-e^t& e^{3t}+e^t\end{array}\right)$$

$t=0$ を入れると対角成分が $1$、非対角成分が $0$ になります。単位行列に戻り、確かに合っている。

計算法 2：スカラーとべき零に分ける

対角化できない場合の定番です。$A=\lambda I+N$ と分け、$N$ が何乗かで零行列になるようにします^[2]。

$\lambda I$ はどんな行列とも可換なので、指数を分けられます。

$$e^{At}=e^{\lambda t}\left(I+tN+\frac{t^2}{2!}{N}^2+\cdots \right)$$

$N^k=O$ なら、級数は $k$ 項で止まります。無限和が有限和になる。

例：ジョルダン細胞

$$A=\left(\begin{array}{cc}3& 1\\ 0& 3\end{array}\right)⟹e^{At}=e^{3t}\left(\begin{array}{cc}1& t\\ 0& 1\end{array}\right)$$

$3$ 次でも同じ要領です。

$$J=\left(\begin{array}{ccc}2& 1& 0\\ 0& 2& 1\\ 0& 0& 2\end{array}\right)⟹e^{Jt}=e^{2t}\left(\begin{array}{ccc}1& t& t^2/2\\ 0& 1& t\\ 0& 0& 1\end{array}\right)$$

$t$ の多項式が現れました。固有値が重なって固有ベクトルが足りないとき、この形になります。

計算法 3：ケイリー・ハミルトンを使う

$n$ 次の行列は自分の特性多項式を満たします。$A^n$ が $A^{n-1}$ 以下で書けるので、級数のすべての項が低次に落ちる。

$$e^{At}=c_0(t)I+c_1(t)A+\cdots +c_{n-1}(t)A^{n-1}$$

係数 $c_k(t)$ は、固有値ごとに次のスカラーの等式を解いて決めます。

$$e^{\lambda t}=c_0+c_1\lambda +\cdots +c_{n-1}{\lambda }^{n-1}$$

固有値が重なるときは、この等式を $\lambda$ で微分した式も使います。重根の数だけ式が増えます。

例：2 次でやってみる

さきほどの対称行列でもう一度計算します。固有値は $3$ と $1$ でした。

$$e^{3t}=c_0+3c_1,e^t=c_0+c_1$$

引き算して $c_1=\frac{e^{3t}-e^t}{2}$、代入して $c_0=\frac{3e^t-e^{3t}}{2}$ です。

$e^{At}=c_{0}I+c_{1}A$ の対角成分を計算します。

$$\frac{3e^t-e^{3t}}{2}+2\cdot \frac{e^{3t}-e^t}{2}=\frac{e^t+e^{3t}}{2}$$

対角化で出した答えと一致しました。固有ベクトルを求めずに済むのが、この方法の利点です。

計算法 4：ラプラス変換で出す

${\mathbf{x}}^{\prime }=A\mathbf{x}$ の両辺を変換すると $s\mathbf{X}-{\mathbf{x}}_0=A\mathbf{X}$ となり、$\mathbf{X}=(sI-A{)}^{-1}{\mathbf{x}}_0$。

つまり $(sI-A{)}^{-1}$ を逆変換すれば $e^{At}$ が出ます。

$$A=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)⟹(sI-A{)}^{-1}=\frac{1}{s^2+1}\left(\begin{array}{cc}s& 1\\ -1& s\end{array}\right)$$

各成分を逆変換します。$\frac{s}{s^2+1}$ が $\cos t$、$\frac{1}{s^2+1}$ が $\sin t$。

$$e^{At}=\left(\begin{array}{cc}\cos t& \sin t\\ -\sin t& \cos t\end{array}\right)$$

回転行列になりました。$A$ が反対称のとき、$e^{At}$ は直交行列になります。長さを変えない変換です。

相図で見る

$e^{At}$ が何をしているかは、平面上の点の動きを追うと見えます。固有値の種類で、軌道の形が変わる。

行列を 1 つ決めると、平面のすべての点の動きが同時に決まります。$e^{At}$ はその流れを $t$ 秒ぶん進める写像です。

使える等式と使えない等式

4 番目は見落とされがちですが、有用です。$\det e^A$ は指数関数の値なので決して $0$ にならない。$e^{At}$ が必ず正則であることが、これで分かります。

証明は三角化して見ます。$A$ を上三角に変換すると、対角には固有値 ${\lambda }_1,\dots ,{\lambda }_n$ が並ぶ。$e^A$ の対角は $e^{{\lambda }_1},\dots ,e^{{\lambda }_n}$ になり、行列式はその積です。

$$\det e^A=e^{{\lambda }_1}\cdots e^{{\lambda }_n}=e^{{\lambda }_1+\cdots +{\lambda }_n}=e^{\mathrm{tr}A}$$

右辺があるとき

${\mathbf{x}}^{\prime }=A\mathbf{x}+\mathbf{f}(t)$ を解きます。$e^{-At}$ を掛けて、積の微分にまとめる。

$$\frac{d}{dt}\left(e^{-At}\mathbf{x}\right)=e^{-At}{\mathbf{x}}^{\prime }-A{e}^{-At}\mathbf{x}=e^{-At}\mathbf{f}(t)$$

$0$ から $t$ まで積分し、$e^{At}$ を掛け戻します。

$$\mathbf{x}(t)=e^{At}{\mathbf{x}}_0+{\int }_0^t{e}^{A(t-\tau )}\mathbf{f}(\tau )d\tau$$

第 2 項をデュアメルの公式と呼びます^[2]。時刻 $\tau$ に加わった力が、残りの $t-\tau$ だけ流されて届く、という読み方ができます。

1 階のスカラー方程式で積分因子を使うのと、まったく同じ手順です。

例：一定の力を加える

回転を生む行列に、一定の力を加えます。初期値は $0$ とする。

$$A=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),\mathbf{f}=\left(\begin{array}{c}0\\ 1\end{array}\right)$$

$e^{A(t-\tau )}\mathbf{f}$ は、回転行列の第 2 列そのものです。

$$e^{A(t-\tau )}\mathbf{f}=\left(\begin{array}{c}\sin (t-\tau )\\ \cos (t-\tau )\end{array}\right)$$

$$\mathbf{x}(t)={\int }_0^t\left(\begin{array}{c}\sin (t-\tau )\\ \cos (t-\tau )\end{array}\right)d\tau =\left(\begin{array}{c}1-\cos t\\ \sin t\end{array}\right)$$

ばねに一定の力を加えたときの応答です。中心がずれた円を描きます。

級数をそのまま計算してはいけない

定義どおりに級数を打ち切る方法は、数値計算では使えません^[3]。

$\Vert A\Vert$ が大きいと項数が足りず、しかも途中の項が巨大になって桁落ちを起こします。答えが小さくても、計算の途中で $10^{16}$ を超える値が出ることがある。

固有ベクトルを使う方法も安全ではありません。固有値が近いと固有ベクトルの行列がほとんど特異になり、$P^{-1}$ で誤差が拡大します。

実際に使われるのはスケーリングと 2 乗です^[3]。

$$e^A={\left(e^{A/2^m}\right)}^{2^m}$$

$m$ を大きくして $\Vert A/2^m\Vert$ を小さくし、そこでパデ近似を当てる。あとは 2 乗を $m$ 回くり返します。

正確な式が欲しいのか、数値が欲しいのか。目的で方法が変わります。

物理での顔

量子力学の時間発展は $U(t)=e^{-\mathrm{i}Ht/\hslash }$ で書けます。$H$ がエルミートなら $U$ はユニタリになる。

確かめは短くて、$U^{†}U=e^{\mathrm{i}{H}^{†}t/\hslash }{e}^{-\mathrm{i}Ht/\hslash }$ です。$H^{†}=H$ なので指数の肩が打ち消し合い、$I$ になります。

指数の肩どうしが可換なので、ここでは $e^A{e}^B=e^{A+B}$ を使ってよい。可換性が効く場面です。

制御工学では $e^{At}$ を状態遷移行列と呼びます。$\mathrm{Re}\lambda <0$ がすべての固有値で成り立てば、系は原点へ収束する。

級数で定義し、微分すると $A$ が出る。ここまではスカラーと同じで、崩れるのは掛ける順序が絡む場所だけです。

参考