フィボナッチ数列と微分方程式の間

フィボナッチ数列の行列表示

フィボナッチ数列とは、1, 1 から出発して、
 1+1=2、1+2=3、2+3=5、3+5=8 のように、
2つの数字を次々と足し合わせてできる数列のこと。
 ※ 0, 1 から出発して、0+1=1、1+1=2、としても、その先は同じになる。
フィボナッチ数列を Fn という記号で表せば、

  F1 = F2 = 1,
  Fn+2 = Fn+1 + Fn  (n = 1, 2, 3, ・・・)

驚きなのは、このフィボナッチ数列の一般項(n番目の数字)が

   Fn = \frac{1}{\sqrt{5}} \left\{ (\frac{1+\sqrt{5}}{2})^n + (\frac{1-\sqrt{5}}{2})^n \right\}

という式になることだ。
なぜこうなるかは、ググれば出てくると思うので、ここでは少し趣を変えて「行列」という見方をしてみよう。

平面上の (a, b) という点を、(b, a + b) という点に移す2×2行列を考える。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} b \\ a+b \end{pmatrix}

この変換行列をAと名付けておこう。
行列Aで、点(1, 1) からスタートして次々に変換を繰り返せば、フィボナッチ数列ができあがる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \end{pmatrix}
   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 1 \\ 2 \end{pmatrix} = \begin{pmatrix} 2 \\ 3 \end{pmatrix}
   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} 2 \\ 3 \end{pmatrix} = \begin{pmatrix} 3 \\ 5 \end{pmatrix}
    ・・・

変換の様子をグラフにプロットしてみよう。

青い点の折れ線が、(1, 1)からスタートしたフィボナッチ数列だ。
オレンジの折れ線は、試しに(1, 4)からスタートした数列を描いたものだ。
見ての通り、どちらの数列もジグザクを繰り返しつつ、とある一本の直線に漸近している。
実はこの直線の傾きが (1+√5)/2 なのである。この傾きを「黄金比」とも言う。

どのようにして、この直線が出てくるのか。
Aという行列は、
 ・点 (1, 0) を (0, 1) に、
 ・点 (0, 1) を (1, 1) に、
移す変換のことである。
グラフの上では、
 ・赤い実線の矢印を、赤い点線の矢印に、
 ・青い実線の矢印を、青い点線の矢印に、
移す変換となる。(グラフ上では長さ1ではなく、5倍の長さの矢印を描いた。)
この変換によって、ピンク色の正方形は斜め45度で裏返された後、水色の平行四辺形に引き伸ばされた形となる。
この、裏返し->引き伸ばし、裏返し->引き伸ばし・・・が繰り返される様子を想像すれば、
平面全体が、引き伸ばされる方向に流れてゆくイメージが浮かぶだろう。
その流れる先が、黄金比の傾きを持つ直線となるのである。

■ 行列の対角化

この流れる先を計算する方法が「行列の対角化」だ。
対角化とは、ざっくばらんに言えば、行列による変換を回転操作と拡縮操作に分解することだ。
行列Aによる変換を、伸び縮みの方向がちょうど縦横(xy座標軸)となるように、うまい具合に回転する。
(正確に言うと、2本の基底ベクトルをそれぞれ回転して、xy座標軸に重なるようにセットする。
 なので、単なる回転ではなく、ひしゃげた回転とでも言うべきものだ。)
この状態で縦横に(xy座標軸の方向に)拡縮した後、再び逆の回転操作を行って元に戻す。

 (変換A) = (回転P)・(拡縮D)・(逆回転P^-1)


このように、とある変換を回転と拡縮に分解できれば、
変換Aの繰り返しは、以下のように拡縮Dの繰り返しに落とし込むことができる。


  A^n = (P・D・P^-1)・(P・D・P^-1)・(P・D・P^-1)・・・(P・D・P^-1)
    = P・D・(P^-1・P)・D・(P^-1・P)・D・(P・D・P^-1)・・・D・P^-1
    = P・D^n・P^-1

