線形代数II/連立線形微分方程式 の履歴(No.25)
更新- 履歴一覧
- 差分 を表示
- 現在との差分 を表示
- ソース を表示
- 線形代数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)$$
である($\lambda_1\ne 0,\lambda_2\ne 0$ の時)。
$A$ の固有値を求め、対角化しよう。
$$|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}\\$$
を確かめられる。
演習問題†
以下の微分方程式を境界条件 $x_1(0)=1,x_2(0)=2$ の下で解け。
$$ \begin{cases} \frac{d}{dt}x_1(t)=-x_1(t)+x_2(t)\\ \frac{d}{dt}x_2(t)=\phantom{-}x_1(t)-x_2(t)\\ \end{cases} $$
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}$$
を得る(三角関数の形が見えてきた)。 変数分離して両辺を積分すると、
$$\pm\frac{x'}{\sqrt{A^2-x^2}}=\omega$$
$$\pm\int \frac{x'}{\sqrt{A^2-x^2}}dt=\pm\int \frac{dx}{\sqrt{A^2-x^2}}=\pm\cos^{-1}\frac xA=\omega t+\theta\hspace{1cm}\theta:\text{ 積分定数}$$
両辺の $\cos$ を取ると $\cos x=\cos(-x)$ より $\pm$ は消えて、
$$ 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)$ を求める式になる。
$\omega=\sqrt\lambda$ であったことを念頭に置きつつ「行列の関数」をイメージすれば、
$$\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$ を含む「行列の関数」は先の授業でやったような「固有値の整数次のみの多項式」の形では表せないが、多項式が分数次の項を含んでも良いとして定義を拡張すれば、ほぼ同じ議論でこの書き方を正当化できる。
例題†
次の連立微分方程式を考える。
$$\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)$$
と書けるから、$A$ の固有値がすべて正であるとき、その解を形式的に
$$\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}>0,\lambda_2=\frac{3k}{m}>0$ であり、
対応する固有ベクトルはそれぞれ $\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\omega_1}\begin{pmatrix} \sin{\omega_1}t+\tfrac1{\sqrt3}\sin{\omega_2}t& \sin{\omega_1}t-\tfrac1{\sqrt3}\sin{\omega_2}t\\ \sin{\omega_1}t-\tfrac1{\sqrt3}\sin{\omega_2}t& \sin{\omega_1}t+\tfrac1{\sqrt3}\sin{\omega_2}t\\ \end{pmatrix}\dot{\bm x}(0) $$
が $\bm x(0),\dot{\bm x}(0)$ で表わした $\bm x(t)$ の表式となる。 ここで、$\omega_1=\sqrt{\lambda_1}=\sqrt{\frac{k}{m}},\omega_2=\sqrt{\lambda_2}=\sqrt3\omega_1$ である。
物理的意味†
計算途中で表れた $y_1(t),y_2(t)$ は物理的にどんな意味を持つだろうか?
これらはそれぞれ個別に調和振動子の方程式
$$ \begin{cases} \ddot y_1=-\omega_1^2y_1\\ \ddot y_2=-\omega_2^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)$ を
$$ \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)} $$
と書きなおすと分かるように、
それぞれ表す。
逆位相モードでは3つのバネがすべて働くのに対して、
同位相モードでは真ん中のバネが働かないため、
逆位相モードの振動数は同位相モードよりも $\sqrt3$ 倍だけ大きくなっている。
この系の任意の運動を、この2つの振動モードの重ね合わせとして表せる、 というのが $\bm x(t)=P\bm y(t)$ の変換の意味である。
解析力学的な話をするならば、$\bm x$ と $\bm y$ との間は双方向に変数変換可能であるから、 $\bm x$ の代りに $\bm y$ を独立変数と取ることができて、このとき2つの質点の連成振動の系は 2つの独立な調和振動子の系に変換可能であると言える。
これは非常に重要な変換で、例えば固体物理で扱う複雑に振動する格子振動を無数の調和振動子に変数変換した際の、それぞれの振動モードの振幅が量子化したものがフォノンであったり、真空中の電磁場をやはり無数の調和振動子に変数変換した際の、それぞれの振動モードの振幅が量子化したものがフォトン(光子)だったりといったように、応用範囲も広いものである。
演習問題†
$x$ 軸上のみを動ける3つの質点が、ばね定数 $k$ を持つ2つのばねで連結された図のような系を考える。ただし、真ん中の質点の質量を $m$、両端の質点の質量をそれぞれ $m/2$ とする。このとき運動方程式は次のようになる。以下の問いに従って、この方程式の一般解を基準運動に分離せよ。
$$ \begin{cases} (m/2)\ddot x_1(t)=\phantom{-}k\{x_2(t)-x_1(t)\}\\ \hspace{6.2mm}m\ddot x_2(t)=-k\{x_2(t)-x_1(t)\}+k\{x_3(t)-x_2(t)\}\\ (m/2)\ddot x_3(t)=\phantom{-k\{x_2(t)-x_1(t)\}}-k\{x_3(t)-x_2(t)\}\\ \end{cases} $$
(1) $\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\\x_3(t)\end{pmatrix}$ に対して、与えられた方程式を $\ddot \bm x=-A\bm x$ と表した際の係数行列 $A$ を求めよ。
(2) $A$ の固有値 $\lambda_1,\lambda_2,\lambda_3$ を求めよ。また、対応する固有ベクトルを1つずつ求め $\bm p_1,\bm p_2,\bm p_3$ とせよ。ただし、$\lambda_1<\lambda_2<\lambda_3$ とする。
(3) $P=\Big(\bm p_1\ \bm p_2\ \bm p_3\Big)$ に対して $P^{-1}$ を求めよ。
(4) $\bm y=\begin{pmatrix}y_1(t)\\y_2(t)\\y_3(t)\end{pmatrix}$ が $\bm x=P\bm y$ を満たすとき、 $y_1,y_2,y_3$ に対する運動方程式を求め、それぞれの一般解 $y_1^*,y_2^*,y_3^*$ を適当なパラメータを含む形で書き下せ。
(5) $\bm x$ の一般解は $\bm x(t)=\bm p_1y_1^*+\bm p_3y_3^*+\bm p_3y_3^*$ で与えられることを念頭に置き、 $y_1,y_2,y_3$ で表される基準運動がどのような運動を表すか説明せよ。
発展:$n$ 粒子の系†
両端を壁で挟まれた質量 $m$ の2つの質点が3つのバネ(ばね定数 $k$)で繋がれた系を上で解いたが、
これを一般化して 両端を壁で挟まれた $n$ 個の質点が $n+1$ 個のバネで繋がれた系について考える。
解くべき方程式†
$n$ 粒子の問題に対する方程式は2粒子の時と同様、
$\bm x(t)=\begin{pmatrix}x_1(t)\\x_2(t)\\\vdots\\x_n(t)\end{pmatrix}$、 $\displaystyle A=\frac{k}{m}\begin{pmatrix}2&-1\\-1&2&-1\\&\ddots&\ddots&\ddots\\&&-1&2&-1\\&&&-1&2\end{pmatrix}$
に対して、
$$ \ddot {\bm x}(t)=-A\bm x(t) $$
となるが、同じ問題が固体物理学においてフォノンのモデル系として表れるため、ここでしっかり見ておこう。
行列の関数を使った解の形†
$A$ は $n$ 次対称行列であり、$n$ 個の互いに直交する固有ベクトルと、対応する正の値を取る $n$ 個の実数固有値を求めることできるため、それらを使って解を
$$\bm x(t)=P\begin{pmatrix}\cos\sqrt{\lambda_1}t\\&\ddots\\&&\cos\sqrt{\lambda_n}t\end{pmatrix}P^{-1}\bm x(0)+ P\begin{pmatrix}\frac{\sin\sqrt{\lambda_1}t}{\sqrt{\lambda_1}}\\&\ddots\\&&\frac{\sin\sqrt{\lambda_n}t}{\sqrt{\lambda_n}}\end{pmatrix}P^{-1}\dot{\bm x}(0)$$
と書けるのは $n=2$ の時と同様である。
固有値と固有関数、分散関係†
上記の $A$ から正攻法で固有値・固有ベクトルを求めるのはちょっと手間がかかるため、 以下では多少天下りな方法を用いて $\lambda_i$ および対応する固有ベクトルを求めることにする。
上記の方程式は位置が固定された仮想的な質点 $x_0(t)=x_{n+1}(t)=0$ を追加することにより、
$$\ddot x_j=-(k/m)(-x_{j-1}+2x_j-x_{j+1})\hspace{1cm}j=1,2,\dots ,n$$
のようにすべての質点に対して同一な条件に書き換えることができる(固定端の条件 = ディリクレ境界条件)。
そこで
$$ \bm x(t)=A\cos(\omega t+\theta)\begin{pmatrix} \sin\alpha\\ \sin2\alpha\\ \vdots\\ \sin n\alpha\\ \end{pmatrix} $$
すなわち $x_j=A\sin j\alpha\cos(\omega t+\theta)$ と置いてみると、
$$\begin{aligned} x_{j\pm 1}&=A\sin (j\pm1)\alpha\cos(\omega t+\theta)\\ &=\pm \sin\alpha \cdot A\cos j\alpha \cos(\omega t+\theta)+\cos\alpha\cdot\underbrace{A\sin j\alpha \cos(\omega t+\theta)}_{x_j}\\ \end{aligned}$$
より、
$$ x_{j-1}+x_{j+1}=2\cos\alpha\cdot x_j $$
に注意して先の条件式を変形すると
$$ \begin{aligned} \ddot x_j&=-(k/m)(-x_{j-1}+2x_j-x_{j+1})\\ {}-\omega^2 x_j&=-\frac{2k}{m}(1-\cos \alpha)x_j\\ \end{aligned} $$
となり、$\omega$ と $\alpha$ との間に「分散関係」
$$ \omega_\alpha^2=\frac{2k}{m}(1-\cos \alpha)\\ $$
が成り立つとき $\bm x(t)$ は与えられた微分方程式の解となる。そしてこのとき、$\bm x$ は固有値 $\lambda_\alpha=\omega_\alpha^2=\frac{2k}{m}(1-\cos \alpha)$ に対する $A$ の固有関数となる。 (物理学において分散関係とは 波数 = $2\pi/$(波長) と (角)振動数あるいはエネルギーとの関係を表した関係式のことである。 今は $\alpha$ が波数を決めており、これと $\omega$ との関係がここに表されている。)
固定端の条件 $x_0=x_{n+1}=0$ は、
$$ \sin\alpha 0=\sin(n+1)\alpha=0 $$
を表すため、ここからある整数 $N$ を用いて
$$ (n+1)\alpha=N\pi $$
と書ける、という $\alpha$ に対する条件が加わる。すなわち $\alpha$ は離散値
$$ \alpha_N=\frac{N\pi}{n+1} $$
を取ることになる。
$N=0$ は $\bm x(t)=\bm 0$ を表し、また、
$$ \sin j\alpha_{(N+2n+2)}=\sin \Big(\frac{j(N+2n+2)\pi}{n+1}\Big) =\sin \Big(\frac{jN\pi}{n+1}+j2\pi\Big)=\sin j\alpha_N $$
であり、
$$ \begin{aligned} \sin j\alpha_{(N+n+1)}&=\sin \Big(\frac{j(N+n+1)\pi}{n+1}\Big)\\ &=-\sin \Big(-\frac{j(N+n+1)\pi}{n+1}\Big)\\ &=-\sin \Big(2\pi-\frac{j(N+n+1)\pi}{n+1}\Big)\\ &=-\sin \Big(\frac{j(n+1-N)\pi}{n+1}\Big)\\ &=-\sin j\alpha_{(n+1-N)} \end{aligned} $$
であるから、$N>n$ なる $N$ は $0\le N'\le n$ を満たす $N'$ と一次従属になる。
したがって、一次独立な $\bm x$ を与えるのは $N=1,\dots,n$ に限られ、 これで $n$ 本の一次独立な固有ベクトル $\bm x_N(t)$ と固有値
$$\lambda_N=-\omega_N^2=-\frac{2k}{m}(1-\cos \alpha_N)$$
が求まったことになる。
つまり、$n$ 粒子の系には $n$ 個の固有振動モードが存在して、それぞれのモードは振動数に対応して決まる特定の波数を持つのである。この話は1次元格子振動のフォノンのモデルとなっている。
蛇足:固有関数の直交性†
$\bm x_N$ はエルミート行列の異なる固有値に対する固有ベクトルであるから、独立であるだけでなく直交する。実際に計算してみると、
$$ \begin{aligned} (\bm x_N,\bm x_{N'}) &\propto\sum_{j=1}^n\sin j\alpha_N\sin j\alpha_{N'}\\ &=-\frac12\sum_{j=1}^n\big[\cos j(\alpha_N+\alpha_{N'})-\cos j(\alpha_{N}-\alpha_{N'})\big]\\ &=-\frac12\sum_{j=1}^n\Big[\cos \frac{jM\pi}{n+1}-\cos \frac{jM'\pi}{n+1}\Big]\\ \end{aligned} $$
ただし、$M=N+N',M'=N-N'$ であるが、$M\ne 0$ に対して
$$ \sum_{j=1}^n\cos \frac{jM\pi}{n+1}= \begin{cases} 0&M\text{が奇数}\\ {}-1&M\text{が偶数} \end{cases} $$
となることと、$M,M'$ の偶奇が等しいことから、
$$ (\bm x_N,\bm x_{N'})=0 $$
が導かれる。
上記の $\cos$ の和は、
$$ \begin{aligned} \sum_{j=1}^n\exp \frac{\pm ijM\pi}{n+1} &=\exp \Big(\frac{\pm iM\pi}{n+1}\Big) \frac{\exp \Big(\frac{\pm inM\pi}{n+1}\Big)-1} {\exp \Big(\frac{\pm iM\pi}{n+1}\Big)-1}\\ &=\exp \Big(\frac{\pm iM\pi}{n+1}\Big) \frac{\exp \Big(\pm iM\pi\pm \frac{-iM\pi}{n+1}\Big)-1} {\exp \Big(\frac{\pm iM\pi}{n+1}\Big)-1}\\ &=\exp \Big(\frac{\pm iM\pi}{n+1}\Big) \frac{(-1)^M\exp \Big(\pm \frac{-iM\pi}{n+1}\Big)-1} {\exp \Big(\frac{\pm iM\pi}{n+1}\Big)-1}\\ &= \frac{(-1)^M-\exp \Big(\frac{\pm iM\pi}{n+1}\Big)} {\exp \Big(\frac{\pm iM\pi}{n+1}\Big)-1}\\ &=\begin{cases} \mp\cot\frac{M\pi}{n+1} &M\text{ が奇数}\\ {}-1&M\text{ が偶数}\\ \end{cases} \end{aligned} $$
から得られる。
微分方程式と固有値問題†
上の問題をもう一度、違う角度から見てみよう。
$$ \ddot {\bm x}(t)=-K\bm x(t) $$
$x_n(t)$ は $n$ と $t$ との関数 $x_n(t)=x(n,t)$ とみなせることに注意して、この方程式を、
$$ \bm x(t)=\underbrace{\bm p\rule[-2.4pt]{0pt}{1pt}}_{n\,\text{依存性}}\cdot\underbrace{ y(t)}_{t\,\text{依存性}} $$
の「変数分離形」を仮定して解くことにする。代入すると、
$$ \bm p\cdot\ddot y(t)=-K\bm p\cdot y(t) $$
$$ \underbrace{\frac{\ddot y(t)}{ y(t)}}_{\text{定数}}\cdot\bm p=-K\bm p $$
右辺は $t$ を含まないから、左辺も $t$ に対して定数となっていなければならない。そこで $\ddot y/ y=-\lambda$ と置くと、
$$ \begin{cases} \frac{d^2}{dt^2} y(t)=-\lambda y(t)\\ K\bm p=\lambda\bm p \end{cases} $$
の2つの固有値問題に分離される。($d^2/dt^2$ は線形変換であることに注目)
$K$ の固有値 $\lambda_i$ と対応する固有ベクトル $\bm p_i$ を求め、それぞれの $\lambda_i$ に対して$\ddot y_i(t)=-\lambda_i y_i(t)$ を解いて $y_i(t)$ を求めると、
$$ \bm x_i(t)=\bm p_i y_i(t) $$
が元の方程式の変数分離形の解となる。
$K$ が対角化可能な場合、任意の解 $\bm x(t)$ を $\bm x_i(t)$ で展開できて、
$$ \bm x(t)=\sum_i c_i\bm x_i(t)=\sum_i c_i\bm p_i y_i(t) $$
これを、
$$ P=\Bigg( \bm p_1\ \ \bm p_2\ \ \dots\ \ \bm p_n \Bigg),\ \ \bm y(t)=\begin{pmatrix} c_1y_1(t)\\ c_2y_2(t)\\ \vdots\\ c_ny_n(t)\\ \end{pmatrix} $$
と置いて書き直すと、
$$ \bm x(t)=P\bm y(t) $$
これがこの系の任意の運動を基準運動(基準振動)へ「分解」する、 あるいは基準運動へ「変数変換」することの意味となる。
特に、「基準運動」というのが粒子の番号 $n$ によらず、 すべての粒子が同じ振動数(時間依存性)をもって運動する運動モードであることが、 この変数分離の解き方からよく理解できたと期待する。
おまけ:グリーン関数法†
この話は「対角化」や「行列の関数」とは関係ないのであるが、微分方程式繋がりということでここで紹介しておく。 カリキュラム上、まだ、「デルタ関数」を学んでいない場合には、 デルタ関数やグリーン関数について学んだ後で戻ってきてざっと目を通せば、 グリーン関数法と線形代数との繋がりを理解できるはずである。
微分方程式:
$$ Tf(x)=g(x) $$
が、「$x$ に対する微分を含むような線形変換 $T$ を未知の関数 $f(x)$ に 適用すると既知の関数 $g(x)$ になる」を表すとき、これを解いて $f(x)$ を求めることを考える。
$T$ の線形性から、恒等的にゼロとなるような $g(x)=0$ を与えれば $f(x)=0$ が解となってしまうため($T$ が正則でなければこれ以外に解はあるかもしれないが)、このような $g(x)$ は $f(x)$ の
このような問題は $T$ が正則であり、
$$TT^{-1}=I$$
を満たすような $T^{-1}$ を求めることさえできれば両辺に $T^{-1}$ をかけることで
$$ f(x)=T^{-1}g(x) $$
として簡単に解けるのであるが、「グリーン関数法」では「グリーン関数」と呼ばれる関数を使って $T^{-1}$ を求め、方程式を解く。
$T$ のグリーン関数とは、$T$ を適用すると ディラックのデルタ関数 を生じるような関数すなわち、
$$ TG(x,x')=\delta(x-x') $$
を満たすような関数 $G(x,x')$ のことである。 (細かいことを言うと この微分方程式に加えて、 元の問題で与えられる境界条件も満たすものでないとならない)
この $G(x,x')$ を使うと、
$$T^{-1}:\rho(x)\mapsto\int G(x,x')\rho(x')dx'$$
と書けるため、
$$ f(x)=\int G(x,x')g(x')dx' $$
が微分方程式の解になる、というのがグリーン関数法のキモである。
実際にこれが $TT^{-1}=I$ を満たすことを、任意の関数 $f(x)$ に対して、
$$ \begin{aligned} TT^{-1}f(x)&=T\int G(x,x')f(x')dx'\\ &=\int TG(x,x')f(x')dx'\\ &=\int \delta(x-x')f(x')dx'\\ &=f(x)\\ \end{aligned} $$
となることから確かめられる。
ここで、
$$ I:f(x)\mapsto\int\delta(x-x')f(x')dx' $$
と書かれたことと、単位行列の要素がクロネッカーのデルタ $\delta_{ij}$ を用いて $I=\big(\delta_{ij}\big)$ と書け、$\bm y=I\bm x$ の要素が
$$ y_{i}=\sum_{j=1}^n\delta_{ij}x_j $$
と書かれたこととが対応している点も良く理解しておくように。
例†
電荷密度分布 $\rho(\bm x)$ の作る電位 $\phi(\bm x)$ を求めたい。
基本方程式はポアソン方程式
$$-\nabla^2\phi(\bm x)=\frac1\epsilon_0\rho(\bm x) $$
で与えられるが、この方程式は
$$ \underbrace{-\epsilon_0\nabla^2}_T\phi(\bm x)=\rho(\bm x) $$
としたときの線形演算子 $T$ の逆 $T^{-1}$ を求めれば $\phi(x)=T^{-1}\rho(x)$ として解ける。
$\bm x'$ に置かれた単位点電荷の作る電位を無限遠で $\phi(\bm x)=0$ となる境界条件の下で解けば、
$$ \phi_1(\bm x;\bm x')=-\frac{1}{4\pi\epsilon_0}\frac{1}{|\bm x-\bm x'|} $$
となるが、これは
$$ T\phi_1(\bm x;\bm x')=\delta^2(\bm x-\bm x') $$
を表しており、$T$ に対する同境界条件を満たすグリーン関数が $G(\bm x,\bm x')=\phi_1(\bm x;\bm x')$ であることを表している。
したがって、任意の電荷分布 $\rho(\bm x)$ の作る電位を
$$ \begin{aligned} \phi(\bm x)&=T^{-1}\rho(\bm x)=\iiint G(\bm x,\bm x')\rho(\bm x')d^3\bm x'\\ &=-\frac{1}{4\pi\epsilon_0}\iiint \frac{\rho(\bm x')}{|\bm x-\bm x'|}d^3\bm x'\\ \end{aligned} $$
のようにして求められるのである。
この式の物理的な意味は、各点における電荷密度が作る電位を重ね合わせることで全体としての電位が求まる、と読めるが、重ね合わせることが可能な理由は $T$ の線形性にあることを理解せよ。
非線形な方程式ではこのような「重ね合わせ」は通用しない。
質問・コメント†
無題†
()
むずかっしー
*1 参考:http://w3e.kanazawa-it.ac.jp/math/physics/category/mechanics/masspoint_mechanics/simple_harmonic_motion/henkan-tex.cgi?target=/math/physics/category/mechanics/masspoint_mechanics/simple_harmonic_motion/shm_solution.html#solution2