連続体力学編 第5章

斉次方程式の解

斉次常微分方程式の初期値問題は、時間発展行列 𝐺𝑡 を求めれば完全に解ける(式(4))。係数行列が定数であれば、 𝐺𝑡 は、ジョルダン分解を用いた解の公式(19)から求められる。

斉次方程式

斉次(常微分)方程式とは、ベクトル 𝒙𝑡 に対する、以下の形の(あるいはこの形に変形できる)微分方程式である: ˙𝒙𝑡=𝐴𝑡⁢𝒙𝑡(1) 添え字の 𝑡 は、 𝑡 の関数であることを表す。 𝐴𝑡 を係数行列という。「時間微分」と「係数行列を掛けること」が等価になるわけである。(「同次」方程式ともいうが、物理だと「同時」と紛らわしい気がする。斉次はせいじと読む。)

代表的な例として、第4章で扱った、バネの運動方程式がある:( 𝜔 は角周波数) ¨𝑥=−𝜔2⁢𝑥(2) これを式(1)の形で書くには、1階正規形にすればよい: 𝑑𝑑⁢𝑡⁢[𝑥˙𝑥]=[01−𝜔20]⏟__⏟__⏟𝐴⁢[𝑥˙𝑥](3) 𝑎2⁢¨𝑥+𝑎1⁢˙𝑥+𝑎1⁢𝑥=0 のような方程式は、 𝑎𝑖 が 𝑥 およびその微分の関数でなければ、斉次になる。このように、1階正規形の形に統一しておくと、一般論を展開し易くなる。連成振動も斉次方程式である。

他にも、力学編第10章で述べた剛体の回転運動の場合、斉次方程式の形になっている: ˙𝒙=𝝎×⏟𝐴⁢𝒙 𝝎 は角速度ベクトル、 𝒙 は「剛体上に固定した点の、重心を基準とした相対位置」である。

係数行列 𝐴𝑡 が 𝑡 に依存しない場合には、一般的な解法が知られている。この章では、斉次方程式(1)の解法について述べる。

5.1行列の指数関数

第3章で、行列の指数関数について触れたが、もう一度考え直す。

5.1.1初期値問題の解は、初期値に比例する:式(4)

斉次方程式(1)を使って、 𝒙𝑡+𝛿⁢𝑡 を1次近似で表すと 𝒙𝑡+𝛿⁢𝑡≐𝒙𝑡+˙𝒙𝑡⁢⁢𝛿⁢𝑡=𝒙𝑡+𝐴𝑡⁢𝒙𝑡⁢⁢𝛿⁢𝑡=(1+𝐴𝑡𝛿⁢𝑡)⏟__⏟__⏟≡𝛿⁢𝐺𝑡⁢𝒙𝑡 よって、行列 𝛿⁢𝐺𝑡 を掛けると、 𝒙 の時刻が 𝛿⁢𝑡 だけ進むわけである。ここからさらに 𝛿⁢𝑡 進めるには 𝒙𝑡+2⁢𝛿⁢𝑡≐𝛿⁢𝐺𝑡+𝛿⁢𝑡⁡𝒙𝑡+𝛿⁢𝑡≈𝛿⁢𝐺𝑡+𝛿⁢𝑡⁡𝛿⁢𝐺𝑡⁡𝒙𝑡 とすればよい。 𝛿⁢𝐺𝑡+𝛿⁢𝑡 を掛けた形になっている。 ≈ としたのは1次近似ではないからである。値はオイラー法(区分求積法)によるものと同じになる。このように、行列の積を繰り返すことで時間を進めることができるのが、斉次方程式の大きな特徴である。