見ての通り、一回前の(逆回転P^-1)と、続く(回転P)が互いにキャンセルし合って、
 (変換A)^n = (回転P)・(拡縮D)^n・(逆回転P^1)
となるわけだ。

拡縮Dをn回繰り返した結果は、少し考えれば分かる。
例えば、横に2倍、縦に3倍拡大する変換をn回繰り返した結果は、
横に 2×2×2・・・ = 2^n倍、縦に3×3×3・・・ = 3^n倍だ。
一般に、横の拡大率をλ1、縦の拡大率をλ2とすれば、

   D^n = \begin{pmatrix} \lambda_1 \quad 0 \\ 0 \quad \lambda_2 \end{pmatrix}^n = \begin{pmatrix} \lambda_1^n \quad 0 \\ 0 \quad \lambda_2^n \end{pmatrix}

※ 縦横の拡縮を行列で書けば対角行列となるので「対角化」と言うわけだ。

問題はどうやって、行列Aをうまい具合に回転と拡縮に分解するか、なのだが、
ここで登場するのが固有ベクトルという概念だ。
固有ベクトルとは、行列の変換によって向きを変えないベクトルのことだ。
上のフィボナッチ数列の例で言えば、点列の流れる先の、黄金比の傾きを持つ直線のことだ。
変換によって向きを変えないベクトルは、それがそのまま平面全体の拡縮の方向を示している。
そのような「向きを変えない不動のベクトル=拡縮方向のベクトル」が見出せれば、
あとはその固有ベクトルをxy座標軸に重なるように持って行くことで、目的を達成できる。

固有ベクトルは必ずしも平面グラフに描けるとは限らないが、複素数まで含めて考えれば、必ずどこかにある。
※ また、複数の固有ベクトルがたまたま重なって一致することもある(重解)。
* 複素数固有ベクトル >> d:id:rikunora:20161208

固有ベクトルの求め方は、やはりググれば出てくるだろうから、簡単に。
固有ベクトルとは、「変換行列Aを施した結果が、やはり元のベクトルのλ倍となっている」ベクトルのことである。

   A \begin{pmatrix} x \\ y \end{pmatrix} = \lambda \begin{pmatrix} x \\ y \end{pmatrix}

0でないベクトル(x, y)に変換を施した結果が0になるためには、変換の行列式が0である必要がある。

   \left\| \begin{pmatrix} a - \lambda \quad b \\ c \quad d - \lambda \end{pmatrix} \right\| = 0
   (a - \lambda)(d - \lambda) - b c = 0   ・・・これを「固有方程式」と言う。

2x2行列の場合、この2次方程式を解くことで、まずλが求まる。
このλのことを固有値と言う。(2次方程式なので、λは2つ出てくる。)
固有値とは、つまり固有ベクトルに沿った拡大率・縮小率のことだ。
得られたλを元の変換に代入して、固有ベクトルが見出せる。
(2つのλに対して、固有ベクトルも2本出てくる。2本の固有ベクトルは互いに直交している。)

こうして得られた固有ベクトルを並べて作った行列が、回転操作Pに相当する。
なぜなら回転操作Pとは、xy座標軸(1,0)と(0,1)を固有ベクトルに重ねる(変換する)操作のことだからだ。
そして、対角行列Dとは、固有値を対角線上に並べた行列になる。
なぜなら対角行列Dは、横にλ1、縦にλ2 拡縮する操作のことだからだ。
最後の逆回転操作 P^-1 は、回転操作Pの逆行列となる。
* 固有ベクトルが直交するのは >> d:id:rikunora:20090307
* 固有ベクトルが直交するのは(2) >> d:id:rikunora:20110203

上の手順でフィボナッチ数列の行列を対角化すれば、結果はこうなる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix}^n = \begin{pmatrix} - \alpha \quad -\beta \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} \beta \quad 0 \\ 0 \quad \alpha \end{pmatrix}^n \left\{ \frac{1}{\alpha-\beta} \begin{pmatrix} -1 \quad -\beta \\ 1 \quad \alpha \end{pmatrix} \right\}
  ここで  \alpha = \frac{1+\sqrt{5}}{2}, \quad \beta = \frac{1-\sqrt{5}}{2} .

