連立微分方程式の3通りの解き方

更新 2022/03/06

この記事では，(1階の)連立線形微分方程式の解法を3通り紹介します。

連立微分方程式

例題

次の連立微分方程式を解け。 $\begin{cases} \dfrac{dx_1(t)}{dt} = x_1(t)-2x_2(t)\\ \dfrac{dx_2(t)}{dt} = x_1(t) +4x_2(t) \end{cases}$ ただし $x_1(0)=1,x_2(0)=-2$ とする。

2つの未知関数 $x_1(t),x_2(t)$ を求める問題です。後半では未知関数 $n$ 個の場合についても述べます。

解き方1：未知関数を1つ消去する方法

まずは， $x_1(t)$ ， $x_2(t)$ のうち1つを消去する方法です。

連立微分方程式の解き方1

$\begin{cases} \dfrac{dx_1}{dt} &= ax_1 + bx_2\\ \dfrac{dx_2}{dt} &= cx_1 + dx_2 \end{cases}$

に対して， $x_1,x_2$ の片方を消去すると，もう片方についての二階線形微分方程式に帰着できる。

例題を解いてみます。

例題の解答

$x_1$ を消去する。2つめの式より $x_1 = \dfrac{dx_2}{dt} - 4x_2$ であり，これを1つめの式
$\dfrac{dx_1}{dt} = x_1-2x_2$

に代入して整理することで，2階の線形微分方程式 $\dfrac{d^2x_2}{dt^2} -5 \dfrac{dx_2}{dt} +6x_2 = 0$ が得られる。特性方程式 $\lambda^2-5\lambda+6=0$ の解は $\lambda=2,3$ である。よって，解は定数 $C_1 , C_2$ を用いて $x_2 = C_1 e^{2t} + C_2 e^{3t}$ と表せる(→補足)。これを元の式に代入すると $x_1 = -2C_1 e^{2t} - C_2 e^{3t}$ となる。初期条件を代入すると， $\begin{cases} -2C_1 - C_2 = 1\\ C_1 + C_2 = -2 \end{cases}$ これを解き $C_1 = 1, C_2 = -3$

以上より $x_1 = -2e^{2t} + 3e^{3t},x_2=e^{2t} -3e^{3t}$

補足：連立ではない「線形微分方程式の解き方」を使いました。 →微分方程式の解法(同次形・線形微分方程式)

解き方2：良い関数を探す方法

「微分しても元の関数の定数倍」になるような良い関数を探します。

連立微分方程式の解き方2

$\dfrac{d}{dt}(x_1-kx_2) = r(x_1 - kx_2)$

となるような $k,r$ が2ペア見つかると解ける。

例題を解いてみましょう。

冒頭の例題の別解

$\begin{cases} \dfrac{dx_1(t)}{dt} = x_1(t)-2x_2(t)\\ \dfrac{dx_2(t)}{dt} = x_1(t) +4x_2(t) \end{cases}$

より，

$\dfrac{d}{dt}(x_1-kx_2)\\ =x_1-2x_2-k(x_1+4x_2)\\ =(1-k)x_1+(-2-4k)x_2$

よって， $1:(1-k)=-k:(-2-4k)$ を解くと $k = -1,-2$ となる。つまり，元の微分方程式から $\begin{cases} \dfrac{d}{dt} (x_1+x_2) = 2(x_1+x_2)\\ \dfrac{d}{dt} (x_1+2x_2) = 3(x_1+2x_2) \end{cases}$ が得られる。それぞれ解くと，定数 $C_1,C_2$ を用いて $x_1+x_2 = C_1 e^{2t} , x_1+2x_2 = C_2 e^{3t}$ となる。初期条件を代入すると， $C_1 = -1,C_2=-3$ となる。 $x_1,x_2$ をそれぞれ消去すると， $x_1 = -2e^{2t} + 3e^{3t}, x_2=e^{2t}-3e^{3t}$

なお， $k$ に関する2次方程式が重解を持つ場合はやや大変なので，さきほどの方法1をオススメします。

連立微分方程式と連立漸化式

3つめの解き方の前に，おもしろい余談です。

実は，線形連立微分方程式の解法は，連立漸化式の解法とほとんど同じです。連立漸化式の3通りの解き方では，以下の3つの方法を紹介しました：
1. 三項間漸化式に帰着させる方法
2. 等比数列を作る方法
3. 行列の $n$ 乗を用いる方法
この記事では，線形連立微分方程式について，この3つに対応する解法をそれぞれ紹介しています。
方法1と2は（そのままでは）2変数のときにしか使えませんが，方法3は $n$ 変数の連立微分方程式・連立漸化式で使えます。

