【微積物理入門】行列を用いた微分方程式の解法

いよいよ行列を用いて微分方程式を解いていく。といっても、現行の学習指導要領では行列を扱わないらしいので次のような執筆方針にしたい。

まずは行列を用いると微分方程式が無理なく”導出”できることを示す
次に線形代数学の基本的な知識について執筆する

線形代数学といっても、これは奥が深い学問である。恥ずかしながら私も”理解した”とは到底言えない。しかし微分方程式を解くのに必要な知識は限られている。本書の目的は微分方程式を解くことで物理を理解することにあるから、必要な数学を提示するにとどめる。

ここでのゴールは固有値問題である。のちに丁寧に解説することにしよう。

行列を用いた解法

二階微分方程式

$\ddot{x} = x$

を解いていく。しかしこのままではうまく扱えないので、速度 $v=\dot{x}$ を挟むことで連立微分方程式に変換することができる。

補助変数（速度）を導入 $v=\dot{x}$

この補助変数を導入すると、もとの微分方程式は次のように表現できる。

$\dot{v} = x$

ところで $v=\dot{x}$ の関係もあるので、整理すると

\begin{eqnarray} \left\{\begin{array}{l}\dot{x} = v \\ \dot{v} = x \end{array} \right. \end{eqnarray}

このような連立方程式を解くことと同じことになる。これは、ベクトルとして書くと：

\begin{eqnarray} \frac{d}{dt} \begin{bmatrix} x \\ v \end{bmatrix} = \begin{bmatrix}0 & 1 \\1 & 0\end{bmatrix}\begin{bmatrix}x \\v \end{bmatrix} \end{eqnarray}

行列表現

ここで次のベクトル $\mathbf{u}$ と行列 $C$ を導入する。 \begin{eqnarray} \mathbf{u} = \begin{bmatrix} x \\ v \end{bmatrix} \end{eqnarray}

\begin{eqnarray} C = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \end{eqnarray}

するとこの微分方程式は次のように書くことができる。

\begin{equation} \frac{d\mathbf{u}}{dt} = C\mathbf{u} \end{equation}

一般解：指数関数解の形

この微分方程式は、ベクトルの形をしているが次のような微分方程式に非常に似ている。

\begin{equation} \frac{dx}{dt} = cx \quad(cは定数) \end{equation}

この微分方程式の一般解は、変数分離を用いて次の一般解を得ることができる。（演習問題にちょうどいいので各自ためしに解いてほしい。）

\begin{equation} x = c_0\exp(ct) \end{equation}

$c_0$ はこの微分方程式の積分定数であるが、 $t=0$ のときの $x$ の値を $x_0$ とすると $c_0 = x_0$ が成り立つので、この値を使おう。

\begin{equation} x = x_0\exp(ct) \end{equation}

では行列を用いた解法に戻ろう。形が似ているのだから、同じ一般解をもつはずである。実際、いま問題としている微分方程式の一般解は指数解をもち

\begin{equation} \mathbf{u}(t) = \exp(Ct) \mathbf{u}(0) \end{equation}

と書くことができる。したがって、必要なのは $\exp{Ct}$ を計算することだ。

行列指数関数の計算

そもそも行列もわからないのに行列の指数関数とは何か？となると混乱するかもしれない。この $\exp(C)$ は、いままでの指数関数の感覚では素朴に定義できない。しかし次の公式（テイラー展開）：

\begin{equation} \exp(c) = \sum _{k=0} ^\infty \frac{c ^k}{k!} = 1+c+\frac{c ^2}{2!} + \frac{c ^3}{3!} + \dots + \frac{c ^k}{k!} + \dots \end{equation}

を用いると、次のように定義することができる。

\begin{equation} \exp(C) = \sum _{k=0} ^\infty \frac{C ^k}{k!} = 1+C+\frac{C ^2}{2!} + \frac{C ^3}{3!} + \dots + \frac{C ^k}{k!} + \dots \end{equation}

つまり行列の $k$ 乗、 $C ^k$ さえ求めることができればよいことがわかる。ところで今扱っている行列 $C$ は次のような形をしていた。

\begin{eqnarray} C = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \end{eqnarray}

これを2乗、3乗、4乗……と求めて足していくのは非常に大変、というより不可能である。しかし運のいいことに、行列指数関数には、次のような性質がある：

対角行列

\begin{eqnarray} D = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \end{eqnarray}

に対して

\begin{eqnarray} \exp(D) = \begin{bmatrix} \exp(\lambda_1) & 0 \\ 0 & \exp(\lambda_2) \end{bmatrix} \end{eqnarray}

と書くことができる。そこで行列 $C$ が $D$ のような対角行列であればいいのだが、実は対角化という手法によってこの形に変換することができる。対角化というのは次の等式を満たす行列 $P$ をみつけることに他ならない。

\begin{equation} C = P D P ^{-1} \end{equation}

この式を $\exp(C)$ に代入しても $\exp( P D P ^{-1})$ となって余計にややこしくなっただけのように見えるかもしれない。しかしまたしても運のいいことにこの式は

\begin{equation} \exp(C) =\exp( P D P ^{-1}) = P \exp(D) P ^{-1} \end{equation}

と変形することができる。これなら計算することができる。したがって、行列指数関数を計算するために必要な数学的技法は対角化にあるということがわかった。

対角化の方法については次の記事で触れた。やりかたがわからない人は参考にしてほしい。

bisebu.jp

連立微分方程式の一般解

話を元に戻そう。もともとは次の方程式、

\begin{equation} \frac{d\mathbf{u}}{dt} = C\mathbf{u} \end{equation}

の一般解が、

\begin{equation} \mathbf{u}(t) = \exp(Ct) \mathbf{u}(0) \end{equation}

と書けるという話だった。ここで $\exp(Ct)$ は次のように対角化することができる。

\begin{equation} \exp( P ^{-1} Ct P) = P ^{-1} \exp(Ct) P = \exp(Dt) \end{equation}

よって

\begin{equation} \exp(Ct) = P \exp(Dt) P ^{-1} \end{equation}

であるから、

\begin{equation} \mathbf{u}(t) = P \exp(Dt) P ^{-1} \mathbf{u}(0) \end{equation}

が得られる。やはり $\exp(Ct)$ を解くのが肝のようだ。

具体的な微分方程式を解く

今回解いている問題は、二階微分方程式

$\ddot{x} = x$

であった。これを連立微分方程式に直すと

\begin{equation} \frac{d\mathbf{u}}{dt} = C\mathbf{u}, \quad \mathbf{u} = \begin{bmatrix} x \\ v \end{bmatrix}, \quad C = \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} \end{equation}