一見複雑に見えるが、これを整理すると、最初に書いたフィボナッチ数列の一般項の式に一致する。

以上はフィボナッチ数列の話だったのだが、同じ方法が他の漸化式一般にも使えることは明らかだろう。
例えば、3つの数字を次々と足し合わせてできるトリボナッチ数列であれば、次の行列を対角化することで一般項が求められる。

   \begin{pmatrix} 0 \quad 1 \quad 0 \\ 0 \quad 0 \quad 1 \\ 1 \quad 1 \quad 1 \end{pmatrix}  \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} b \\ c \\ a+b+c \end{pmatrix}

 ※ トリボナッチ数列の場合、固有方程式が3次方程式となるので、一般項はかなり複雑になる。
 ※ 数式処理ソフトMathematicaを使って計算した結果を、記事の最後に貼り付けておいた。

■ 漸化式から微分方程式

さて、漸化式をわざわざ行列に書き直したのは、次につなげたかったからだ。
 『漸化式のステップを、うんと細かくしたもは、微分方程式となる。』
漸化式は整数ステップで進行するが、これを極めて小さなステップで進めれば、微分方程式となる。
つまり漸化式の整数nを、実数tに置き換えたものが微分方程式なのである。
具体的に、フィボナッチ数列を細かく刻んだ微分方程式は、こうなる。

  Fn+2 = Fn+1 + Fn    ・・・フィボナッチ数列
  X(t)''= X(t)'+ X(t)   ・・・対応する微分方程式

この微分方程式の解は、以下の通り。

   X(t) = C_1 \, e^{ \frac{\sqrt{5}+1}{2} t } + C_2 \, e^{ - \frac{\sqrt{5}-1}{2} t }

この解は明らかに、フィボナッチ数列の一般項の類似物だ。

フィボナッチ数列と同じように、微分方程式を行列に書き直してみよう。
X(t) の微分 X(t)' を V(x) という記号に置き換える。
すると2階の微分方程式は、このような2本の連立微分方程式に書き直せる。

      V(t) = X'(t)
  X(t) + V(t) = V'(t)

この2本の式をさらに、行列でまとめれば、こうなる。

   \begin{pmatrix} 0 \quad 1 \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} X(t) \\ V(t) \end{pmatrix} = \frac{d}{dt} \begin{pmatrix} X(t) \\ V(t) \end{pmatrix}

結局のところ、漸化式も(線形)微分方程式も、コアとなる部分は「行列の繰り返し操作」だったのだ。

この微分方程式の解法については、やはりググれば見つかると思うので、繰り返さない。
 ※「定数係数 2階線形同次微分方程式」などでググってみよう。
この微分方程式を解くためには、まず以下のような二次方程式を解く。

   \lambda^2 - \lambda - 1 = 0   ・・・特性方程式

この「微分方程式を代数方程式に置き換えたような式」のことを特性方程式と言うのだが、
この特性方程式の正体は、実は先ほど行列の対角化に見た固有方程式(の同等物)なのである。
コアとなる行列が同じなのだから、解き方も同じなのだ。

特性方程式を解くと、平面全体の拡縮率であるλが得られる。
微分方程式の場合、このλまで求まれば、それで一般解を組み立てることができる。
というのも微分方程式の一般解は、2つの独立な解の組み合わせ(線形結合)でカバーできるからだ。
数列のときに考えた回転操作Pも、所詮は線形結合なので、わざわざ回転させる必要が無い。
どの2つの独立な解を選んでもかまわないのであれば、通常は最も単純な「λ1倍拡縮」と「λ2倍拡縮」の線形結合を選ぶだろう。
 ※ 回転操作は積分定数 C1, C2 が吸収しているわけだ。

なぜ微分方程式の一般解は、線形結合でカバーできるのか。
元を正せば、今行っている操作が全て、ベクトルの足し算と掛け算の組み合わせ(1次式)の範囲内だからだ。
なので、平面であれば、とにかく2本の元になるベクトル(特殊解)さえ見つかれば、あとはその組み合わせで作り出せるわけだ。



■ なぜ  e^\lambda なのか

漸化式と微分方程式を比べたとき、漸化式では、繰り返し操作が対角行列  D^n となっているのに対し、
微分方程式では、特性方程式の解λを  e^{\lambda t} のように、e の累乗とするところが著しく異なっている。
この e (自然対数の底) は、どこから湧いて出てきたのだろうか。

実は、 e という数は、整数ステップを無限小ステップに組み直す過程で立ち現れる数なのだ。
e の素性を知るために、単純な漸化式 An+1 = 2 An を考えてみよう。
1, 2, 4, 8, 16 ・・・ のように、自分自身を2倍にする数列は、
「次の項までの差分が自分自身の値に等しい」という性質を持つ。
  An+1 - An = An
たとえば3番目の項「4」と、その次の「8」の差は、8−4=4となっている。
ただ、これは漸化式が整数ステップだからであって、同じことを小数ステップで行うと結果が違ってくる。

2倍に増える数列を、1年間で100%の利子が付く銀行に例えてみよう。
1年で100%の利子が付く銀行に、1000円の預金を1年間預けたら、2000円になって戻ってくる。
これを半年ステップに分割し、最初の半年で増えた預金を全ていったん引き落とし、再び預け直したら、どうなるか。
 ・まず最初の半年で 1000円×(1+(1/2)) = 1500円に増え、
 ・残りの半年で 1500円×(1+(1/2)) = 2250円となる。
同じように、これを3ヶ月ごとのステップに4分割したならば、
 ・1000×(1+(1/4)) = 1250
 ・1250×(1+(1/4)) = 1562.5
 ・1562.5×(1+(1/4)) = 1953.125
 ・1953.125×(1+(1/4)) = 2441.40625
となる。
1年間をn分割した場合、結果は  \left( 1 + \frac{1}{n} \right)^n 倍に増えることになる。

このnをどこまでも増やした極限  \lim_{n \to \infty} \left( 1 + \frac{1}{n} \right)^n = 2.71828\cdots のことを、我々は e と呼んでいたのである。
e とは、「次の項までの差分が自分自身の値に等しい」数列 An+1 - An = An の、n を無限に小さく刻んだ極限のことだ。
* eについて >> http://miku.motion.ne.jp/precomipo/ (要Flash)

An+1 - An = An に対応する「傾きがその場の値に等しい」微分方程式は、X(t)' = X(t) となる。
この微分方程式の解は  X(t) = C e^t となる。
というより、この微分方程式の答になるような数を探し出してきて、それを e と名付けたのである。

ここまでくれば、なぜ漸化式で  \lambda^n だったものが、
微分方程式 e^{\lambda t} に置き換わるのか、はっきりするのではなかろうか。
整数ステップを無限小ステップに拡張するには、 \lambda^n \left( e^\lambda \right) ^t に置き換える必要がある。
 ※ なお、上の具体例を当てはめると分かるのだが、
 ※ 漸化式の 2^n に対応するものが、微分方程式では (e^1)^t となる。
 ※ つまり漸化式の「λ=2」が、微分方程式の「λ=1」に対応する。
 ※ この違いはどこから来るのかと言えば、微分の X(t)' が、漸化式の差分 An+1 - An に対応するからである。

■ 単振動運動の行列表示

さて、漸化式を細かく刻んで微分方程式ができるなら、逆に微分方程式を細かい漸化式として解釈することも可能となる。
ここでは単振動運動の微分方程式

   X(t)'' = - k X(t)
    (k は振動の固さを表す比例定数、物体の質量 m は1とした)

を例に取り上げてみよう。
この微分方程式を、単純に漸化式に直そうとすると、
  An+2 = - k An
 「現在の値をマイナスにしたものが、2つ後にくる?」