よって、初期値問題は、 𝑡=0 から任意の時刻 𝑡 まで 𝛿⁢𝐺 を繰り返しかけて、 𝛿⁢𝑡→0 の極限を取ればよい: 𝒙𝑡=(lim𝛿⁢𝑡→0⁡𝛿⁢𝐺𝑡−𝛿⁢𝑡⁡𝛿⁢𝐺𝑡−2⁢𝛿⁢𝑡⋯𝛿⁢𝐺𝛿⁢𝑡⁡𝛿⁢𝐺0)⏟_______⏟_______⏟≡𝐺𝑡⁢𝒙0(4) (区分求積による解と、極限を取る前の値も含めて等しい。式の表現が異なるだけである。)この行列 𝐺𝑡 を時間発展行列と呼ぶことにする。時間発展行列 𝐺𝑡 さえ求めれば、任意の初期値 𝒙0 における解 𝒙𝑡 が、上式のように得られるわけである。よって、「初期値問題を解析的に解く」ことは、「 𝐺𝑡 を求める」ことに他ならない。 𝐺𝑡 の各列は斉次方程式の解になっており、任意の 𝒙0 での解は、それらの線形結合になっていることが分かる。 𝒙0 が 𝒆𝑖 (=第 𝑖 成分だけが 1 でそれ以外が 0 のベクトル)である時の解が、 𝐺𝑡 の 𝑖 列になる。

式(4)を、元の斉次方程式(1)に代入すると、 𝐺𝑡 は以下の方程式から決まることが分かる: ˙𝐺𝑡=𝐴𝑡⁢𝐺𝑡𝐺0=1}(5)

5.1.2係数行列 𝐴𝑡 が定数の場合、時間発展行列 𝐺𝑡 は式(6)

特に、係数行列 𝐴𝑡 が時刻 𝑡 に依存しない場合、 𝛿⁢𝐺𝑡 の 𝑡 によらなくなるので、 𝐺𝑡 は、べき乗で書ける: 𝐺𝑡=lim𝛿⁢𝑡→0⁡(1+𝐴⁢𝛿⁢𝑡)𝑡/𝛿⁢𝑡 右辺において、 𝐴 が行列でなく実数であれば、指数関数 𝑒𝑡⁢𝐴 の定義そのものである。そこでこの記法を、 𝐴 が一般の行列の場合にも使うことにする。即ち、 𝑒𝑡⁢𝐴≡exp⁡(𝐴⁢𝑡)≡lim𝛿⁢𝑡→0⁡(1+𝐴⁢𝛿⁢𝑡)𝑡/𝛿⁢𝑡(6)

なお、式(6)のべき乗を2項展開すると、以下のようになる: 𝑒𝑡⁢𝐴=1+𝑡11!⁢𝐴+𝑡22!⁢𝐴2+𝑡33!⁢𝐴3+⋯ これは常に収束する。これを定義としてもよい。

5.1.3参考 𝐴𝑡 が定数でない場合

係数行列 𝐴𝑡 が時刻 𝑡 に依存する場合、式(6)のように書くことはできない。ただし、 𝛿⁢𝐺𝑡 は、指数関数を使って 𝛿⁢𝐺𝑡=1+𝛿⁢𝑡⋅𝐴𝑡≐𝑒𝐴𝑡⁢𝛿⁢𝑡 と表せるので、式(4)に代入することで、 𝐺𝑡 を、以下のように書くことはできる: 𝐺𝑡=lim𝛿⁢𝑡→0⁡𝑒𝐴𝑡−𝛿⁢𝑡⁢𝛿⁢𝑡⁢𝑒𝐴𝑡−2⁢𝛿⁢𝑡⁢𝛿⁢𝑡⋯𝑒𝐴𝛿⁢𝑡⁢𝛿⁢𝑡⁢𝑒𝐴0⁢𝛿⁢𝑡(7) もし、「任意の時刻 𝑡,𝑡′ に対して 𝐴𝑡,𝐴𝑡′ が可換」であれば、指数部分をまとめることができる: 𝐺𝑡=lim𝛿⁢𝑡→0⁡𝑒(𝐴𝑡−𝛿⁢𝑡+𝐴𝑡−2⁢𝛿⁢𝑡+⋯+𝐴𝛿⁢𝑡+𝐴0)⁢𝛿⁢𝑡=exp⁡(∫𝑡𝑡′=0𝐴𝑡′)(8)