解き方3：行列を用いる方法

3つめの解き方は，一般の nnn 変数の場合で説明します。以下の1階線形連立微分方程式を考えます：
{dx1dt=a11x1+a12x2+⋯+a1nxndx2dt=a21x1+a22x2+⋯+a2nxn⋮dxndt=an1x1+an2x2+⋯+annxn\begin{cases}
\dfrac{dx_1}{dt} = a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n\\
\dfrac{dx_2}{dt} = a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n\\
\quad\vdots\\
\dfrac{dx_n}{dt} = a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n
\end{cases}⎩⎨⎧​dtdx1​​=a11​x1​+a12​x2​+⋯+a1n​xn​dtdx2​​=a21​x1​+a22​x2​+⋯+a2n​xn​⋮dtdxn​​=an1​x1​+an2​x2​+⋯+ann​xn​​
この式を次のように書き換えます。
ddt(x1x2⋮xn)=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)(x1x2⋮xn)
\dfrac{d}{dt} \begin{pmatrix} x_1\\x_2\\ \vdots \\x_n \end{pmatrix}=
\begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn} \end{pmatrix}
\begin{pmatrix} x_1\\x_2\\ \vdots \\x_n \end{pmatrix}
dtd​⎝⎛​x1​x2​⋮xn​​⎠⎞​=⎝⎛​a11​a21​⋮an1​​a12​a22​⋮an2​​⋯⋯⋱⋯​a1n​a2n​⋮ann​​⎠⎞​⎝⎛​x1​x2​⋮xn​​⎠⎞​
ここで x=(x1x2⋮xn)\mathbb{x} = \begin{pmatrix} x_1\\x_2\\ \vdots \\x_n \end{pmatrix}x=⎝⎛​x1​x2​⋮xn​​⎠⎞​，A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)A = \begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn} \end{pmatrix}A=⎝⎛​a11​a21​⋮an1​​a12​a22​⋮an2​​⋯⋯⋱⋯​a1n​a2n​⋮ann​​⎠⎞​ とおくと，
ddtx=Ax
\dfrac{d}{dt} \mathbb{x} = A \mathbb{x}
dtd​x=Ax
と書けます。実は，上の解は x=etAc\mathbb{x} = e^{tA} \mathbb{c}x=etAc になることが知られています。ただし，etAe^{tA}etA は行列の指数関数です。詳細は，以下を参照してください。
行列の指数関数とその性質
行列の指数関数の計算法
また，c\mathbb{c}c は要素数 nnn の行のベクトルです。初期値が x0\mathbb{x}_0x0​ で与えられている場合，x=etAx0\mathbb{x} = e^{tA} \mathbb{x}_0x=etAx0​ となります。
行列 AAA が対角化できる場合は etAe^{tA}etA の計算は簡単です。対角行列 DDD を用いて A=PDP−1A = PDP^{-1}A=PDP−1 と書けるので，etA=PetDP−1e^{tA} = P e^{tD} P^{-1}etA=PetDP−1 となります。
冒頭の例題の別解
与式は
ddt(x1x2)=(1−214)(x1x2)
\dfrac{d}{dt} \begin{pmatrix}x_1\\x_2\end{pmatrix}
= \begin{pmatrix} 1&-2\\1&4\end{pmatrix} \begin{pmatrix} x_1\\x_2 \end{pmatrix}
dtd​(x1​x2​​)=(11​−24​)(x1​x2​​)
と変形できる。上の式の 2×22\times 22×2 行列を AAA とおく。AAA は行列 P=(2−1−11)P = \begin{pmatrix} 2&-1\\-1&1\end{pmatrix}P=(2−1​−11​) により対角化される。実際 D=(2003)=P−1APD = \begin{pmatrix} 2&0\\0&3\end{pmatrix} = P^{-1} APD=(20​03​)=P−1AP が成立する。よって
etA=PetDP−1=(2−1−11)(e2t00e3t)(1112)=(2e2t−e3t2e2t−2e3t−e2t+e3t−e2t+2e3t)\begin{aligned}
e^{tA} &= P e^{tD} P^{-1}\\
&= \begin{pmatrix} 2&-1\\-1&1\end{pmatrix}
\begin{pmatrix} e^{2t}&0\\0&e^{3t}\end{pmatrix}
\begin{pmatrix}1&1\\1&2\end{pmatrix}\\
&= \begin{pmatrix} 2e^{2t}-e^{3t} & 2e^{2t} - 2e^{3t} \\ -e^{2t} + e^{3t} &-e^{2t} +2e^{3t}\end{pmatrix}
\end{aligned}etA​=PetDP−1=(2−1​−11​)(e2t0​0e3t​)(11​12​)=(2e2t−e3t−e2t+e3t​2e2t−2e3t−e2t+2e3t​)​
となる。初期値 x0=(1−2)\mathbb{x}_0 = \begin{pmatrix} 1\\-2\end{pmatrix}x0​=(1−2​) より
(x1x2)=etAx0=(2e2t−e3t2e2t−2e3t−e2t+e3t−e2t+2e3t)(1−2)=(−2e2t+3e3te2t−3e3t)\begin{aligned}
\begin{pmatrix} x_1\\x_2 \end{pmatrix} &= e^{tA} \mathbb{x}_0\\
&=\begin{pmatrix} 2e^{2t}-e^{3t} & 2e^{2t} - 2e^{3t} \\ -e^{2t} + e^{3t} &-e^{2t} +2e^{3t}\end{pmatrix} \begin{pmatrix} 1\\-2\end{pmatrix}\\
&= \begin{pmatrix} -2e^{2t} + 3e^{3t} \\ e^{2t} -3e^{3t}\end{pmatrix}
\end{aligned}(x1​x2​​)​=etAx0​=(2e2t−e3t−e2t+e3t​2e2t−2e3t−e2t+2e3t​)(1−2​)=(−2e2t+3e3te2t−3e3t​)​
である。よって，x1=−2e2t+3e3t,x2=e2t−3e3tx_1 = -2e^{2t} + 3e^{3t}, x_2=e^{2t}-3e^{3t}x1​=−2e2t+3e3t,x2​=e2t−3e3t
なお，対角化できない場合はジョルダン標準形を用います。ジョルダン標準形の場合は，各ジョルダン細胞について考えれば良いです。ジョルダン細胞は，λiI+N\lambda_i I + Nλi​I+N と分解できたのでした。そのため etλiIetNe^{t \lambda_i I} e^{tN}etλi​IetN を計算することになります。etλiIe^{t\lambda_i I}etλi​I は対角成分が teλite^{\lambda_i}teλi​ の対角行列ですので問題はありません。しかし，etNe^{tN}etN は少々大変です。計算してみると，
etN=I+tN+12t2N2+⋯=(101⋱101)+(0t00t⋱⋱0t00)+12(0t200⋱⋱t2000)+⋯+1(m−1)!(0tm−10⋱000)=(1t12t2⋯1(m−1)!tm−11t⋱⋮⋱⋱12t21t01)\begin{aligned}
e^{tN} &= I + tN + \dfrac{1}{2} t^2 N^2 + \cdots\\
&= \begin{pmatrix} 1&&&&0\\&1&&&\\&&\ddots&&\\&&&1&\\0&&&&1\end{pmatrix} + \begin{pmatrix}0&t&&&0\\&0&t&&\\&&\ddots&\ddots&\\&&&0&t\\0&&&&0\end{pmatrix} + \dfrac{1}{2} \begin{pmatrix} 0&&t^{2}&&0\\&0&&\ddots&\\&&\ddots&&t^{2}\\&&&0&\\0&&&&0\end{pmatrix} +\cdots +\dfrac{1}{(m-1)!} \begin{pmatrix} 0&&&&t^{m-1}\\&0&&&\\&&\ddots&&\\&&&0&\\0&&&&0\end{pmatrix}\\
&= \begin{pmatrix} 1 &t &\dfrac{1}{2} t^2 &\cdots& \dfrac{1}{(m-1)!} t^{m-1}\\ &1&t &\ddots&\vdots \\ &&\ddots&\ddots&\dfrac{1}{2} t^2\\&&&1&t\\0&&&&1\end{pmatrix}
\end{aligned}etN​=I+tN+21​t2N2+⋯=⎝⎛​10​1​⋱​1​01​⎠⎞​+⎝⎛​00​t0​t⋱​⋱0​0t0​⎠⎞​+21​⎝⎛​00​0​t2⋱​⋱0​0t20​⎠⎞​+⋯+(m−1)!1​⎝⎛​00​0​⋱​0​tm−10​⎠⎞​=⎝⎛​10​t1​21​t2t⋱​⋯⋱⋱1​(m−1)!1​tm−1⋮21​t2t1​⎠⎞​​
となります。こうして微分方程式の解が計算できます。
「漸化式の解き方と微分方程式の解き方が似ている」のは意識しておくとよいです。