となり、これだけではうまくいかない。
そこで、もう1つ、
 「1つ後の項は、現在の値がそのままコピーされる」※
というルールを付け加えると、漸化式はこのように書ける。
  An+1 = An
  An+2 = - k An
    ※このルールは「慣性の法則」、できるだけ直前の状態を維持することの現れなのだと思う。
この漸化式を行列に直すと、このようになる。

   \begin{pmatrix} 0 \quad 1 \\ -k \quad 0 \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} = \begin{pmatrix} X_2 \\ -k X_1 \end{pmatrix}

この行列を対角化すると、

   \begin{pmatrix} 0 \quad 1 \\ -k \quad 0 \end{pmatrix} = \begin{pmatrix} \frac{i}{\sqrt{k}} \quad -\frac{i}{\sqrt{k}} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -\sqrt{k}i \quad 0 \\ 0 \quad \sqrt{k}i \end{pmatrix} \begin{pmatrix} -\frac{\sqrt{k}i}{2} \quad \frac{1}{2} \\ \frac{\sqrt{k}i}{2} \quad \frac{1}{2} \end{pmatrix}

いちいち  \sqrt{k} が出てくるので、最初から  k = \omega^2 と置いておけば、結果がすっきりする。

   X(t)'' = -\omega^2 X(t)

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad 0 \end{pmatrix} = \begin{pmatrix} \frac{i}{\omega} \quad -\frac{i}{\omega} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -\omega i \quad 0 \\ 0 \quad \omega i \end{pmatrix} \begin{pmatrix} -\frac{\omega i}{2} \quad \frac{1}{2} \\ \frac{\omega i}{2} \quad \frac{1}{2} \end{pmatrix}

この微分方程式の解は、対角行列にある固有値を e の肩に乗せたものの、線形結合となる。

   X(t) = C_1 \, e^{ - \omega i t } + C_2 \, e^{ \omega i t }

固有値虚数の場合、固有ベクトルは実数の平面上には描けない。
その場合、解X(t)の流れは回転となり、結果は正弦波 Sin, Cos の組み合わせとなる。


そして、このように行列に分解してみると、なぜ単振動の角速度が ω^2 のように二乗となるのか、理由が見えてくる。
そもそも2階微分とは「小さく2回変換を繰り返す」ことだったのだから、
1回目の変換でω、2回目の変換でω^2となるわけだ。

力学では、位置Xの2階微分が力に等しい(ニュートンの運動法則)とされている。
これを数列としてとらえるなら、力とは「隣り合う3つの項の差分」ということになる。
  Force = (Xn+2 - Xn+1) - (Xn+1 - Xn)
図に描くと、こういうことだ。

力は、傾きが変化するその点にかかる。

この考え方をもとに、単振動のメカニズムを図解すると、こうなる。


オイラーの公式の行列としての意味

ところで、なぜ固有値虚数の場合は回転となるのだろうか。
理由は「オイラーの公式」でググれば出てくると思うが、ここでも行列流の解釈をしてみよう。

オイラーの公式」とは、 e という数が無限小ステップで立ち現れたのと同じ考え方を、2×2行列にあてはめたものなのである。
90度回転  \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix} を対角化すると、

   \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix} = \begin{pmatrix} -i \quad i \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -i \quad 0 \\ 0 \quad i \end{pmatrix} \begin{pmatrix} -\frac{i}{2} \quad \frac{1}{2} \\ \frac{-i}{2} \quad \frac{1}{2} \end{pmatrix}

90度回転を無限小ステップに刻んだものは、これまでと同じ考え方で、対角行列に現れた固有値を e の肩の上に乗せたものになるだろう。

   \begin{pmatrix} 0 \quad -1 \\ 1 \quad 0 \end{pmatrix}^t = \begin{pmatrix} -i \quad i \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} e^{-i t} \quad 0 \\ 0 \quad e^{i t} \end{pmatrix} \begin{pmatrix} -\frac{i}{2} \quad \frac{1}{2} \\ \frac{-i}{2} \quad \frac{1}{2} \end{pmatrix}

