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

更新


  • 追加された行はこの色です。
  • 削除された行はこの色です。
[[前の単元 <<<>線形代数II/行列の関数]]
               [[線形代数II]]
               [[>>> 全体のまとめ>線形代数II/まとめ]]

#contents

* 概要 [#rd73f160]

行列の対角化や行列の関数を利用して、連立線形微分方程式を解く。
&katex();
* 1階線形微分方程式 [#f6af6a12]

** 基本 [#wd466ab6]

通常の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)$$

と表すことができる。

** 連立方程式 [#j6598f91]

未知の関数 $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}$ を固有値を用いてどのように書いたらよいるかを理解するには 
[[発展:ジョルダン標準形>線形代数I/広義固有空間の構造とジョルダン標準形]] 
について学ぶ必要がある。

** 例題 [#c8f38900]

$\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$ の下で解け。

#hr();

$\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}$$

が求める解である。

*** 検算 [#e851f146]

上記の $\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階線形微分方程式 [#b5cf68c1]

** 基本 [#a9451aba]

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

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

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

これを解くには((参考: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))、
まず全体に $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$$

と表せる。

** 連立方程式 [#c1b36c65]

未知の関数 $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$ を含む「行列の関数」は先の授業でやったような「固有値の整数次のみの多項式」の形では表せないが、多項式が分数次の項を含んでも良いとして定義を拡張すれば、ほぼ同じ議論でこの書き方を正当化できる。

** 例題 [#ha0f47df]

#ref(coupled_oscillation1.png,right,around,50%);

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

$$\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)
$$

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


*** 物理的意味 [#q4312501]

計算途中で表れた $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)} 
$$

と書きなおすと分かるように、

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

それぞれ表す。

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

この系の任意の運動を、この2つの振動モードの重ね合わせとして表せる、
というのが $\bm x(t)=P\bm y(t)$ の変換の意味である。

解析力学的な話をするならば、$\bm x$ と $\bm y$ との間は双方向に変数変換可能であるから、
$\bm x$ の代りに $\bm y$ を独立変数と取ることができて、このとき2つの質点の連成振動の系は
2つの独立な調和振動子の系に変換可能であると言える。

** 発展:$n$ 粒子の系 [#cbfb66a5]

*** 解くべき方程式 [#waef2077]

$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)
$$

となるが、同じ問題が固体物理学においてフォノンのモデル系として表れるためここで見ておこう。

*** 行列の関数を使った解の形 [#u91dbff9]

$A$ は $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$ の時と同様である。

*** 固有値と固有関数、分散関係 [#a2528248]

上記の $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$$

のようにすべての質点に対して同一な条件に書き換えることができる(固定端の条件 = ディリクレ境界条件)。

ここで $x_j=\sin j\alpha\cos\omega t$ すなわち、

$$
\bm x(t)=\cos\omega t\begin{pmatrix}
\sin\alpha\\
\sin2\alpha\\
\vdots\\
\sin n\alpha\\
\end{pmatrix}
$$

と置いてみると、

$$\begin{aligned}
x_{j\pm 1}&=\sin (j\pm1)\alpha\cos\omega t\\
&=\pm\sin\alpha\cos j\alpha \cos\omega t+\cos\alpha\,\underbrace{\sin j\alpha \cos\omega t}_{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}
$$

となる。これは、

$$
\begin{aligned}
&\ddot {\bm x}_\alpha(t)=-\omega^2\bm x_\alpha(t)=\\
&-A\bm x_\alpha(t)=-\frac{2k}{m}(1-\cos\alpha)\bm x_\alpha(t)
\end{aligned}
$$

を表しており、ここから $\omega$ と $\alpha$ との間の「分散関係」

$$
\lambda_\alpha=\omega_\alpha^2=\frac{2k}{m}(1-\cos \alpha)\\
$$

が得られる。(物理学において分散関係とは 波数 = $2\pi/$(波長) と
(角)振動数あるいはエネルギーとの関係を表した関係式のことである。
今は $\alpha$ が波数を決めており、これと $\omega$ との関係がここに表されている。)

境界条件 $x_0=x_{j+1}=0$ は、

$$
\sin\alpha 0=\sin(n+1)\alpha=0
$$

を表すため、ここから $\alpha$ に対する条件が任意の整数 $N$ を用いて

$$
(n+1)\alpha=N\pi
$$

と書ける。すなわち $\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)$$

が求まったことになる。

*** 固有関数の直交性 [#leff2001]

$\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}
$$

から得られる。

* おまけ:グリーン関数法 [#g65dc6f9]

この話は対角化や行列の関数とは関係ないのであるが、微分方程式繋がりということでここで紹介しておく。
カリキュラム上、またデルタ関数などを学んでいないかもしれないのでここでは深入りしないが、
カリキュラム上、まだ、「デルタ関数」を学んでいないかもしれないのでここでは深入りしないが、
デルタ関数やグリーン関数について学んだ後で ざっと目を通してもらえると、
線形代数との繋がりを理解できるはずである。
グリーン関数法と線形代数との繋がりを理解できるはずである。

未知の関数 $f(x)$ に $x$ での微分を含むような線形変換 $T$ を適用すると既知の関数 $g(x)$ となるという、

$$
Tf(x)=g(x)
$$

のような微分方程式を解いて $f(x)$ を求めることを考える。

$T$ の線形性から、恒等的にゼロとなるような $g(x)=0$ を与えれば解は $f(x)=0$ となってしまうため、このような $g(x)$ は $f(x)$ の&ruby(みなもと){源};であるという意味でソース関数と呼ばれる。

このような問題は、

$$TT^{-1}=I$$

を満たすような $T^{-1}$ を求めることさえできれば両辺に $T^{-1}$ をかけることで

$$
f(x)=T^{-1}g(x)
$$

として簡単に解けるのであるが、
実際にこのような解き方を行う手法を「グリーン関数法」と呼ぶ。
として簡単に解けるのであるが、「グリーン関数法」とは「グリーン関数」と呼ばれる関数を使うことで $T^{-1}$ を求め、方程式を解く手法である。

というのも、このような $T^{-1}$ を以下のようにして求められるためだ。

$T$ のグリーン関数とは、$T$ を適用すると [[ディラックのデルタ関数>Wikipedia:ディラックのデルタ関数]] を生じるような関数すなわち、

$$
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')\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
$$

と書かれたこととが対応している点も良く理解しておくように。

** 例 [#r6f1a545]

電荷密度分布 $\rho(\bm x)$ の作る電位 $\phi(\bm x)$ を求めたい。

基本方程式はポアソン方程式

$$-\nabla^2\phi(\bm x)=\frac1\epsilon_0\rho(\bm x)
$$

で与えられるが、この方程式は

$$
\underbrace{-\epsilon_0\nabla^2}_T\phi(x)=\rho(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$ の線形性にあることを理解せよ。

非線形な方程式ではこのような「重ね合わせ」は通用しない。

[[前の単元 <<<>線形代数II/行列の関数]]
               [[線形代数II]]
               [[>>> 全体のまとめ>線形代数II/まとめ]]

* 質問・コメント [#j4c326c0]

#article_kcaptcha
**無題 [#abcb204a]
> (&timetag(2018-08-23T06:02:03+09:00, 2018-08-23 (木) 15:02:03);)~
~
むずかっしー~

//
- 勉強になります --  &new{2019-06-08 (土) 18:11:21};
- 特性根が重解 と 複素数解の場合も言ってほしいけど --  &new{2019-07-15 (月) 00:39:43};
- 参考になります -- [[雑魚]] &new{2020-07-26 (日) 23:25:58};

#comment_kcaptcha


Counter: 235158 (from 2010/06/03), today: 13, yesterday: 0