(スマホからだと数式が表示できない場合があるようなので、その場合はPCなどからご覧ください。申し訳ありません。)
トリボナッチ数列とは？
フィボナッチ数列は、直前の2つの項を足して次の項を作る数列です。式で表すと、
\begin{cases}
a_{n+2} = a_{n+1} + a_{n}\\ a_0 = 0\\a_1 = 1
\end{cases}
となります。これを拡張して直前の3つの項を足した数を次の項にしたものがトリボナッチ数列です。式で表すと、
\begin{cases}
a_{n+3} = a_{n+2} + a_{n+1} + a_{n}\\ a_0 = 0\\a_1 = 0\\a_2 = 1
\end{cases}
となります。少し計算してみると、
a_3 = 1+0+0=1\\
a_4 = 1+1+0=2\\
a_5 = 2+1+1=4\\
a_6 = 4+2+1=7\\
\vdots
となります。
前準備
漸化式を満たす実数の数列のなす空間を$V$とします。数列${ a_n }$の最初の3項 $a_0,a_1,a_2$ が与えられているので、$n\geq3$において数列 ${ a_n }$ は一意に定まります。
\boldsymbol{u}=\{1,0,0,1,\dots \},\ \boldsymbol{v}=\{0,1,0,1,\dots \},\ \boldsymbol{w}=\{0,0,1,1,\dots\}
とします。$c_1,c_2,c_3$に対して、$c_1\boldsymbol{u}+c_2\boldsymbol{v}+c_3\boldsymbol{w}=\boldsymbol{o}$が成り立ったとすると、
c_1\{1,0,0,\dots \}+c_2\{0,1,0,\dots \}+c_3\{0,0,1,\dots \}=\{c_1,c_2,c_3,\dots \}=\{0,0,0,\dots \}
となるので、$c_1=c_2=c_3=0$です。よって、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$は一次独立とわかります。
次に、
\boldsymbol{a}= \{ a_0,a_1,a_2,\dots \}
を$V$の勝手な元とすると、
\begin{align}
\boldsymbol{a}&=\{ a_0,0,0,\dots \}+\{0,a_1,0,\dots \} +\{0,0,a_2,\dots \} \\
&=a_0\{1,0,0,\dots \}+a_1\{0,1,0,\dots \}+a_2\{0,0,1,\dots \} \\
&=a_0\boldsymbol{u}+a_1\boldsymbol{v}+a_2\boldsymbol{w}
\end{align}
となり、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ の一次結合で表せます。よって、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ は$V$を生成します。
以上より、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ は一次独立かつ$V$を生成するので、$V$の基底となります。
ここで、写像 $f:V\rightarrow V$を
f(\{ a_n\}_{n=0}^{\infty})=\{ a_n\}_{n=1}^{\infty}
とします。$\boldsymbol{a}={ a_0,a_1,a_2,\dots } \in V$のとき、$f(\boldsymbol{a})={ a_1,a_2,a_3,\dots }$も漸化式を満たすので、$f(\boldsymbol{a})\in V$です。
\boldsymbol{a}=\{ a_n\}_{n=0}^{\infty}\in V \\
\boldsymbol{b}=\{ b_n\}_{n=0}^{\infty}\in V \\
とすると、定数 $c$ に対して、
f(\boldsymbol{a}+\boldsymbol{b})=f(\{ a_n+b_n\}_{n=0}^{\infty})=\{ a_n+b_n\}_{n=1}^{\infty}=\{ a_n\}_{n=1}^{\infty}+\{ b_n\}_{n=1}^{\infty}=f(\boldsymbol{a})+f(\boldsymbol{b})\\
f(c\boldsymbol{a})=f(c\{ a_n\}_{n=0}^{\infty})=c\{ a_n\}_{n=1}^{\infty}=cf(\boldsymbol{a})
が成り立つので、$f$は$V$の線形変換です。
漸化式を行列で表現
$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$に関して、
\begin{align}
f(\boldsymbol{u})&=\{0,0,1,\dots \}=\boldsymbol{w}\\
f(\boldsymbol{v})&=\{1,0,1,\dots \}=\boldsymbol{u}+\boldsymbol{w}\\
f(\boldsymbol{w})&=\{0,1,1,\dots \}=\boldsymbol{v}+\boldsymbol{w}
\end{align}
なので、
[f(\boldsymbol{u})\quad f(\boldsymbol{v})\quad f(\boldsymbol{w})]=
[\boldsymbol{u}\ \boldsymbol{v}\ \boldsymbol{w}]
\begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 1 & 1
\end{bmatrix}
と表せます。よって、$f$の基底 $\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$に関する表現行列は、
\begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 1 & 1
\end{bmatrix}
です。この表現行列を$A$とします。
固有値と公比の関係
\boldsymbol{p}=\{ r^{n} \}_{n=0}^{\infty}=\{1,r,r^2,\dots \}
が$V$に属したとすると、
f(\boldsymbol{p})=f(\{ r^{n} \}_{n=0}^{\infty})=\{ r^{n} \}_{n=1}^{\infty}=\{ r^{n+1} \}_{n=0}^{\infty}=r\{ r^{n} \}_{n=0}^{\infty}=r\boldsymbol{p}
より、$\boldsymbol{p}$は固有値$r$の固有ベクトルになります。逆に、上の式から$f$の固有値 $r$の固有ベクトルは、公比$r$の等比数列になることも分かります。よって、公比と固有値は等しいとわかります。
固有値を求める
$f$の固有値を求めます。$A$の固有多項式は、$I$を単位行列として、
\begin{align}
\varphi_A(\lambda)&=|A-\lambda I|\\
&=\begin{vmatrix}
0-\lambda & 1 & 0 \\
0 & 0-\lambda & 1 \\
1 & 1 & 1-\lambda
\end{vmatrix}\\
&=-\lambda^3+\lambda^2+\lambda+1
\end{align}
となります。しかし、$\varphi_A(\lambda)=0$は簡単には解くことはできません。3つの解をそれぞれ $\alpha,\beta,\gamma$とします。$\alpha,\beta,\gamma$は、あとの章で求めることにします。解 $\alpha,\beta,\gamma$に対応する固有ベクトルである等比数列の
\boldsymbol{a}=\{ \alpha^{n} \}_{n=0}^{\infty}\\
\boldsymbol{b}=\{ \beta^{n} \}_{n=0}^{\infty}\\
\boldsymbol{c}=\{ \gamma^{n} \}_{n=0}^{\infty}
は$V$に属します。これらは $\alpha,\beta,\gamma$が全て異なるとすると、$f$の相異なる固有値に対応する固有ベクトルなので一次独立です。また、$\dim V=3$です。よって、$\boldsymbol{a},\boldsymbol{b},\boldsymbol{c}$は$V$の基底となります。
係数の決定
以上より、ある $c_1,c_2,c_3$が存在して、
\boldsymbol{a_n}=c_1\boldsymbol{a}+c_2\boldsymbol{b}+c_3\boldsymbol{c}
と表せるので、
\begin{align}
\boldsymbol{a_n}&=c_1\boldsymbol{a}+c_2\boldsymbol{b}+c_3\boldsymbol{c}\\
\Leftrightarrow \{a_0,a_1,a_2,\dots \}&=c_1\{ \alpha^{n} \}_{n=0}^{\infty}+c_2\{ \beta^{n} \}_{n=0}^{\infty}+c_3\{ \gamma^{n} \}_{n=0}^{\infty}\\
\Leftrightarrow \{0,0,1,\dots \}&=c_1\{ 1,\alpha,\alpha^2,\dots \}+c_2\{1,\beta,\beta^2,\dots \}+c_3\{ 1,\gamma,\gamma^2,\dots\} \\
\Leftrightarrow \{0,0,1,\dots \}&=\{ c_1+c_2+c_3,c_1\alpha+c_2\beta+c_3\gamma,c_1\alpha^2+c_2\beta^2+c_3\gamma^2, \dots\}
\end{align}
これを行列で表現すると、
\begin{bmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{bmatrix}
\begin{bmatrix}
c_1 \\ c_2 \\ c_3
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0 \\ 1
\end{bmatrix}
です。左の行列の行列式は、
\begin{align}
\begin{vmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{vmatrix}
&=
\begin{vmatrix}
1 & 0 & 0 \\
\alpha & \beta-\alpha & \gamma-\alpha \\
\alpha^2 & \beta^2-\alpha^2 & \gamma^2-\alpha^2
\end{vmatrix}\\
&=
\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
\beta^2-\alpha^2 & \gamma^2-\alpha^2
\end{vmatrix}\\
&=\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
\beta^2-\alpha^2 - (\beta-\alpha)(\beta + \alpha) & \gamma^2-\alpha^2 - (\gamma-\alpha)(\beta + \alpha)
\end{vmatrix}\\
&=\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
0 & (\gamma-\alpha)(\gamma+\alpha) - (\gamma-\alpha)(\beta + \alpha)
\end{vmatrix}\\
&=(\beta-\alpha)(\gamma-\alpha)(\gamma - \beta)\\
&=(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
\end{align}
となります。これはヴァンデルモンド行列式としても知られています。クラメルの公式を用いると、
\begin{align}
c_1 &= \frac{
\begin{vmatrix}
0 & 1 & 1 \\
0 & \beta & \gamma \\
1 & \beta^2 & \gamma^2
\end{vmatrix}
}{
\begin{vmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{vmatrix}
}\\
&=
\frac{
-\begin{vmatrix}
1 & \beta^2 & \gamma^2 \\
0 & \beta & \gamma \\
0 & 1 & 1
\end{vmatrix}
}{
(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=
\frac{
-\begin{vmatrix}
\beta & \gamma\\
1 & 1
\end{vmatrix}
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=\frac{
\gamma-\beta
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=\frac{1}{(\alpha-\beta)(\alpha-\gamma)
}
\end{align}
と求められます。同様に、
c_2 = \frac{
\alpha-\gamma
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}=\frac{1
}{(\beta-\alpha)(\beta-\gamma)
}\\
c_3 = \frac{
\beta-\alpha
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}
=\frac{1
}{(\gamma-\alpha)(\gamma-\beta)
}
と求められます。
以上より、
\begin{align}
\boldsymbol{a_n} &= c_1 \boldsymbol{a} + c_2 \boldsymbol{b} + c_3 \boldsymbol{c}\\
&=\left\{ \frac{\alpha^n}{(\alpha-\beta)(\alpha-\gamma)} + \frac{\beta^n}{(\beta-\alpha)(\beta-\gamma)} + \frac{\gamma^n}{(\gamma-\alpha)(\gamma-\beta)}\right\}_{n=0}^{\infty}
\end{align}
なので、
a_n = \frac{\alpha^n}{(\alpha-\beta)(\alpha-\gamma)} + \frac{\beta^n}{(\beta-\alpha)(\beta-\gamma)} + \frac{\gamma^n}{(\gamma-\alpha)(\gamma-\beta)}
が導けました。
具体的な数値の計算
以上までで、一般項の形はわかりましたが、具体的な $\alpha,\beta,\gamma$の値はわかっていません。3次方程式の解の公式を使って求める方法もありますが、今回はPythonで求めてみましょう。$a[i]$をトリボナッチ数列の定義通りに計算したもの、$b[i]$を求めた一般項で計算したものとします。sympyを使って代数的に解く方法もありますが、今回はnumpyを使って数値的に求めます。
ソースコード
import numpy as np
import warnings
warnings.filterwarnings('ignore') # 複素数の計算時にWarningが発生するので消しておく
alpha, beta, gamma = np.roots([1,-1,-1,-1]) # solve 1*x^3 + (-1)*x^2+ (-1)*x + (-1) = 0
print("alpha =", alpha)
print("beta =", beta)
print("gamma =", gamma)
a = np.zeros(101)
b = np.zeros(101)
a[2] = 1
print(" i "," a[i] "," b[i] ", " a[i]-b[i] ")
for i in range(101):
if i > 2:
a[i] = a[i-1] + a[i-2] + a[i-3]
b[i] = np.power(alpha,i)/((alpha-beta)*(alpha-gamma)) + np.power(beta,i)/((beta-alpha)*(beta-gamma)) + np.power(gamma,i)/((gamma-alpha)*(gamma-beta))
print('{:3}'.format(i), '{:26}'.format(int(a[i])), '{:26}'.format(int(round(abs(b[i])))), '{:12}'.format(int(a[i]) - int(round(abs(b[i])))))
出力結果
alpha = (1.839286755214161+0j)
beta = (-0.41964337760708065+0.6062907292071997j)
gamma = (-0.41964337760708065-0.6062907292071997j)
i a[i] b[i] a[i]-b[i]
0 0 0 0
1 0 0 0
2 1 1 0
3 1 1 0
4 2 2 0
5 4 4 0
6 7 7 0
7 13 13 0
8 24 24 0
9 44 44 0
10 81 81 0
11 149 149 0
12 274 274 0
13 504 504 0
14 927 927 0
15 1705 1705 0
16 3136 3136 0
17 5768 5768 0
18 10609 10609 0
19 19513 19513 0
20 35890 35890 0
21 66012 66012 0
22 121415 121415 0
23 223317 223317 0
24 410744 410744 0
25 755476 755476 0
26 1389537 1389537 0
27 2555757 2555757 0
28 4700770 4700770 0
29 8646064 8646064 0
30 15902591 15902591 0
31 29249425 29249425 0
32 53798080 53798080 0
33 98950096 98950096 0
34 181997601 181997601 0
35 334745777 334745777 0
36 615693474 615693474 0
37 1132436852 1132436852 0
38 2082876103 2082876103 0
39 3831006429 3831006429 0
40 7046319384 7046319384 0
41 12960201916 12960201916 0
42 23837527729 23837527729 0
43 43844049029 43844049029 0
44 80641778674 80641778674 0
45 148323355432 148323355432 0
46 272809183135 272809183135 0
47 501774317241 501774317241 0
48 922906855808 922906855808 0
49 1697490356184 1697490356184 0
50 3122171529233 3122171529233 0
51 5742568741225 5742568741225 0
52 10562230626642 10562230626642 0
53 19426970897100 19426970897100 0
54 35731770264967 35731770264967 0
55 65720971788709 65720971788709 0
56 120879712950776 120879712950775 1
57 222332455004452 222332455004451 1
58 408933139743937 408933139743935 2
59 752145307699165 752145307699160 5
60 1383410902447554 1383410902447546 8
61 2544489349890656 2544489349890641 15
62 4680045560037375 4680045560037345 30
63 8607945812375585 8607945812375530 55
64 15832480722303616 15832480722303512 104
65 29120472094716576 29120472094716384 192
66 53560898629395776 53560898629395416 360
67 98513851446415968 98513851446415296 672
68 181195222170528320 181195222170527072 1248
69 333269972246340096 333269972246337792 2304
70 612979045863284352 612979045863280000 4352
71 1127444240280152832 1127444240280144512 8320
72 2073693258389777408 2073693258389762048 15360
73 3814116544533214720 3814116544533186560 28160
74 7015254043203144704 7015254043203092480 52224
75 12903063846126137344 12903063846126039040 98304
76 23732434433862496256 23732434433862311936 184320
77 43650752323191775232 43650752323191439360 335872
78 80286250603180408832 80286250603179769856 638976
79 147669437360234692608 147669437360233480192 1212416
80 271606440286606852096 271606440286604623872 2228224
81 499562128250021937152 499562128250017808384 4128768
82 918838005896863416320 918838005896855945216 7471104
83 1690006574433492271104 1690006574433477853184 14417920
84 3108406708580377427968 3108406708580351213568 26214400
85 5717251288910732984320 5717251288910684749824 48234496
86 10515664571924601634816 10515664571924511457280 90177536
87 19341322569415712047104 19341322569415540080640 171966464
88 35574238430251047714816 35574238430250737336320 310378496
89 65431225571591361396736 65431225571590774194176 587202560
90 120346786571258121158656 120346786571257030639616 1090519040
91 221352250573100547047424 221352250573098500227072 2046820352
92 407130262715950037991424 407130262715946279895040 3758096384
93 748829299860308689420288 748829299860301844316160 6845104128
94 1377311813149359375122432 1377311813149346221785088 13153337344
95 2533271375725618236751872 2533271375725593540689920 24696061952
96 4659412488735286167076864 4659412488735240533049344 45634027520
97 8569995677610263510515712 8569995677610179758653440 83751862272
98 15762679542071167914344448 15762679542071011148038144 156766306304
99 28992087708416715444453376 28992087708416427681644544 287762808832
100 53324762928098150090539008 53324762928097265327276032 884763262976
このように、$i$ が小さいときは一致しています。$i$ が大きいときは浮動小数点数による誤差が出てきていますが、数値に比べて誤差が小さいので一般項は正しいと考えられます。
計算量
Pythonでは累乗に用いるpow関数は繰り返し二乗法で実装されているため、ある数の$N$乗は$O(\log N)$ で高速に求められます。そのため、$a[N]$を 求めるのに必要な計算量は $O(N)$ ですが、$b[N]$ は $O(\log N)$ で求めることができます。計算機上では、$a[i]$では値が正確に求められるが計算量は大きい、$b[i]$では高速に計算できるが誤差が発生する、というそれぞれのメリットとデメリットがあります。
トリボナッチ数列をさらに拡張
拡張フィボナッチ数列(k-ボナッチ数列 : 直前のk項の和)の一般項を線形代数とPythonで求めるという記事を書きました。こちらもよければご覧ください。
参考文献
こちらのpdfを参考にさせていただきました。ありがとうございました。