ジョルダン標準化の応用
今回はジョルダン標準化の応用です。
ジョルダン標準化が一般で、対角化は特殊です。
だから、ジョルダン標準化のステップをふむことの中に対角化は含まれてしまいます。
また、ジョルダン標準化で行列のベキ計算がカンタンになります。
それに、数列、差分、微分などの連立方程式を行列に直して処理するときも、標準化は見通しを与えてくれます。
なお、ｋ次のJB（ジョルダンブロック）はｋ乗するとゼロ行列になるベキゼロ行列です。
＜対角化＞
「n次の行列Aの固有多項式が互いに素なn個の因数の積に分解できればよい。」
n個の異なる固有値が対角成分として並ぶと、Jの直和分解の成分が自動的に次数1になるしかないから、
当然の条件と言えるね。
ｎ個の固有値に対する固有ベクトルがすべて独立しているということになる。
Jが固有空間（箱）がすべて線状になるように、全体空間を分割する行列になっているということだね。
＜行列のベキ＞
Aの変換行列をPとすると、A^nの計算は、
J＝P^(-1)APから、
A＝PJP^(-1)なので、A^n=PJ^nP^(-1)と変形できるけど、
J=S＋N
のように
Jは、対角化行列（伸び縮み）＋情報伝播と分解できるイメージだったことを思い出そう。
N部分はベキゼロで減衰するので、n乗の過程でほぼS部分の骨格だけが繰り返されるという構造だった。
それに、「行列〜表のデータを取り出したもの」
くわしくはこちらにもあるように、
A＝λE＋Nとおくと、
λEが単位行列の倍数だからどんな行列相手でもかけ算の順番を入れ替えることができる。
だから、２項定理が使える。
N= {{0,1,0},
{0,0,1},
{0,0,0}},
N^2={{0,0,1},
{0,0,0},
{0,0,0}},
N^3=N^4=....=O
A^n=(λE)^n + nC1(λE)^{n-1}N + nC2(λE)^{n-2}N^2 +nC3(λE)^{n-3}N^3 +.....
=λ^n E + nλ^{n-1} N + n(n-1)/2λ^{n-2} N^2
={{λ^n, nλ^{n-1}, n(n-1)/2 λ^{n-2} },
{0 , λ^n, nλ^{n-1} },
{0 , 0 , λ^n }}
Nの3乗以降が減衰によってゼロ行列になるため、
たったの3次行列に収まってしまっていたよね。
＜線形差分方程式＞
数列 x(m)(m=O, 1, ，…）に関する方程式
x(m+n) = a_0 x(m）＋ a_1 x(m+l）＋… a_{n-1} x(m+n-l)
(m=O, 1, 2,…)
をn階の1変数線形差分方程式という。
たとえば、2階の線形差分方程式
x(m+2)=x(m)+x(m+1) (m=0,1,2,...)において、
初期値x(0)=0,x(1)=1を与えると、
{x(i)}={0,1,1,2,3,5,8,13,...}というフィボナッチ数列を表すことがわかるね。
また、1階の3変数線形差分方程式は
x1(m+1)=a_11 x1(m)+a_12 x2(m)+a_13 x3(m)
x2(m+1)=a_21 x1(m)+a_22 x2(m)+a_23 x3(m)
x3(m+1)=a_31 x1(m)+a_32 x2(m)+a_33 x3(m) (m=0,1,2,...)
m=0のときのxi(m)の値を初期値という。
この式のイメージは状態xiが３つあり、それぞれの遷移確率がa_jiで表されているときの、先の遷移状態を求めたりすることができる。
この係数を取り出すとカンタンに行列ができるから、
行列方程式や、行列のｎ乗の置き換えられるはずだ。
だから、それをジョルダン標準化につなげられる。
最初のn階の1変数線形差分方程式
x(m+n) = a_0 x(m）＋ a_1 x(m +1）＋… a_{n-1} x(m +n-1)
は、ただの多項式に見えるが、
x_1(m)=x(m)
x_2(m)=x(m+1)
....
x_n(m)=x(m+n-1)
のようにｘの引数をxの番号としての添え字に移し替えると、
x_{n+1}(m) = a_0 x_1(m）＋ a_1 x_2(m）＋… a_{n-1} x_n(m)
のように多項漸化式に変わる。
ここで、
X(m)=t[x_1(m),x_2(m),...x_n(m)]とおくと、
X(m+1)=t[x_2(m),x_3(m),...x_n+1(m)]だから、
どの成分もまとめて、１つ先に進むのだけど、最後のx_n+1(m)だけは、
x_{n+1}(m)=[a_0,a_1,...a_{n-1}] X(m)
と内積になる。
ということは、
A={{0,1,0,0,....,0},
{0,0,1,0,....,0},
.....
{0,0,.......,1,0},
{a_0.a_1,...,a_{n-1}}
とおくと、
X(m+1)=AX(m)
というベクトルX(m)の漸化式ができるね。
この行列Aをコンパニオン（友伴）行列といいます。
だから、ベクトルX(0)を初期値として、X(m)=A^m X(0)によって、
最初の解空間Wを定めることが、最初のn階の1変数線形差分方程式を解くことになる。
ここで、ジョルダン標準化によって、A^mの計算を見通しよくすることができるはずだ。
X_k = A^k x_0 = P J^k P^{-1} X_0
ジョルダン行列のｎ乗では、対角成分がλ^nで、
そのとなりが nλ^{n-1}、右上かどがn(n-1)/2 λ^{n-2}の上三角行列だった。
対角化できれば、対角成分だけになるので、λ^nの指数的な爆発だけを引き起こすはずだ。
しかし、対角化できないときはnλ^{n-1}成分が効いてくる。
λ＜1ならば、nλ^{n-1}は一時的に大きくなる（過渡的な増加）ものの、
最後は λ^n の衰退パワーに負けて0に収束します。
λ=1 のときはλ^n = 1なので、n 1^{k-1} = n となり、時間経過とともに状態が直線的（線形）に増大し続けるという現象が起きます。
ジョルダン標準化の行列のｎ乗の上三角行列であること、対角部分以外があることによって、
爆発か減衰かという骨格部分の挙動だけではない、
過渡的現象や直線現象という中間的な動きを表現できます。
＜行列の極限・級数・微分＞
Aの最小多項式φ(x)の固有値λiのリストΣAをAのスペクトルという。
多項式の因数(x-λi)の指数がmiだとする。iが1からs。m=m1+m2+...+ms
関数列f1(x),f2(x),...,ｆp(x),...に対して、kを1からsまで動かすとき、p→∞のとき、fp(λk),f'p(λk),f''p(λk),...f^(mk-1)p(λk)の極限値がすべてあるとき、fp(x)はΣAで極限値を持つという。
関数f(x)について、p→∞のとき、fp(λk),f'p(λk),f''p(λk),...f^(mk-1)p(λk)の極限値がすべてあるとき、f(x)に収束するといい、fp(ΣA)→f(ΣA)とかく。
p→∞のとき、fp(ΣA)→f(ΣA)ならば、
p→∞のとき、fp(A)→f(A)である。
このように、１つの変数ｘの関数と同様に収束が定義できるので、
１つの行列Aの関数についても収束が定義できて、級数展開も定義できるね。
指数関数のベキ級数
e^x = Σ1/p! x^p
の収束半径が∞だから、xをAに置き換えて、
どんな正方行列Aでもe^A=Σ1/p! A^pが成り立ち、これを行列Aの指数関数と呼ぶ。
ちなみに、e^0=E（Oがゼロ行列、Eは単位行列）
行列に極限が定義できるのだから、微分も定義できるはず。
パラメータtの成分がならぶ行列A(t)={aij(t)}がm行n列だとする。d/dt aij(t)を(i,j)成分とする行列のことを行列A(t)の微分といい、d/dt A(t)とかく。
たとえば、F(t)=e^{tA}の微分を考えてみよう。
e^{tA}=Σ1/p! t^p A^p のΣされる(i,j)成分をt^p aij^pとかこう。
e^{iA}の(i,j)成分はそのΣだから、aij(t)＝Σt^p aij^pとかける。
だから、その微分は項別微分により、
d/dt aij(t)=Σpt^{p-1} aij^p .
1/p!A^pをaij^pとおいたので、p-1=qとおくと
d/dt e^{tA}=Σpt^{p-1} 1/p!A^p＝AΣ1/q! t^q A^q
=Ae^{tA}
他にもd/dtを'印でかくと
(A(t)B(t))'=(A(t))'B(t)+A(t)(B(t))'
と積の微分法則が行列でも成り立つ。
＜微分方程式＞
連立微分方程式
d/dt x1(t)=a_11 x1(t)+a_12 x2(t)....+a_1n xn(t),
d/dt x2(t)=a_21 x1(t)+a_22 x2(t)....+a_2n xn(t),
....................
....................
d/dt xn(t)=a_n1 x1(t)+a_n2 x2(t)....+a_3n xn(t)
x(t)=t[x1(t),....,xn(t)]とおき、aijを(i,j)成分とすると、
d/dt x(t)=Ax(t)と行列の微分方程式で、初期値はx(0)となるね。
この解x(t)の全体Wがベクトル空間となるから、
これが解空間（解箱）だ。さっきの
積の微分法則と指数の微分から、
x(t)=e^{tA}x(0)
となるね。
ジョルダン標準形Jが求まれば、A＝PJP^(-1)。
指数関数のテイラー展開はe^{tA}=Σ1/p! t^p A^p だった。
ここで各項のA^pをPJ^pP^(-1)におきかえよう。
e^{tA}=Σ1/p! t^p A^p
＝Σ1/p! t^p (PJ^pP^(-1))
＝P(Σ1/p! t^p J^p)P^(-1)
このカッコの中はe^{tA}のAがJに置き換わった形だね。
だから、
e^{tA}=Pe^{tJ}P^(-1)となる。
Aのジョルダン標準化
J=P^(-1)AP=JB1⊕JB2⊕...⊕JBr
だから、
P^(-1)e^{tA}P
=e^{t(P^(-1)AP}=e^{tJB1}⊕e^{tB2}...⊕e^{tBr}=e^{tJ}
最後にe^{tA}=P(e^{tJ})P^(-1)となるので、これに初期値x(0)をかければ、解の箱が決まるね。
ところで、Jと指数関数は仲がよいことに気づきまたか？
J=λE+Nと和分解できるね。
ところがJはeの指数にくるから、和が積に直せる。
e^{tJ}=e^t{λE+N}=e^{tλE+tN}
=e^λtE * e^tN となるね。
e^{λtE}は単位行列Eのe^λt倍だから、対角成分がe^λtだ。この成分は、λに応じた回転、振動を表しているね。
一方でe^tNはどうだろうか。
N^3=Oだから、e^{tN}のテイラー展開がN^2で止まる。
e^{tA}=Σ1/p! t^p A^pだったので、
e^{tN}=1/0! t^0 N^0 +1/1! t^1 N^1 +1/2! t^2 N^2
= I + t N + 1/2 t^2 N^2
だから、対角成分が１、となりがt, 右上角が1/2 t^2
ではいよいよ積の計算だ。
e^{tJ}=e^{λtE}*e^{tN}は、e^{tN}の成分がe^λt倍される。
e^{tJ}=
{{e^λt, t e^λt, 1/2! t^2 e^λt},
{{ 0, e^λt, t e^λt},
{{ 0, 0, e^λt}}
最後にe^{tA}=P(e^{tJ})P^(-1)となるので、
これに初期値x(0)をかければ、解の箱が決まる。
e^{tJ}の対角成分指数関数なので、λに応じた回転に応じた振動・波でしかない。
しかし、残りの三角成分にｔの多項式がかけられている分、変化の仕方に加速、減速の要素が加味される。
それが、ばねの共振や臨界減衰を作りだしていることがわかる。
物理現象の裏にジョルダン細胞のN成分あり
というところだね。
＜Geogebraでコード化＞
// 固有値 λをlamとして、0.7を初期値にする。
// 右上成分 N (0: OFF / 1: ON)の切り替えにswの01をアニメーションにする。
// 初期値ベクトル (x1_0, x2_0)としてc1,c2を設定する。
// 微分方程式の解f(x)= 指数関数×（定数＋多項式）にして、多項式部分がswでONOFF切り替えになる。
// 差分方程式の点列ptsも同様に作る。
//表示用のテキストをif文で作ることによって、swの変化にあわせた表示になる。
lam = Slider(0.1, 1.5, 0.05)
sw = Slider(0, 1, 1)
c1 = 0
c2 = 1
f(x) = exp(-lam * x) * (c1 + sw * x * c2)
pts = Sequence(Point({k, exp(-lam * k) * (c1 + sw * k * c2)}), k, 0, 15)
If(sw ≟ 1, "{sw=1ジョルダン型では過渡的な「山」が発生}", "{対角型 では単調減衰}")