可換でない場合には、式(8)のように書くこともできない。このような場合であっても、 𝐺𝑡 を以下のように表すことがある: 𝐺𝑡=T⁢exp⁡(∫𝑡𝑡′=0𝐴𝑡′)(9) T という記号が付いただけである。 T⁡(⋯) は、時間順序積と呼ばれる線形作用素であり、 (⋯) 部分に含まれる行列の積を、時刻がより小さいものがより右になるように並べ変える操作を表す。例えば T⁡(𝐴𝑡⁢𝐴𝑡′)={𝐴𝑡⁢𝐴𝑡′,𝑡>𝑡′𝐴𝑡′⁢𝐴𝑡,𝑡′>𝑡T⁡(𝐴𝑡+𝐴𝑡′)=𝐴𝑡+𝐴𝑡′ この操作により、式(9)は、式(7)に一致する(行列の積への分解は一意的ではないが、式(9)は一意的に決まる)。

5.2時間発展行列 𝐺𝑡 の例

時間発展行列 𝐺𝑡 の例を、2つ示す。

5.2.1【例1】バネの運動方程式:式(12)

バネの運動方程式(2):(再掲) ¨𝑥=−𝜔2⁢𝑥(10) の初期値問題の解は以下のようになることを既に見た(力学編【2.3-注2】): 𝑥𝑡=𝑥0⁢cos⁡𝜔⁢𝑡+˙𝑥0𝜔⁢sin⁡𝜔⁢𝑡(11)

運動方程式を、1階正規形に変形する1階正規形にしたものは式(3)である:(再掲) 𝑑𝑑⁢𝑡⁢[𝑥˙𝑥]=[01−𝜔20]⏟__⏟__⏟𝐴⁢[𝑥˙𝑥] 式(11)を1階正規形に合わせて書き換えると、時間発展行列 𝐺𝑡 (=𝑒𝐴⁢𝑡) が得られる: [𝑥˙𝑥]=[cos⁡𝜔⁢𝑡𝜔−1⁢sin⁡𝜔⁢𝑡−𝜔⁢sin⁡𝜔⁢𝑡cos⁡𝜔⁢𝑡]⏟_____⏟_____⏟𝐺𝑡⁢[𝑥0˙𝑥0](12) ( 𝐺𝑡 の第2行は、第1章を微分したものである。)この 𝐺𝑡 は、確かに定義式(5):(再掲) ˙𝐺𝑡=𝐴⁢𝐺𝑡𝐺0=1}(13) を満たしている。

あるいは、別の導出もできる。 𝐺𝑡 の2階微分を取ると ¨𝐺𝑡=−𝜔2⁢𝐺𝑡 となる。これは、元の運動方程式(10)において、 𝑥 を 𝐺 で置き換えたものになっている。よって、 𝐺𝑡 は、初期値問題の解(11)に対して同様の置き換えをすればよい: 𝐺𝑡≡𝑒𝐴⁢𝑡=𝐺0⁡cos⁡𝜔⁢𝑡+˙𝐺0𝜔⁢sin⁡𝜔⁢𝑡(14) 𝐺0=1 および ˙𝐺0=𝐴 を代入すれば、式(12)と同じ 𝐺𝑡 が得られる。

上の議論は、すでに分かっている初期値問題の解から、 𝐺𝑡 を構成しているので変則的である。本来は、 𝐺𝑡 を求めることによって、初期値問題を解くことになる。

5.2.2【例2】角速度一定の回転:式(17)

もう1つの例は、剛体の回転である。剛体上に固定された点の(重心を基準とする)相対ベクトル 𝒙 について、その速度 ˙𝒙 は、角速度 𝝎 を用いて ˙𝒙=𝝎×⏟𝐴⁢𝒙 となるのであった(力学編第10章)。この方程式の初期値問題の解は、回転行列 𝑅𝑡 を用いて 𝒙𝑡=𝑅𝑡⁢𝒙0 の形で書けるのであった。 𝑅𝑡 は、時間発展行列に他ならない。