この行列の繰り返し操作をまとめると(そのまま掛け算すると)、こうなる。

   \frac{1}{2} \begin{pmatrix} e^{-i t} + e^{i t} \quad -i ( e^{-i t} - e^{i t} ) \\ i ( e^{-i t} - e^{i t} ) \quad e^{-i t} + e^{i t} \end{pmatrix}

この行列は「90度回転の無限小ステップ」という意味から考えて、小数の角度tの回転変換に一致するはずだろう。

   \begin{pmatrix} cos(t) \quad -sin(t) \\ sin(t) \quad cos(t) \end{pmatrix}

2つの行列の要素を比較して、以下のオイラーの公式を得る。

   cos(t) = \frac{1}{2} ( { e^{-i t} + e^{i t} } )
   sin(t) = \frac{i}{2} ( { e^{-i t} - e^{i t} } )
   e^{i t} = cos(t) + i sin(t)


■ 減衰振動運動の行列表示

単振動が分かったついでに、さらに、抵抗のある減衰振動についても行列化してみよう。

   X(t)'' = -\omega^2 X(t) - 2 r X'(t)
      (ωは角速度、r は抵抗の大きさを表す係数)

この微分方程式を行列に直すと、

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad -2r \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} = \begin{pmatrix} X_2 \\ -\omega^2 X_1 - 2r X_2 \end{pmatrix}

対角化すると、

   \begin{pmatrix} 0 \quad 1 \\ -\omega^2 \quad -2r \end{pmatrix} = \begin{pmatrix} \frac{ -r + \sqrt{r^2-\omega^2} }{\omega^2} \quad -\frac{ +r + \sqrt{r^2-\omega^2} }{\omega^2} \\ 1 \quad 1 \end{pmatrix} \begin{pmatrix} -r + \sqrt{r^2-\omega^2} \quad 0 \\ 0 \quad -r + \sqrt{r^2-\omega^2} \end{pmatrix} \begin{pmatrix} \frac{\omega^2}{2 \sqrt{r^2-\omega^2}} \quad \frac{1}{2} (1 + \frac{r}{ \sqrt{r^2-\omega^2} }) \\ - \frac{\omega^2}{2 \sqrt{r^2-\omega^2}} \quad  \frac{1}{2} (1 - \frac{r}{ \sqrt{r^2-\omega^2} })  \end{pmatrix}

この微分方程式の解は、対角行列にある固有値を e の肩に乗せたものの、線形結合となる。

   X(t) = C_1 \, e^{ -r - \sqrt{r^2-\omega^2} t } + C_2 \, e^{ -r + \sqrt{r^2-\omega^2} t }

ここでも抵抗の係数 r を、なぜ最初から 2r と置くとすっきりするのか、理由が分かってくる。
抵抗は2回の変換について、1回目で r、2回目でもう +r だけかかる。
なので、合わせて 2r となるわけだ。
(ωが掛け算でかかっていたのに対し、抵抗rは足し算でかかる。)

■ まとめ

・漸化式のステップを、うんと細かくしたもは、微分方程式となる。
・漸化式も(線形)微分方程式も、コアとなる部分は「行列の繰り返し操作」
・整数ステップを無限小ステップに拡張するには、 \lambda^n e^{\lambda t} に置き換える。

数列、行列、微分方程式、この3つは合わさって1つの「線形」という世界を形作っている。
というより、もともと世界は1つなので、3者は1つの世界を別の側面から切り取った姿なのである。
大学や専門学校では1つ1つを学ぶのに息切れして、なかなか3つ合わせた1つの世界までたどり着けない。
それでも線形世界は、1つだけ、あるいはバラバラなまま終わってしまうには余りにも惜しいほどの視点を提供する。
3つは1つなのだと心してほしい。

■ おまけ:トリボナッチ数列の対角化