フィボナッチ数列

1. 数列の定義と生成

フィボナッチ数列(Fibonacci sequence)は以下のように、再帰的に定義されています。

  • $F_1 = 1$
  • $F_2=1$
  • $F_n = F_{n-1} + F_{n-2}\quad \quad (n>2)$

フィボナッチ数列と兎の問題

フィボナッチ数は以下の兎の問題で始まってを考えることが多いです。

  • 1つがいの兎は、産まれて2か月後から毎月1つがいずつの兎を産む。
  • 兎が死ぬことはない。
  • この条件のもとで、産まれたばかりの1つがいの兎は1年の間に何つがいの兎になるか?
つがいの数について、どの月のつがいの合計も、その前の2つの月での合計の和となり、フィボナッチ数列が現れていることがわかります。

演習1

  • プログラで$N$までのフィボナッチ数を作りなさい。。
  • 数列$F$のグラフを描いてみます。
  • $N=10,15,20,25$に対して数列を考察しなさい。

以下のことを考えてください。

  • 数列は発散しますか?
  • グラフはどちらの関数のグラフに似ってますか? $f(x)=x^a$, $f(x)=a^x$, $f(x)=\log(x)$ ?

プログラミングの準備

  • F=[1,1] は長さが2である配列を作成します。配列の成分は1、1です。
In [2]:
N=10;
F=[1,1];
for k=3:N
   F(k) = F(k-1) + F(k-2);
end
%printf("%15d \n",F)
F;
plot(F,'-+')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 10 20 30 40 50 60 0 2 4 6 8 10 gnuplot_plot_1a

フィボナッチ数列のグラフを描いてみます。

2: フィボナッチ数列の式を推定する

数列Fのlogを計算してグラフで考察してください。$\log(F)$のグラフはどちらの関数に近いでしょうか?

演習 2

  • $\log(F)$に当てはまる直線(フィッティング直線)を求めてください。即ち、$\log(F)$のグラフに近い直線$y=ax+b$の係数$a$と$b$を求めてください。解く方法の使用は自由です。
$$ \log(F_n)\approx a\times n+b $$
  • $n$ ~ $\log(F_n)$ のグラウは直線に近いなので,$F_n$を以下のように近似できます。
$$ \log(F_n)\approx a \times n+b \Rightarrow F_n \approx \exp(a)^n*\exp(b) $$
  • $c:=exp(b), r:=exp(a)$とおく。$N$の値が大きくなる時、$c$と$r$の値を計算してください。$r$の値と黄金比との関係を検討しなさい。

演習2の結果と結論

(ここに書いてください)

In [9]:
plot(1:N, log(F),'-o')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 1 2 3 4 5 0 2 4 6 8 10 gnuplot_plot_1a

補足1:直線 y=ax+b の係数 a と bの求める方法1

良さそうなa,bに対応する直線のグラフ書いて、それと$n-log(F_n)$のグラフと比較してから、より良い近似になる係数の値を見つけます。

例えば、「a=0.45,b=-0.7」に対応する直線のグラフを描いてみます。

In [12]:
plot(log(F),'.')
x=1:N;
a=0.45;b=-0.7;

my_log_F = a*x + b;
hold on 
plot(my_log_F)
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -1 0 1 2 3 4 5 0 2 4 6 8 10 gnuplot_plot_1a gnuplot_plot_2a

補足2:直線 y=ax+b の係数 a と bの求める方法2

(x1,y1), (x2,y2)を通っている直線を求めます。

$$ y = y1 + (x-x1)\frac{y2-y1}{x2-x1} $$

よって、2点の座標を使って、$y=ax+b$の係数a,bを求めることができます。

In [13]:
x1=3,x2=10;
y1=log(F(3));y2=log(F(10));
a=(y2-y1)/(x2-x1)
b= y1 - x1*a

plot(log(F),'.')

my_log_F = a*(1:N)+b;
hold on 
plot(my_log_F)
grid on
r=exp(a)
x1 =  3
a =  0.47346
b = -0.72722
r =  1.6055
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -1 0 1 2 3 4 5 0 2 4 6 8 10 gnuplot_plot_1a gnuplot_plot_2a

3. [オプション] 最小二乗法による当てはまる直線の計算方法

最小二乗法とは

与えられる点$(x_i,y_i)_{i=1}^N$について、一次関数$y=ax+b$を用い、点列に対してよい近似となるように、残差の二乗和を最小とするような係数を決定する方法です。

定式化

数式で考えるとき、以下の式を最小化する$a,b$を求める問題となります。

$$ V(a,b)=\sum_{i=1}^N (ax_i + b - y_i)^2 $$

解法

行列を利用して、上記の最小化問題を考えます。

$$ \mathbf{x}=\left(\begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_N \end{array}\right),\quad \mathbf{y}=\left(\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array}\right) $$

上記の$\mathbf{x},\mathbf{y}$を利用して、 $$ V(a,b)= (a \mathbf{x} + b - \mathbf{y})^T\cdot ( a \mathbf{x} + b - \mathbf{y}) $$ さらに、以下の$A$を定義します。 $$ \mathbf{A}=\left(\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_N & 1 \end{array}\right),\quad \alpha = \left(\begin{array}{c} a \\ b \end{array}\right) $$ よって, $a \mathbf{x} + b = A\alpha $,

$$ V(\alpha)=(A \alpha - \mathbf{y})^T\cdot ( A \alpha- \mathbf{y}) $$