この 𝑅𝑡 が満たすべき微分方程式は ˙𝑅𝑡=𝐴⁢𝑅𝑡=𝝎×𝑅𝑡 である。もう一度微分すると ¨𝑅𝑡=(𝝎×)2⁢𝑅𝑡=−|𝝎|2⁢𝑃⟂⁢𝑅𝑡(15) となる。ただし、 𝑃⟂ は 𝝎 と垂直な方向への正射影行列である: 𝑃⟂≡1−̂𝝎⁢̂𝝎T 式(15)は、2つの方程式に分離できる:( 𝜔≡|𝝎| ) ¨𝑅⁢⟂𝑡=−𝜔2⁢𝑅⟂𝑡¨𝑅⁢∥𝑡=0}(16) ただし、 𝑅⟂𝑡,𝑅∥𝑡 は、以下のように、 𝑅𝑡 に正射影行列を掛けたものである:( 𝑃∥≡1−𝑃⟂ ) 𝑅⟂𝑡≡𝑃⟂⁢𝑅𝑡,𝑅∥𝑡≡𝑃∥⁢𝑅𝑡

式(16)の第1式は、式(14)がそのまま使える形になっているので、以下の解が得られる: 𝑅⟂𝑡=𝑅⟂0⁢cos⁡𝜔⁢𝑡+˙𝑅⁢⟂0𝜔⁢sin⁡𝜔⁢𝑡𝑅∥𝑡=˙𝑅⁢∥0⁢𝑡+𝑅∥0⎫{ {⎬{ {⎭ これに、初期条件 𝑅0=1 と ˙𝑅0=𝝎× を代入して、辺々足し合わせると、 𝑅𝑡 が求まる: 𝑅𝑡≡𝑒𝑡⁢𝝎×=̂𝝎⁢̂𝝎T+(1−̂𝝎⁢̂𝝎T)⁢cos⁡𝜔⁢𝑡+̂𝝎×sin⁡𝜔⁢𝑡(17) これは、力学編第10章で述べたロドリゲスの回転公式である。同章では幾何学的に説明したが、このように導くこともできるわけである。

5.3ジョルダン分解による解法

時間発展行列 𝐺𝑡 を求めることは、一般には困難である。しかし、係数行列 𝐴𝑡 が定数 𝐴 の場合、 𝐺𝑡 の一般的な導出方法が存在する。この節ではこれを述べる。

5.3.1ジョルダン分解による斉次微分方程式の解の公式:式(19)

時間発展行列 𝐺𝑡 を求めるにあたって、解くべき方程式は、式(5)である:(再掲) ˙𝐺𝑡=𝐴⁢𝐺𝑡𝐺0=1}(18) 係数行列 𝐴 は、定数である。

まず、 𝐴 をジョルダン分解(以下の【5.3-注1】)する: 𝐴=𝑃⁢⎡⎢ ⎢⎣𝐽1𝐽2⋱⎤⎥ ⎥⎦⁢𝑃−1 式(18)にこれを代入して、 𝑃−1 を左乗し、 𝑃 を右乗すると、 𝐺′𝑡≡𝑃−1⁢𝐺𝑡⁡𝑃 として ˙𝐺′𝑡=⎡⎢ ⎢⎣𝐽1𝐽2⋱⎤⎥ ⎥⎦⁢𝐺′𝑡𝐺′0=1⎫{ { {⎬{ { {⎭ となる。これの解は、以下で与えられる: 𝐺′𝑡=⎡⎢ ⎢ ⎢⎣𝑒𝐽1⁢𝑡𝑒𝐽2⁢𝑡⋱⎤⎥ ⎥ ⎥⎦ ただし、ジョルダン細胞 𝐽𝑖 の指数関数 𝑒𝐽𝑖⁢𝑡 を求める必要があるが、以下の【5.3-注2】で与えられる。

後は、 𝐺𝑡 の式に戻せば、最終的に以下を得る: 𝐺𝑡=𝑃⁢⎡⎢ ⎢ ⎢⎣𝑒𝐽1⁢𝑡𝑒𝐽2⁢𝑡⋱⎤⎥ ⎥ ⎥⎦⁢𝑃−1(19) ジョルダン分解によって 𝑃,𝐽𝑖 を求めれば、この式によって解が得られるわけである。

【5.3-注1】ジョルダン分解

任意の(実または複素)正方行列 𝐴 は、ある行列 𝑃 を用いて、必ずジョルダン標準形: 𝐴=𝑃⁢⎡⎢ ⎢⎣𝐽1𝐽2⋱⎤⎥ ⎥⎦⁢𝑃−1(20) の形に変形(=ジョルダン分解)できる。ただし、 𝐽𝑖 はジョルダン細胞と呼ばれ、固有値と呼ばれる複素数 𝜆𝑖 を用いて、以下のいずれかの形になる: 𝐽𝑖=𝜆𝑖,[𝜆𝑖10𝜆𝑖],⎡⎢ ⎢⎣𝜆𝑖100𝜆𝑖100𝜆𝑖⎤⎥ ⎥⎦,⎡⎢ ⎢ ⎢ ⎢⎣𝜆𝑖1000𝜆𝑖1000𝜆𝑖1000𝜆𝑖⎤⎥ ⎥ ⎥ ⎥⎦,…(21)

補足

式(20)の形にできることの証明、および、ジョルダン分解の方法については、適当な線形代数の教科書を参照のこと。

全てのジョルダン細胞が 1×1 行列であるとき、 𝐴 は対角化可能であるという。

【5.3-注2】ジョルダン細胞の時間発展行列

ジョルダン細胞 𝐽 に関する時間発展行列 𝑒𝐽⁢𝑡 は、固有値を 𝜆 として以下のように書ける: 𝑒𝐽⁢𝑡=𝑒𝜆⁢𝑡,𝑒𝜆⁢𝑡⁢[1𝑡01],𝑒𝜆⁢𝑡⁢⎡⎢ ⎢ ⎢ ⎢ ⎢⎣1𝑡1!𝑡22!01𝑡1!001⎤⎥ ⎥ ⎥ ⎥ ⎥⎦,𝑒𝜆⁢𝑡⁢⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1𝑡1!𝑡22!𝑡33!01𝑡1!𝑡22!001𝑡1!0001⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,…(22) 𝑒𝜆⁢𝑡 は、 𝜆 が複素数 𝜆=𝛼+𝑖⁢𝜔 の時 𝑒𝜆⁢𝑡=𝑒𝛼⁢𝑡+𝑖⁢𝜔⁢𝑡=𝑒𝛼⁢𝑡⁢𝑒𝑖⁢𝜔⁢𝑡=𝑒𝛼⁢𝑡⁢(cos⁡𝜔⁢𝑡+𝑖⁢sin⁡𝜔⁢𝑡)(23) である。

証明

定義式 𝑑𝑑⁢𝑡⁢𝑒𝐽⁢𝑡=𝐽⁢𝑒𝐽⁢𝑡𝑒𝐽⁢𝑡∣𝑡=0=1 を使って示す。式(22)がこの第2式を満たすことは明らかである。よって、第1式が成り立つことを示せばよい。まず、 𝐽 が 3×3 行列の場合を考える。式(22)の両辺を 𝑒𝜆⁢𝑡 で割ってから微分すると 𝑑𝑑⁢𝑡⁢(左辺)=𝑒−𝜆⁢𝑡⁢(−𝜆⁢𝑒𝐽⁢𝑡+𝑑𝑑⁢𝑡⁢𝑒𝐽⁢𝑡)𝑑𝑑⁢𝑡⁢(右辺)=⎡⎢ ⎢ ⎢⎣01𝑡1!001000⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢⎣010001000⎤⎥ ⎥⎦⁢⎡⎢ ⎢ ⎢ ⎢ ⎢⎣1𝑡1!𝑡22!01𝑡1!001⎤⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢⎣010001000⎤⎥ ⎥⎦⁢𝑒−𝜆⁢𝑡⁢𝑒𝐽⁢𝑡∴𝑑𝑑⁢𝑡⁢𝑒𝐽⁢𝑡=⎛⎜ ⎜ ⎜⎝𝜆+⎡⎢ ⎢⎣010001000⎤⎥ ⎥⎦⎞⎟ ⎟ ⎟⎠⁢𝑒𝐽⁢𝑡=𝐽⁢𝑒𝐽⁢𝑡 青字部分の行列を掛けることにより、掛けられた行列の成分が1つ上に移動することに着目すると、 3×3 以外の場合も同様に成り立つことが分かる。

また、式(23)についても、右辺を微分すると 𝜆⁢𝑒𝜆⁢𝑡 となることから分かる。 ◼

5.4例題

5.4.1【例題】減衰振動の解:【5.4-注1】

バネの振動(2)に、摩擦などによる減衰を加えたものを、減衰振動という。ここでは、速度に比例するような抵抗を受ける場合を考える。運動方程式は ¨𝑥=−𝜔2⁢𝑥⁢−2⁢𝛼⁢˙𝑥⏟抵抗,𝛼≥0(24) である( 𝛼 とせず 2⁢𝛼 としたのは後の式を簡単にするためである)。

これを1階微分方程式の形にすると 𝑑𝑑⁢𝑡⁢[𝑥˙𝑥]=[01−𝜔2−2⁢𝛼]⏟__⏟__⏟≡𝐴⁢[𝑥˙𝑥] となる。この方程式の初期値問題を解くには、上述の議論に従って 𝐴 をジョルダン分解すればよい( 𝜔 と 𝛼 が等しいかどうかで場合分けが必要): 𝐴=⎧{ { {⎨{ { {⎩𝑃≠⁢[𝜆+00𝜆−]⁢𝑃−1≠,𝜔≠𝛼𝑃=⁢[−𝛼10−𝛼]⁢𝑃−1=,𝜔=𝛼 ただし、 𝜆±,𝑃≠,𝑃= は以下で与えられる: 𝜆±≡−𝛼±√𝛼2−𝜔2𝑃≠≡[11𝜆+𝜆−],𝑃=≡[𝛼1−𝛼20]

よって、時間発展行列 𝑒𝐴⁢𝑡 は、式(22)より 𝑒𝐴⁢𝑡=⎧{ { {⎨{ { {⎩𝑒−𝛼⁢𝑡⁢𝑃≠⁢[𝑒√𝛼2−𝜔2𝑡00𝑒−√𝛼2−𝜔2𝑡]⁢𝑃−1≠,𝜔≠𝛼𝑒−𝛼⁢𝑡⁢𝑃=⁢[1𝑡01]⁢𝑃−1=,𝜔=𝛼 となる( 𝜆± を展開した)。後は、 𝑒𝐴⁢𝑡 に含まれる行列の積を計算して、初期値問題の解: [𝑥˙𝑥]𝑡=𝑒𝐴⁢𝑡⁢[𝑥˙𝑥]0 の第1成分を求めれば、 𝑥𝑡 が求まる。実際に計算すると、以下の【5.4-注1】のようになる( 𝜔≠𝛼 については、固有値 𝜆± が実数か複素数かでさらに場合分けする必要があることに注意)。

【5.4-注1】減衰振動の解

減衰振動の運動方程式:( 𝛼≥0 ) ¨𝑥=−𝜔2⁢𝑥−2⁢𝛼⁢˙𝑥 の解 𝑥𝑡 は、初期値 𝑥0,˙𝑥0 のもとで、以下のようになる:

  • 【減衰振動】抵抗が小さい場合: 𝛼<𝜔 𝑥𝑡=𝑒−𝛼⁢𝑡⁢[𝑥0⁢cos⁡(𝜔′⁢𝑡)+𝛼⁢𝑥0+˙𝑥0𝜔′⁢sin⁡(𝜔′⁢𝑡)]𝜔′≡√|𝛼2−𝜔2| 𝑒−𝛼⁢𝑡 項のため、指数関数的に減衰する振動となる。 𝛼=0 とすれば、減衰しないバネの振動になる。
  • 【過減衰】抵抗が大きい場合: 𝛼>𝜔 𝑥𝑡=𝑒−𝛼⁢𝑡⁢[𝑥0⁢cosh⁡(𝛼′⁢𝑡)+𝛼⁢𝑥0+˙𝑥0𝛼′⁢sinh⁡(𝛼′⁢𝑡)]𝛼′≡√𝛼2−𝜔2 ほぼ指数関数的に、振動せず減衰するが、 𝑒−𝛼⁢𝑡 よりは減衰が弱い。
  • 【臨界減衰】中間の場合: 𝛼=𝜔 𝑥𝑡=𝑒−𝛼⁢𝑡⁢[𝑥0+𝑡⁢(𝛼⁢𝑥0+˙𝑥0)] ほぼ 𝑒−𝛼⁢𝑡 に比例して指数関数的に減衰する。振動はしない。減衰振動や過減衰の場合の式で極限 𝛼→𝜔 を取ったものに一致する。