行列計算を使わない線形代数 #11 〜 おまけ：線形常微分方程式の解（行列の指数関数とLie群の視点から）

今回は線形常微分方程式の解を取り上げます。定係数の場合は、よく知られているように行列の指数関数で計算が可能になります。しかし、係数行列が時間依存する場合はその方法では解は計算できません。計算方法としては、逐次近似法などいくつかありますが、ここではLie群の視点から解の級数展開を考えます。
ただし、すべての用語を説明していると書ききれないので、たいぶ省略した説明になってます。

$${n}$$次元ユークリッド空間$${\mathbb{R}^n}$$上の微分方程式

$$
\displaystyle
\frac{d\bm{x}}{dt} = A(t) \bm{x}, \quad \bm{x}(0)=\bm{x}_0, \quad \quad \quad (1) 
$$

を考えます。ここで、$${A(t) \in \mathbb{R}^{n\times n},t\in\mathbb{R}, }$$は$${n}$$次正方行列の1径数族（one-parameter family）で、$${t}$$に関して滑らかであるとします。

いま、$${\mathbb{R}^n}$$の標準基底$${\bm{e}_1,\cdots, e_n}$$を初期値$${\bm{x}(0)}$$とする解を、それぞれ$${\bm{u}_1, \cdots, \bm{u}_n}$$とします。つまり、$${d\bm{u}_k/dt = A \bm{u}_k, \bm{u}_k(0)=\bm{e}_k}$$です。$${\bm{u}_1, \cdots, \bm{u}_n}$$をベクトル値関数とみなして、並べて作られる行列値関数$${X(t) = (\bm{u}_1(t), \cdots, \bm{u}_n(t))}$$を考えると、$${X(t)}$$は以下の微分方程式を満たします：

$$
\displaystyle
\frac{dX(t)}{dt} = A(t)X(t) , \quad X(0) = I_n, \quad \quad \quad (2).  
$$

ここで、$${I_n}$$は$${n}$$次の単位行列です。つまり、微分方程式として(1)と(2)は等価になるので、以降では(2)を扱っていきます。

(2)の解をどのようにして求めたらいいでしょうか？まずは常微分方程式の解の存在を示すときに利用する、Picard（ピカール）の逐次近似法を考えてみましょう。逐次近似法を(2)に適用すると、関数列

$$
\displaystyle 
X_{0}(t) = I_{n}, \quad X_{k}(t) = \int_{0}^{t} A(\tau)X_{k-1}(\tau)d\tau ,
$$

を考えることになります。これを計算することで、

$$
\displaystyle
X_{k}(t) = I_n+ \sum_{m=1}^{k}\int_{0}^{t} A(t_{m}) \int_{0}^{t_{m}} A(t_{m-1}) \cdots\int_{0}^{t_2} A(t_{1}) dt_{1}\cdots dt_{m-1}dt_{m}
$$

となり、(2)の解は$${\displaystyle X(t)=\lim_{k\to\infty}X_{k}(t)}$$として得ることができます。ここでの収束はPicardの逐次近似法の証明で保証されています。

ただ、この関数列では収束するわけですが、もう少し別の表示を考えてみることにしましょう。そうすることで、その背景にある理論が少し見えてきます。話を進める前にいくつか記号を導入しておきます。

$$
G:=GL_{+}(n,\mathbb{R}) := \{ X \in \mathbb{R}^{n\times n} \,\,|\,\, \det A > 0 \},
$$

$$
\frak{g} := \frak{gl}(n,\mathbb{R}) := \mathbb{R}^{n\times n},
$$

とします。行列式$${\det A}$$は本編ではまだ出てきていないですが、これはおまけということで許してください。

■ケース1：定数行列の場合

最初に$${A(t)}$$が時間$${t}$$に依存しない場合を考えましょう。いま、$${A(t)\equiv A}$$と書きます。このとき、上記の$${X_k(t)}$$は

$$
\displaystyle
X_{k}(t) = I_n + \sum_{m=1}^{k}A^k\int_{0}^{t} \int_{0}^{t_{m}} \cdots\int_{0}^{t_2} 1 dt_{1}\cdots dt_{m-1}dt_{m}
$$

となるので、積分を実行し$${k\to\infty}$$とすることで、(2)の解は$${t}$$のべき級数として

$$
\displaystyle 
X(t)  = \sum_{k=0}^{\infty}\frac{t^k}{k!} A^k= I_n + tA + \frac{t^2}{2!}A^2 + \frac{t^3}{3!}A^3 + \cdots , \quad \quad (3), 
$$

で得られます。先程も述べたように、この級数は逐次近似法の証明で収束が保証されていますが、行列のノルム$${ \| A \|^2 = \mathrm{tr}(A^T A)}$$を導入することで、このべき級数が$${t}$$に関して一様収束することを直接的に示せます。

定数行列$${A\in \frak{g}}$$に対して、(3)のべき級数で定まる行列を$${e^{tA}}$$と書き、行列の指数関数と言います。行列の指数関数$${X(t)=e^{tA}}$$は、1変数の指数関数$${e^{at}, (a\in\mathbb{R}\text{は定数}),  }$$と同じく、微分に関して

$$
\displaystyle
\frac{dX}{dt} = AX = XA
$$

が成り立ちます。また、$${\det e^{tA}=e^{t\,\mathrm{tr}(A)}>0}$$が成り立つので、$${e^{tA}\in G}$$になります。つまり、(2)の解は$${X(t) = e^{tA}}$$であって、(2)は$${G}$$上の微分方程式を与えているということが言うことができます。

さらに$${\displaystyle \frac{dX}{dt}(0)=A}$$であるので、$${G}$$上の曲線$${\mathbb{R}\to G ; t\mapsto A(t)}$$上の点$${X(0)=I_n}$$の接ベクトルを行列$${A\in\frak{g}}$$が与えています。つまり、$${G}$$の点$${I_n\in G}$$における接ベクトルのなすベクトル空間（接空間：tangent space）$${T_{I_n}G}$$は、$${\frak{g}}$$と同一視できます。したがって、行列の指数関数は以下のような写像として考えることができます：

$$
\exp \,\, : \,\, T_{I_n}G\cong\frak{g} \to G \,\, ;\,\, A \mapsto e^{A}. 
$$

この写像はLie群$${G}$$の指数写像（exponential map）と呼ばれます。

■ケース2：定数行列ではない場合

次に$${A(t)}$$が$${t}$$に関して一定ではない場合を考えます。ケース1の議論を延長することで、$${A(t)}$$が$${t}$$に関して一定でない場合でも、微分方程式(2)は$${G}$$の微分方程式としてみなすことができます。つまり、$${X(t)\in G.}$$ しかし、状況はケース1よりも複雑になります。

まず、Picardの逐次近似法で導出した列$${X_k}$$を考えます。ここでの難点は、必ずしも$${X_k(t)\in G, ( t \in\mathbb{R}) }$$とは限らないということです。つまり、Picardの逐次近似法では$${G}$$上の微分方程式(2)の$${G}$$上の近似列を与える訳ではありません。

次に定数行列の場合のように、行列の指数関数で解を表現することを考えてみましょう。しかし、ここでも難点があります。それは定数行列である場合を一般化したような$${\displaystyle X(t)=e^{\int_{0}^{t}A(\tau)d\tau}}$$が、(2)の解を与えないのです。実際に$${\displaystyle \Sigma(t)=\int_{0}^{t}A(\tau)d\tau}$$とおいて、$${X(t)=e^{\Sigma(t)}}$$の微分を考えると、

$$
\displaystyle 
\frac{de^{\Sigma}}{dt} =\sum_{k=0}^{\infty}\frac{1}{n!}\left(\frac{d\Sigma}{dt}\Sigma^{n-1}+\Sigma\frac{d\Sigma}{dt}\Sigma^{n-2}+\cdots+\Sigma^{n-1}\frac{d\Sigma}{dt}\right)
$$

となるので、$${\displaystyle \Sigma(t)=\int_{0}^{t}A(\tau)d\tau}$$と$${\displaystyle\frac{d\Sigma(t)}{dt}=A(t)}$$が交換可能でなければ、$${e^{\Sigma(t)}}$$は微分方程式(2)を満たしません。ここで交換可能であるとは、

$$
\displaystyle
\left[ \Sigma, \frac{d\Sigma}{dt}\right] := \Sigma \frac{d\Sigma}{dt}-\frac{d\Sigma}{dt}\Sigma = 0
$$

を満たすことを指します。

そこで、(2)の解が$${X(t)=e^{\Omega(t)}}$$であると仮定して、$${\Omega}$$が満たすべき条件を考えてみることにしましょう。まず、$${X(t)=e^{\Omega(t)}}$$が(2)の解であることから、

$$
\displaystyle
\frac{d e^{\Omega(t)}}{dt} e^{-\Omega(t)} = A(t), \quad \Omega(0)=O, \quad\quad\quad (4) 
$$

を満たします。ただ、この式では$${d\Omega/dt}$$は直接には現れてきません。そこで次のようなトリックを使って、$${d\Omega/dt}$$が現れる式を導出します。

いま、

$$
\displaystyle U(s,t) := e^{-s\Omega(t)}\frac{\partial}{\partial t} e^{s\Omega(t)}, \quad s,t \in\mathbb{R}, 
$$

とおきます。この$${U(s,t)}$$の変数$${s}$$による偏微分を考えると、

$$
\displaystyle \begin{align*}\frac{\partial U}{\partial s} &= e^{-s\Omega(t)}(-\Omega)\frac{\partial}{\partial t} e^{s\Omega(t)} + e^{-s\Omega(t)}\left(\frac{d\Omega}{dt} e^{s\Omega(t)}+\Omega\frac{\partial}{\partial t} e^{s\Omega(t)} \right)\\ &= e^{-s\Omega(t)}\frac{d\Omega}{dt} e^{s\Omega(t)}\end{align*}
$$

となります。したがって、$${U(0,t)=O}$$なので、

$$
\displaystyle 
U(1,t) = \int_{0}^{1} \frac{\partial U}{\partial s} ds = \int_{0}^{1} e^{-s\Omega(t)}\frac{d\Omega}{dt} e^{s\Omega(t)} ds. 
$$

よって、$${U(s,t)}$$の定義と(4)を使うことで、

$$
\displaystyle 
A(t) =  \int_{0}^{1} e^{s\Omega(t)}\frac{d\Omega}{dt} e^{-s\Omega(t)} ds. 
$$

$${e^{\pm s\Omega(t)}}$$をべき級数に展開して、整理すると、

$$
\displaystyle 
A(t) =  \int_{0}^{1} \sum_{k=0}^{\infty}\frac{s^{k}}{k!}\sum_{\ell=0}^{\infty} \frac{k!}{\ell!(k-\ell)!} \Omega^{\ell} \frac{d\Omega}{dt} (-\Omega)^{k-\ell} ds
$$

となります。さらに、上式の右辺は次のように表すことができます（演習問題【2】も参照）。

$$
\displaystyle 
A(t) =  \sum_{k=0}^{\infty}\frac{(\mathrm{ad}_{\Omega})^{k}}{(k+1)!}\left( \frac{d\Omega}{dt} \right). \quad \quad \quad (5) 
$$

ここで、$${\mathrm{ad}_{\Omega} : \frak{g}\to\frak{g}}$$は$${\Omega\frak{g}}$$の随伴作用と呼ばれ、$${\mathrm{ad}_{\Omega} (Y) := [\Omega,Y] = \Omega Y - Y\Omega}$$で定義されます。さらに、

$$
(\mathrm{ad}_{\Omega})^{k}(\bullet) := \underbrace{(\mathrm{ad}_{\Omega}\circ\cdots\circ\mathrm{ad}_{\Omega})}_{k\,\text{回}}(\bullet) = \underbrace{[\Omega, [\Omega, [ \cdots ,[\Omega}_{k\, \text{回}}, \bullet]]]]
$$

です。(5)の右辺は形式的に$${\displaystyle \frac{e^{\mathrm{ad}_{\Omega}}-1}{\mathrm{ad}_{\Omega}}\left(\frac{d\Omega}{dt}\right)}$$とも書かれます。これは1変数関数$${(e^{x}-1)/x}$$のべき級数展開が$${\sum_{k=0}^{\infty}x^k/(k+1)!}$$となることからです。

