線形代数II/連立線形微分方程式 のバックアップ(No.14)

更新


前の単元 <<<                線形代数II                >>> 全体のまとめ

概要

行列の対角化を利用して、連立線形微分方程式を解く。

1階連立線形微分方程式

基本

通常の1階線形微分方程式

$$\frac{dx(t)}{dt}=ax(t)$$

を解くには、まず両辺を $x(t)$ で割り、

$$\frac{\ \frac{dx(t)}{dt}\ }{x(t)}=\frac{d}{dt}\log x(t)=a$$

中辺、右辺を積分して($C'$: 積分定数)、指数関数を取ることにより

$$\log x(t)=at + C'$$

$$x(t)=Ce^{at}$$

とすればよい。ここで $C=e^{C'}$ であるが、$t=0$ を代入すれば

$$x(0)=C$$

となることに注意すれば、

$$x(t)=e^{at}\,x(0)$$

と表すことができる。

連立方程式

未知の関数 $x_1(t),x_2(t)$ に対して 1階の線形微分方程式が連立方程式として以下のように与えられているとする。

$$\left\{ \begin{matrix} \displaystyle\frac{dx_1(t)}{dt}=\dot x_1(t)=ax_1(t)+bx_2(t) \\ \displaystyle\frac{dx_2(t)}{dt}=\dot x_2(t)=cx_1(t)+dx_2(t) \end{matrix} \right.$$

この方程式は

$$\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\end{pmatrix}, \ \ A=\begin{pmatrix}a&b\\c&d\end{pmatrix}$$

を用いれば

$$\frac{d\bm x(t)}{dt}=\dot{\bm x}(t)=A\bm x(t)$$

と表せるから、行列の関数を用いることで方程式の解を

$$ \bm x(t)=e^{At}\bm x(0) $$

として求められる!!!(ドヤッ

で終わってしまっても良いのだけれど念のため、$A$ が対角化可能な場合について行列の関数を用いない方法との対比を確認しておこう。

$A$ を $P^{-1}AP=D=\begin{pmatrix}\lambda_1&\\&\lambda_2\end{pmatrix}$ の形で対角化できるなら、1年生の二次形式のところでやったのと同様に、

$$ \bm x(t)=P\bm y(t)=P\begin{pmatrix}y_1(t)\\y_2(t)\end{pmatrix} $$

を満たす $\bm y(t)$ へと変数変換することで、

$$ \dot{\bm x}(t)=\frac{d}{dt}\{P\bm y(t)\}=P\dot{\bm y}(t)=A\bm x(t)=AP\bm y(t) $$

すなわち、

$$\dot{\bm y}(t)=P^{-1}AP\bm y(t)=D\bm y(t)$$

のように非対角項がゼロとなる方程式が得られる。これは

$$\left\{\begin{matrix} \dot y_1(t)=\lambda_1y_1(t)\\ \dot y_2(t)=\lambda_2y_2(t) \end{matrix} \right.$$

を表わすから、それぞれ $y_1(t),y_2(t)$ に対して独立に解くことができて、

$$\left\{\begin{matrix} y_1(t)=y_1(0)e^{\lambda_1t}\\ y_2(t)=y_2(0)e^{\lambda_2t}\end{matrix} \right.$$

あるいは $$\bm y(t)=\begin{pmatrix}e^{\lambda_1t}\\&e^{\lambda_2t}\end{pmatrix}\bm y(0)$$

のように解が得られる。

このように、方程式を単純化する基底を取って解を表わした座標を「基準座標」と呼ぶ。

最終的な問題は基準座標に対する解 $\bm y(t)$ を得ることではなく、 $\bm x(t)$ を求めることなので、 $\bm y(0)=P^{-1}\bm x(0)$ を使って変形すると、

$$\bm x(t)=P\bm y(t) =P\begin{pmatrix}e^{\lambda_1t}\\&e^{\lambda_2t}\end{pmatrix}\bm y(0) =P\begin{pmatrix}e^{\lambda_1t}\\&e^{\lambda_2t}\end{pmatrix}P^{-1}\bm x(0) $$

が得られ、これが初期条件 $\bm x(0)$ が与えられた際に $\bm x(t)$ を求めるための式となる。

先にやったように $e^{At}=P\begin{pmatrix}e^{\lambda_1t}\\&e^{\lambda_2t}\end{pmatrix}P^{-1}$ であることを用いてこの式を、

$$\bm x(t)=e^{At}\bm x(0)$$

と書き直すと、通常の線形微分方程式の解とほとんど同じに解ける、というわけだ。

$$ \begin{array}{lll} \dot x(t)=ax(t) &\to& x(t)=x(0)e^{at} \\ \dot{\bm x}(t)=A\bm x(t) &\to& \bm x(t)=e^{At}\,\bm x(0) \\ \end{array} $$

またこの $\bm x(t)$ が与えられた微分方程式を満たすことも容易に確認できる。

$$\dot{\bm x}(t)=\frac{d}{dt}\left\{e^{At}\bm x(0)\right\} =\left\{\frac{d}{dt}e^{At}\right\}\bm x(0) =Ae^{At}\bm x(0) =A\bm x(t) $$

$A$ が対角化可能でない場合にも $\bm x(t)=e^{At}\bm x(0)$ と書いて問題ないが、 $e^{At}$ を固有値を用いてどのように書いたらよいるかを理解するには 発展:ジョルダン標準形 について学ぶ必要がある。

例題

$\left\{\begin{matrix} \frac{dx_1(t)}{dt}&=3x_1(t)+\phantom 3x_2(t)\\[3mm] \frac{dx_2(t)}{dt}&=\phantom 3x_1(t)+3x_2(t) \end{matrix}\right.$  を初期条件 $x_1(0)=6,x_2(0)=2$ の下で解け。


$\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\end{pmatrix}$、 $A=\begin{pmatrix}3&1\\1&3\end{pmatrix}$ と置けば、 $\dot{\bm x}(t)=A\bm x(t)$ と書けるから、

$\bm x(t)=e^{tA}\bm x(0)$ と表せる。

ただし、$P^{-1}AP=\begin{pmatrix}\lambda_1\\&\lambda_2\end{pmatrix}$ と対角化できるならさらに変形して、

$$\bm x(t)=P\begin{pmatrix}e^{\lambda_1t}\\&e^{\lambda_2t}\end{pmatrix}P^{-1}\bm x(0)$$

である。今の場合、

$$|A-\lambda E|=\begin{vmatrix}3-\lambda&1\\1&3-\lambda\end{vmatrix}=(3-\lambda)^2-1=(4-\lambda)(2-\lambda)$$

より、$\lambda_1=4,\lambda_2=2$

それぞれに対応する固有ベクトルは、

$(A-4 E)\bm x=\begin{pmatrix}-1&1\\1&-1\end{pmatrix}\bm x=\bm 0$ より $\bm x=s\begin{pmatrix}1\\1\end{pmatrix}$

$(A-2 E)\bm x=\begin{pmatrix}1&1\\1&1\end{pmatrix}\bm x=\bm 0$ より $\bm x=s\begin{pmatrix}1\\-1\end{pmatrix}$

したがって、$P=\begin{pmatrix}1&1\\1&-1\end{pmatrix}$、$P^{-1}=\frac{1}{2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}$、$P^{-1}AP=\begin{pmatrix}4\\&2\end{pmatrix}$ となって、

$$\begin{aligned} \bm x(t)&=\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}e^{4t}\\&e^{2t}\end{pmatrix}\frac{1}{2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}6\\2\end{pmatrix}\\ &=\begin{pmatrix}e^{4t}&e^{2t}\\e^{4t}&-e^{2t}\end{pmatrix}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}3\\1\end{pmatrix}\\ &=\begin{pmatrix}e^{4t}&e^{2t}\\e^{4t}&-e^{2t}\end{pmatrix}\begin{pmatrix}4\\2\end{pmatrix}\\ &=\begin{pmatrix}4e^{4t}+2e^{2t}\\4e^{4t}-2e^{2t}\end{pmatrix}\\ \end{aligned}$$

が求める解である。

検算

上記の $\bm x(t)$ より

$$\bm x(0)=\begin{pmatrix}6\\2\end{pmatrix}$$

および、

$$\dot{\bm x}(t)=\begin{pmatrix}16e^{4t}+4e^{2t}\\16e^{4t}-4e^{2t}\end{pmatrix}\\$$

$$A\bm x(t)=\begin{pmatrix}3&1\\1&3\end{pmatrix}\begin{pmatrix}4e^{4t}+2e^{2t}\\4e^{4t}-2e^{2t}\end{pmatrix}=\begin{pmatrix}16e^{4t}+4e^{2t}\\16e^{4t}-4e^{2t}\end{pmatrix}\\$$

を確かめられる。

2階連立線形微分方程式

基本

通常の2階線形微分方程式が次のように与えられたとする。

$$\frac{d^2x(t)}{dt^2}=-\omega^2x(t)\hspace{10mm}(\omega>0)$$

$$x''=-\omega^2x$$

これを解くには*1参考:http://w3e.kanazawa-it.ac.jp/m...、 まず全体に $x'$ をかけてから積分して、

$$x'x''+\omega^2xx'=0$$

$$\frac12x'^2+\frac12\omega^2x^2=C_0\hspace{1cm}C_0:\text{ 積分定数}$$

を得る(この関係は物理学ではエネルギー保存則と関係が深い)。 $x$ が実数関数であれば明らかに $C_0\ge 0$ であるから $2C_0=\omega^2A^2\ge 0$ と書くことにすると、

$$x'^2=\omega^2A^2-\omega^2x^2$$

すなわち、

$$x'=\pm\omega\sqrt{A^2-x^2}$$

を得る(三角関数の形が見えてきた)。 変数分離して両辺を積分すると、

$$\mp\frac{x'}{\sqrt{A^2-x^2}}=\omega$$

$$\mp\int \frac{x'}{\sqrt{A^2-x^2}}dt=\mp\int \frac{dx}{\sqrt{A^2-x^2}}=\mp\cos^{-1}\frac xA=\omega t+\theta\hspace{1cm}\theta:\text{ 積分定数}$$

両辺の $\cos$ を取ると $\cos x=\cos(-x)$ より $\mp$ は消えて、

$$ x(t)=A\cos(\omega t+\theta)=C\cos\omega t+C'\sin\omega t $$

ただし、$C=A\cos\theta,C'=-A\sin\theta$ と表せる。

$$\dot x(t)=-\omega C\sin\omega t+\omega C'\cos\omega t$$

より、

$$x(0)=C, \dot x(0)=\omega C'$$

であるから $t=0$ における境界条件を用いれば

$$x(t)=x(0)\cos\omega t+\frac{\dot x(0)}{\omega}\sin\omega t$$

と表せる。

連立方程式

未知の関数 $x_1(t),x_2(t)$ が次の2階線形連立微分方程式を満たす時、

$$\left\{ \begin{matrix} \displaystyle\frac{d^2x_1(t)}{dt^2}=\ddot x_1(t)=-ax_1(t)-bx_2(t) \\ \displaystyle\frac{d^2x_2(t)}{dt^2}=\ddot x_2(t)=-cx_1(t)-dx_2(t) \end{matrix} \right.$$

この方程式はベクトル関数 $\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\end{pmatrix}$、 行列 $A=\begin{pmatrix}a&b\\c&d\end{pmatrix}$ を用いて、

$$\frac{d^2\bm x(t)}{dt^2}=\ddot{\bm x}(t)=-A\bm x(t)$$

と表せる。

もし $A$ を $P^{-1}AP=D=\begin{pmatrix}\lambda_1&\\&\lambda_2\end{pmatrix}$ の形で対角化できるなら、

$$ \bm x(t)=P\bm y(t)=P\begin{pmatrix}y_1(t)\\y_2(t)\end{pmatrix} $$

の変数変換により、

$$\ddot{\bm y}(t)=-P^{-1}AP\bm y(t)=-D\bm y(t)$$

が得られる。

$\lambda_1>0,\ \lambda_2>0$ のとき、$\omega_1=\sqrt{\lambda_1},\ \omega_2=\sqrt{\lambda_2}$ と置けば、これは

$$\left\{\begin{matrix} \ddot y_1(t)=-\lambda_1y_1(t)=-\omega_1^2y_1(t)\\ \ddot y_2(t)=-\lambda_2y_2(t)=-\omega_2^2y_2(t) \end{matrix} \right.$$

という方程式を表わすから、その解を

$$\left\{\begin{matrix} y_1(t)=y_1(0)\cos\omega_1t+\frac{\dot y_1(0)}{\omega_1}\sin\omega_1t\\ y_2(t)=y_2(0)\cos\omega_2t+\frac{\dot y_2(0)}{\omega_2}\sin\omega_2t \end{matrix}\right.$$

あるいは

$$\bm y(t)=\begin{pmatrix}\cos\omega_1t\\&\cos\omega_2t\end{pmatrix}\bm y(0)+ \begin{pmatrix}\frac{\sin\omega_1t}{\omega_1}\\&\frac{\sin\omega_1t}{\omega_1}\end{pmatrix}\dot{\bm y}(0)$$

と表せる。$\bm y(0)=P^{-1}\bm x(0), \dot{\bm y}(0)=P^{-1}\dot{\bm x}(0)$ を用いて書き直すと、

$$\bm x(t)=P\bm y(t)=P\begin{pmatrix}\cos\omega_1t\\&\cos\omega_2t\end{pmatrix}P^{-1}\bm x(0)+ P\begin{pmatrix}\frac{\sin\omega_1t}{\omega_1}\\&\frac{\sin\omega_1t}{\omega_1}\end{pmatrix}P^{-1}\dot{\bm x}(0)$$

これが、初期条件 $\bm x(0),\dot{\bm x}(0)$ が与えられた際に $\bm x(t)$ を求める式になる。

行列の関数をイメージすれば、

$$\bm x(t)=\Bigg[\cos\left(\sqrt{A}\,t\right)\Bigg]\bm x(0)+\Bigg[\frac{\sin\left(\sqrt{A}\,t\right)}{\sqrt{A}}\Bigg]\dot{\bm x}(0)$$

$$x(t)=x(0)\cos\omega t+\frac{\dot x(0)}{\omega}\sin\omega t$$

の形であり、通常の調和振動の式との類似性が際立つ。
(ここで表れた $\sqrt A$ を含む「行列の関数」はこの授業でやった整数次のみの多項式の形では表せないが、多項式が分数次の項を含んでも良いとして定義を拡張すればほぼ同じ議論でこの書き方を正当化できる)

例題

coupled_oscillation1.png

次の連立微分方程式を考える。

$$\left\{\begin{matrix} m\ddot x_1(t)=-kx_1(t)-k\{x_1(t)-x_2(t)\}\\ m\ddot x_2(t)=-kx_2(t)-k\{x_2(t)-x_1(t)\}\\ \end{matrix}\right.$$

これは右図のように、3つのバネで繋がれた2つの質点の運動を記述する方程式を考えているのであるが、 物理モデルから上記方程式を組み立てる部分は数学の範疇ではないので、線形代数学においては 上記方程式をどのように解くかに焦点を絞る。

$\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\end{pmatrix}$、 $\displaystyle A=\frac{k}{m}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}$ と置けば、上式は

$$\ddot{\bm x}(t)=-A\bm x(t)$$

と書けるから、その解を形式的に

$$\bm x(t)=\Bigg[\cos\sqrt{A}t\Bigg]\bm x(0)+\Bigg[\frac{\sin\sqrt{A}t}{\sqrt A}\Bigg]\dot{\bm x}(t)$$

と書ける。

今の場合、$|(m/k)A-\lambda' E|=\begin{vmatrix}2-\lambda'&-1\\-1&2-\lambda'\end{vmatrix}=(2-\lambda')^2-1=(1-\lambda')(3-\lambda')=0$ より、

$\lambda'=1,3$ すなわち $\lambda_1=\frac{k}{m},\lambda_2=\frac{3k}{m}$ であり、

対応する固有ベクトルはそれぞれ $\begin{pmatrix}1\\1\end{pmatrix},\begin{pmatrix}1\\-1\end{pmatrix}$ であるから、

$$P=\begin{pmatrix}1&1\\1&-1\end{pmatrix}$$ と置けば、

$$P^{-1}AP=\begin{pmatrix}\lambda_1&0\\0&\lambda_2\end{pmatrix}$$

として対角化できる。

このとき、

$$\bm x(t)= P\begin{pmatrix}\cos\sqrt{\lambda_1}t&0\\0&\cos\sqrt{\lambda_2}t\end{pmatrix}P^{-1}\bm x(0)+ P\begin{pmatrix}\frac{\sin\sqrt{\lambda_1}t}{\sqrt{\lambda_1}}&0\\0&\frac{\sin\sqrt{\lambda_2}t}{\sqrt{\lambda_2}}\end{pmatrix}P^{-1}\dot{\bm x}(0)$$

と表せ、

$$ \begin{aligned} P\begin{pmatrix}a\\&b\end{pmatrix}P^{-1}&=\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}a\\ &b\end{pmatrix}\frac{1}{2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\\ &=\frac{1}{2}\begin{pmatrix}a&b\\a&-b\end{pmatrix}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\\ &=\frac{1}{2}\begin{pmatrix}a+b&a-b\\a-b&a+b\end{pmatrix} \end{aligned} $$

を用いて計算を進めると、

$$ \bm x(t)= \frac{1}{2}\begin{pmatrix} \cos{\omega_1}t+\cos{\omega_2}t&\cos{\omega_1}t-\cos{\omega_2}t\\ \cos{\omega_1}t-\cos{\omega_2}t&\cos{\omega_1}t+\cos{\omega_2}t \end{pmatrix}\bm x(0)+ \frac{1}{2}\begin{pmatrix} \sin{\omega_1}t+\sin{\omega_2}t&\sin{\omega_1}t-\sin{\omega_2}t\\ \sin{\omega_1}t-\sin{\omega_2}t&\sin{\omega_1}t+\sin{\omega_2}t \end{pmatrix}\dot{\bm x}(0) $$

ただし $\omega_1=\sqrt{\frac{k}{m}},\omega_2=\sqrt{\frac{3k}{m}}$ が $\bm x(0),\dot{\bm x}(0)$ で表わした $\bm x(t)$ の表式となる。

物理的意味

以下、$y_1(t),y_2(t)$ ただし、

$$ \begin{cases} \ddot y_1=-\omega_1^2y_1\\ \ddot y_2=-\omega_1^2y_2\\ \end{cases} $$

$\bm x(t)=P\bm y(t)=\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}y_1(t)\\y_2(t)\end{pmatrix}$ の物理的意味について考察する。

$$ \bm x(t)=\begin{pmatrix}1\\1\end{pmatrix}\underbrace{A_1\cos(\omega_1 t+\theta_1)}_{y_1(t)} +\begin{pmatrix}1\\-1\end{pmatrix}\underbrace{A_2\cos(\omega_2 t+\theta_2)}_{y_2(t)} $$

と書けることから明らかなように、

  • $y_1(t)$ は2つの質点が同位相で振動する振動モードの揺れを、in-phase.gif
  • $y_1(t)$ は2つの質点が逆位相で振動する振動モードの揺れを、opposite-phase.gif

それぞれ表す。

逆位相モードでは3つのバネがすべて働くのに対して、
同位相モードでは真ん中のバネが働かないため、
逆位相モードの振動数は同位相モードよりも $\sqrt3$ 倍だけ大きくなっている。

前の単元 <<<                線形代数II                >>> 全体のまとめ

質問・コメント




無題

()

むずかっしー

  • 勉強になります --
  • 特性根が重解 と 複素数解の場合も言ってほしいけど --
  • 参考になります -- 雑魚?

Counter: 250103 (from 2010/06/03), today: 2, yesterday: 0