ぱたへね

はてなダイアリーはrustの色分けができないのでこっちに来た

差分方程式

courseraでPractical Time Series Analysisという時系列データを分析する講義を受けています。時系列のデータが、過去のデータから導かれるとき、順番に計算せずに一般解を求めるやり方があります。

この計算方法を知らず、授業を聞いていてびっくりしたのでまとめました。

問題

解くべき問題は、プログラミングの演習なんかでよく見る漸化式です。

a_0 = 3
a_1 = 8
a_n = 5a_{n-1} - 6a_{n-2}

gaucheで書くとこんな感じ。

(define diff-eq
  (lambda (n)
    (cond
     ((eq? n 0) 3)
     ((eq? n 1) 8)
     (else 
      (+ (* 5 (diff-eq (- n 1)))
         (* -6 (diff-eq (- n 2))))))))
      
(map diff-eq (iota 10))

10項まで計算するとこうなります。

(3 8 22 62 178 518 1522 4502 13378 39878)

簡単ですね。

auxiliary equation, characteric equation

aをλに置き換えて、表記を変えます。

a_n = \lambda ^n
\lambda^n = 5\lambda ^{n-1} - 6\lambda^{n-2}

最後の式が、auxiliary equation, characteristic equationと呼ばれるそうです。日本語だと特性方程式でしょうか。その式のnを2ずらし、二次方程式と見なして解を求めます。
\lambda^2 - 5\lambda + 6 = 0

これをλについて解くと、λ=2, λ=3になります。この2と3を使って
a_n = c_12^n + c_2 3^n

の関係式を作ります。あとは、a0=3, a1=8を使ってc1, c2を求めると、c1=1, c2=2になり、これがanの一般解になります。
a_n = 2^n + 2 * 3^n

試してみた

このアルゴリズムが初見だったので、いまいち納得いかずgaucheで実装して比べてみました。

(define diff-eq-2o
  (lambda (n)
    (+ (expt  2 n)
       (* 2 (expt  3 n)))))

(map diff-eq-2o (iota 10))

これで実行すると、全く同じ結果になります。

(3 8 22 62 178 518 1522 4502 13378 39878)

すっごーい

まとめ

差分方程式を使うことで、二次の漸化式について最初から計算しなくても好きな値を計算する方法があります。今回は二次の漸化式でしたが、授業では任意のk次についても説明があります。もちろんみんな大好きフィボナッチも解くことができます。