さて、$${(e^{\mathrm{ad}_{\Omega}}-1)/{\mathrm{ad}_{\Omega}}}$$は次のような線形写像と考えることができます：

$$
\displaystyle 
\frak{g}\rightarrow\frak{g} \,\,;\,\, Y \mapsto \frac{e^{\mathrm{ad}_{\Omega}}-1}{\mathrm{ad}_{\Omega}}(Y)= \sum_{k=0}^{\infty}\frac{(\mathrm{ad}_{\Omega})^{k}}{(k+1)!}( Y ). 
$$

長くなるので証明は省略しますが、この線形写像は次のような逆写像を持つことがわかっています：

$$
\displaystyle
\left( \frac{e^{\mathrm{ad}_{\Omega}}-1}{\mathrm{ad}_{\Omega}} \right)^{-1} = \sum_{k=0}^{\infty}\frac{B_k}{k!}(\mathrm{ad}_\Omega)^{k}.
$$

ここで、$${B_k, k=0,1,2,\cdots,}$$はベルヌーイ数と呼ばれる実数で、1変数関数$${x/(e^x-1)}$$のTaylor展開することで以下のように定義されます：

$$
\displaystyle 
\frac{x}{e^{x}-1} = \sum_{k=0}^{\infty}\frac{B_k}{k!}x^k = 1 -\frac{1}{2}x +\frac{1}{6}x^{2} - \frac{1}{30}x^{4} + \cdots .
$$

この事実を認めたとすると、$${d\Omega/dt}$$は(5)より次のように書き下すことができます：

$$
\displaystyle \begin{align*}
\frac{d \Omega }{dt} &=  \sum_{k=0}^{\infty}\frac{B_k}{k!}(\mathrm{ad}_{\Omega})^{k}( A) \\ &= A - \frac{1}{2}[\Omega, A] +\frac{1}{12}[\Omega,[\Omega,A]]-\frac{1}{720}[\Omega,[\Omega, [\Omega,A]]] + \cdots . \end{align*}
$$

ようやく$${\Omega(t)}$$が満たすべき微分方程式を（級数の形として）導くことができました。この微分方程式を解くために逐次近似法をもう一度使ってみましょう。つまり、関数列$${\Omega_{k}(t)}$$を

$$
\displaystyle
\begin{align*}\Omega_0(t) & \equiv O \\
\Omega_{k+1}(t) & := \sum_{\ell=0}^{\infty}\frac{B_\ell}{\ell !} \int_{0}^{t}(\mathrm{ad}_{\Omega_k(\tau)})^\ell(A(\tau))d\tau.\quad\quad\quad (6) \end{align*}
$$

で定義します。この計算を実行して$${k\to\infty}$$とすることで、$${\Omega(t)}$$を得ることができました：

$$
\displaystyle
\begin{align*}& \Omega(t) \\ = &\int_{0}^{t}A(t_1)dt_1 + \frac{1}{2}\int_{0}^{t}dt_1\int_{0}^{t_1}dt_2 [A(t_1), A(t_2)] \\ &+\frac{1}{6}\int_{0}^{t}dt_1 \int_{0}^{t_1} dt_2 \int_{0}^{t_2}dt_3 \Big( [A(t_1), [A(t_2), A(t_3)]] + [[A(t_1), A(t_2)], A(t_3)] \Big) \\ &+ \cdots .  \end{align*}
$$

だいぶ長くなってしまいましたが、当初考えていた微分方程式(2)の解は、上式を満たす$${\Omega(t)}$$を用いて、$${X(t)=e^{\Omega(t)}}$$として得ることができました。この解の表示をマグヌス展開（Magnus expansion）と呼びます。ただし、$${\Omega}$$の収束性は議論していませんでした。さすがに収束性までは書ききれないので、詳細は関連の文献などを見てください。

参考文献：S. Blanes, F. Casas, J.A. Oteo, J. Ros, "The Magnus expansion and some of its applications", preprint, https://arxiv.org/abs/0810.5488.

本文は以上で、記事自体は完結しています。演習問題以降は有料になってますが、支援いただけるという意味で、この記事を購入いただけましたら幸いです。

■演習問題