となる。一般解を求めるために、行列 $C$ の対角化を試みる。固有方程式は、

\begin{equation} \begin{vmatrix} \lambda & -1 \\ -1 & \lambda \end{vmatrix} = 0 \end{equation}

よって $\lambda ^2 -1 = 0 \therefore \lambda = \pm 1$ を得る。 $\lambda_1 = 1$ に対応する固有ベクトルは $\mathbf{p_1} = (1,1)$ で、固有値 $\lambda_2 = -1$ に対応する固有ベクトルは $\mathbf{p_2} = (1,-1)$ である（ $\mathbf{p_2} = (-1,1)$ としてもよい）。対角化行列の逆行列は

\begin{equation} P ^{-1} = \frac{1}{2} \begin{bmatrix} -1 & 1 \\ 1 & 1 \end{bmatrix} \end{equation}

対角化した結果は

\begin{equation} P ^{-1} C P = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \end{equation}

となるはずなので、ここから、

\begin{equation} C = P \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} P ^{-1} =\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}\frac{1}{2} \begin{bmatrix} -1 & 1 \\ 1 & 1 \end{bmatrix} \end{equation}

よって

\begin{equation} \exp(Ct) = P \exp(Dt) P ^{-1} = \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} \exp(t) & 0 \\ 0 & \exp(-t) \end{bmatrix} \frac{1}{2}\begin{bmatrix} -1 & 1 \\ 1 & 1 \end{bmatrix} \end{equation}

計算して、

\begin{equation} \exp(Ct) =\frac{1}{2}\begin{bmatrix} -\exp(t)+\exp(-t) & \exp(t)+\exp(-t)\\ -\exp(t)-\exp(-t) & \exp(t)-\exp(-t) \end{bmatrix}\end{equation}

結果的にこの微分方程式の解は、

\begin{eqnarray} \begin{bmatrix} x(t) \\ v(t) \end{bmatrix} = \begin{bmatrix} -\exp(t)+\exp(-t) & \exp(t)+\exp(-t) \\ -\exp(t)-\exp(-t) & \exp(t)-\exp(-t) \end{bmatrix} \begin{bmatrix} x(0) \\ v(0) \end{bmatrix} \\ = \begin{bmatrix} \left(\frac{-\exp(t)+\exp(-t)}{2} \right) x(0) + \left(\frac{\exp(t)+\exp(-t)}{2}\right) v(0) \\ \left(\frac{-\exp(t)-\exp(-t)}{2} \right) x(0) + \left(\frac{\exp(t)-\exp(-t)}{2}\right) v(0) \end{bmatrix} \end{eqnarray}

これでもいいのだが、すこし式変形しよう。 $\exp(t) と \exp(-t)$ でそれぞれくくると、

\begin{equation} \begin{bmatrix} x(t) \\ v(t) \end{bmatrix}
= \begin{bmatrix} \left(\frac{-x(0)+v(0)}{2} \right) \exp(t) + \left(\frac{x(0)+v(0)}{2}\right) \exp(-t) \\ \left(\frac{-x(0)+v(0)}{2} \right) \exp(t) - \left(\frac{x(0)+v(0)}{2}\right) \exp(-t) \end{bmatrix} \end{equation}

$x(0),v(0)$ は任意定数であるから、係数をそれぞれ $A,B$ とまとめると、

\begin{equation} \begin{bmatrix} x(t) \\ v(t) \end{bmatrix}
= \begin{bmatrix} A \exp(t) + B \exp(-t) \\ A \exp(t) - B \exp(-t) \end{bmatrix} \end{equation}

となる。求めたかった $x(t)$ だけ取り出すと、

\begin{equation} x(t) = A \exp(t) + B \exp(-t) \end{equation}

が得られた。

今後の方針

ここまで本当に大変な道のりだった。しかし微分方程式が天下り式ではなく一歩ずつ順を追って解けることが確認できただろう。特に $\exp(t) と \exp(-t)$ の定数倍と和で書けることがわかったのは大きな成果である。なぜかというと、どんな理由であれ答え（のようなもの）を必要な個数見つけさえしてしまえばあとは係数を調整するだけで済むからだ。

もちろん、実際には一筋縄に行かない微分方程式もたくさんある……というよりそちらのほうがはるかに多い。例えば特性方程式が重解を持った場合は一歩難しくなる。このような場合について解説する必要でてきたら、天下り式の解き方と丁寧な導出を併記することにする。

それでもまだ「なぜ特性方程式の方法でうまくいくのか」に疑問を持つ人もいるかもしれない。（そういう人は物理学科より数学科のほうが向いている気がする。）本書はあくまでも微積物理の本であるから深入りはしたくないのだが、もし時間があれば線形代数学の少し深い部分について触れることもあるかもしれない。