斉次方程式
斉次(常微分)方程式とは、ベクトル
𝒙𝑡
に対する、以下の形の(あるいはこの形に変形できる)微分方程式である:
˙𝒙𝑡=𝐴𝑡𝒙𝑡(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)のように書くこともできない。このような場合であっても、
𝐺𝑡
を以下のように表すことがある:
𝐺𝑡=Texp(∫𝑡𝑡′=0𝐴𝑡′)(9)
T
という記号が付いただけである。
T(⋯)
は、時間順序積と呼ばれる線形作用素であり、
(⋯)
部分に含まれる行列の積を、時刻がより小さいものがより右になるように並べ変える操作を表す。例えば
T(𝐴𝑡𝐴𝑡′)={𝐴𝑡𝐴𝑡′,𝑡>𝑡′𝐴𝑡′𝐴𝑡,𝑡′>𝑡T(𝐴𝑡+𝐴𝑡′)=𝐴𝑡+𝐴𝑡′
この操作により、式(9)は、式(7)に一致する(行列の積への分解は一意的ではないが、式(9)は一意的に決まる)。
5.2時間発展行列 𝐺𝑡 の例
時間発展行列
𝐺𝑡
の例を、2つ示す。
5.2.1【例1】バネの運動方程式:式(12)
バネの運動方程式(2):(再掲)
¨𝑥=−𝜔2𝑥(10)
の初期値問題の解は以下のようになることを既に見た(力学編【2.3-注2】):
𝑥𝑡=𝑥0cos𝜔𝑡+˙𝑥0𝜔sin𝜔𝑡(11)
運動方程式を、1階正規形に変形する1階正規形にしたものは式(3)である:(再掲)
𝑑𝑑𝑡[𝑥˙𝑥]=[01−𝜔20]⏟__⏟__⏟𝐴[𝑥˙𝑥]
式(11)を1階正規形に合わせて書き換えると、時間発展行列
𝐺𝑡
(=𝑒𝐴𝑡)
が得られる:
[𝑥˙𝑥]=[cos𝜔𝑡𝜔−1sin𝜔𝑡−𝜔sin𝜔𝑡cos𝜔𝑡]⏟_____⏟_____⏟𝐺𝑡[𝑥0˙𝑥0](12)
(
𝐺𝑡
の第2行は、第1章を微分したものである。)この
𝐺𝑡
は、確かに定義式(5):(再掲)
˙𝐺𝑡=𝐴𝐺𝑡𝐺0=1}(13)
を満たしている。
あるいは、別の導出もできる。
𝐺𝑡
の2階微分を取ると
¨𝐺𝑡=−𝜔2𝐺𝑡
となる。これは、元の運動方程式(10)において、
𝑥
を
𝐺
で置き換えたものになっている。よって、
𝐺𝑡
は、初期値問題の解(11)に対して同様の置き換えをすればよい:
𝐺𝑡≡𝑒𝐴𝑡=𝐺0cos𝜔𝑡+˙𝐺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)がそのまま使える形になっているので、以下の解が得られる:
𝑅⟂𝑡=𝑅⟂0cos𝜔𝑡+˙𝑅⟂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
のもとで、以下のようになる:
- 【減衰振動】抵抗が小さい場合:
𝛼<𝜔
𝑥𝑡=𝑒−𝛼𝑡[𝑥0cos(𝜔′𝑡)+𝛼𝑥0+˙𝑥0𝜔′sin(𝜔′𝑡)]𝜔′≡√|𝛼2−𝜔2|
𝑒−𝛼𝑡
項のため、指数関数的に減衰する振動となる。
𝛼=0
とすれば、減衰しないバネの振動になる。
- 【過減衰】抵抗が大きい場合:
𝛼>𝜔
𝑥𝑡=𝑒−𝛼𝑡[𝑥0cosh(𝛼′𝑡)+𝛼𝑥0+˙𝑥0𝛼′sinh(𝛼′𝑡)]𝛼′≡√𝛼2−𝜔2
ほぼ指数関数的に、振動せず減衰するが、
𝑒−𝛼𝑡
よりは減衰が弱い。
- 【臨界減衰】中間の場合:
𝛼=𝜔
𝑥𝑡=𝑒−𝛼𝑡[𝑥0+𝑡(𝛼𝑥0+˙𝑥0)]
ほぼ
𝑒−𝛼𝑡
に比例して指数関数的に減衰する。振動はしない。減衰振動や過減衰の場合の式で極限
𝛼→𝜔
を取ったものに一致する。