2変数関数$V(\alpha)$が最小値を取るとき、$\alpha=(a,b)$に関する微分が$0$になることがあります。$V$の$\alpha=(a,b)$に関する微分式を整理すると、以下の式になります。

$$ (A^TA) \alpha - A^T \mathbf{y} = 0 $$

$A^TA$が正定値であれば、$\alpha$を次のように算出できます。

$$ \alpha = (A^TA)^{-1} (A^T\mathbf{y}) $$

演習 4

1) 以下の点列に対して、点列に当てはまる直線$y=ax+b$を求めなさい。

  • x=1, y=4.44357
  • x=2, y=5.07900
  • x=3, y=6.53377
  • x=4, y=8.57194

2) $n$~$\log(F_n)$のグラフに当てはまる直線を求めなさい。

以下は演習4.1の解答例です。

In [8]:
x=[ 1 2 3 4]';
y=[4.444, 5.079, 6.534, 8.572]';
plot(x,y,'o')
grid on

%% Aを作って、以下のコードを使って alpha=(a,b)を計算しなさい。
A=[x, ones(4,1)]

alpha=(A'*A)\(A'*y)

a=alpha(1)
b=alpha(2)
my_y = a.*x + b;

hold on
plot(x,my_y,'r-+')
A =

   1   1
   2   1
   3   1
   4   1

alpha =

   1.3839
   2.6975

a =  1.3839
b =  2.6975
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 4 5 6 7 8 9 1 1.5 2 2.5 3 3.5 4 gnuplot_plot_1a gnuplot_plot_2a
In [7]:
(sqrt(5)+1)/2
ans =  1.6180

4. [オプション] フィボナッチ数列の式を求める 

$F_n$の式を推定(1)

$F_n = cr^n$を仮定して、この式を$F_n=F_{n-1}+F_{n-2}$に代入して、$r$に関する方程式が分かります。

$$ r^2=r+1 $$

よって、 $$ r=\frac{1+\sqrt{5} }{2}, \quad r=\frac{1-\sqrt{5} }{2} $$

数列$F_n$は無限に発散するので、$r$を正の値を取ります。よって、 $$ F_n = c \left[\frac{1+\sqrt{5} }{2}\right]^n $$ 上記の式はフィボナッチ数列にならないですが、フィボナッチ数列の主要項となります。

$F_n$の式を推定(2)

$z_n:=F_n - c \left[ \frac{1+\sqrt{5} }{2}\right]^n$を考察します。

$Z_n$は次ぎの式を満たす。

$$ Z_n =Z_{n-1} + Z_{n-2} $$

よって、 $Z_n$も以下の形式を持っていると推定します。 $$ Z_n = \tilde{c} \tilde{r}^n $$ この式を満たす$\tilde{r}$は$\tilde{r}^2=\tilde{r}+1$の解となりますので、 $$ \tilde{r} = \frac{1-\sqrt{5} }{2} $$ よって、$F_n$に関して、新しい式が考えられます。 $$ F_n = c r^n + \tilde{c} \tilde{r}^n $$ 条件$F_1=1,F_2=1$を利用して、以下のことが分かります。

$$ c = \frac{1}{\sqrt{5}}, \quad \tilde{c} = - \frac{1}{\sqrt{5}} $$

即ち、

$$ F_n = \frac{1}{\sqrt{5}} ( \left[ \frac{1+\sqrt{5} }{2} \right] ^n - \left[\frac{1-\sqrt{5} }{2}\right]^n ) $$

5. コーヒーブレーク

フィボナッチ数列によるゴールデン・スパイラルの図を描く。

In [28]:
N=8;
F=[1,1];
for k=3:N
   F(k) = F(k-1) +F(k-2);
end

x_indx=[-(-1).^(1:(N/2));(-1).^(1:(N/2))];
y_indx=[-(-1).^(1:(N/2));-(-1).^(1:(N/2))];

x=[0];y=[0];

hold on
text( -1,0.3, '1');
for k=1:N-1
  x(k+1)=x(k) + x_indx(k)*F(k+1);
  y(k+1)=y(k) + y_indx(k)*F(k+1);
  x1 = min(x(k), x(k+1));  x2 = max(x(k), x(k+1));
  y1 = min(y(k), y(k+1));  y2 = max(y(k), y(k+1));
  x_list = [x1,x2,x2,x1,x1];
  y_list = [y1,y1,y2,y2,y1];
  fill(x_list(1:4), y_list(1:4),k-1);
  plot(x_list, y_list,'-');
  plot([x(k),x(k+1)],[y(k),y(k+1)],'-');
  text( (x1+x2)/2, (y1+y2)/2, num2str(F(k+1)));
end
daspect ([1 1])
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -5 0 5 10 15 20 -25 -20 -15 -10 -5 0 5 10 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a gnuplot_plot_6a gnuplot_plot_7a gnuplot_plot_8a gnuplot_plot_9a gnuplot_plot_10a gnuplot_plot_11a gnuplot_plot_12a gnuplot_plot_13a gnuplot_plot_14a gnuplot_plot_15a gnuplot_plot_16a gnuplot_plot_17a gnuplot_plot_18a gnuplot_plot_19a gnuplot_plot_20a gnuplot_plot_21a gnuplot_plot_22a gnuplot_plot_23a gnuplot_plot_24a gnuplot_plot_25a gnuplot_plot_26a gnuplot_plot_27a gnuplot_plot_28a 1 1 2 3 5 8 13 21

Created by Xuefeng LIU, 2017/04/23, modified on 2018/04/22.

In [ ]: