フィボナッチ数列

授業担当:劉雪峰

1. 数列の定義と生成

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

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

フィボナッチ数列の最初の数字を書いてみます。 $$ 1,1,2,3,5,8,13, \cdots $$

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

フィボナッチ数は以下の兎の問題で考えることが多いです。
  • 1つがいの兎は、産まれて2か月後から毎月1つがいの兎を産む。
  • 兎が死ぬことはない。
  • この条件のもとで、産まれたばかりの1つがいの兎は1年の間に何つがいの兎になるか?
毎月のつがいの数(F(n))について、前月のつがい数(これまでにあるもの=F(n-1))+2ヶ月前の兎が新しく産む兎(2ヶ月前の兎のつがい数=F(n-2))となっています。即ち、$F(n)=F(n-1)+F(n-2)$です。よって、フィボナッチ数列が現れていることがわかります。

フィボナッチ数列を作成するプログラム

フィボナッチ数列の第10項まで作成してみます。

  • F=[1,1] は成分2つのある配列を作成します。配列の各成分は1です。

  • 配列の最初の値と2番目の値はそれぞれ $F[0], F[1]$ で現しています。

In [1]:
#演習1
import matplotlib.pyplot as plt

N=10
n_list=range(1,N+1)
F=[1,1]
print("第%d月のつがい数は%dです。"%(1, F[0]))
print("第%d月のつがい数は%dです。"%(2, F[1]))

# 第3月からのつがい数(F[2])はforループで作成します。
for n in range(2,N):
    F.append(F[n-1]+F[n-2])
    print("第%d月のつがい数は%dです。"%(n+1, F[n]))
    
plt.plot(n_list,F,'r.-')
plt.grid()
plt.show()
第1月のつがい数は1です。
第2月のつがい数は1です。
第3月のつがい数は2です。
第4月のつがい数は3です。
第5月のつがい数は5です。
第6月のつがい数は8です。
第7月のつがい数は13です。
第8月のつがい数は21です。
第9月のつがい数は34です。
第10月のつがい数は55です。
<matplotlib.figure.Figure at 0x7f60d84ffe10>

演習1

$N=10,15,20,25$に対して数列を考察しなさい。

  • $N$までのフィボナッチ数列$F$を作りなさい。
  • 数列$F$のグラフを描いてみます。

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

  • 数列は発散しますか?
  • フィボナッチ数列のグラフはどちらの関数のグラフに似ていますか? 

$$f(x)=x^a, f(x)=a^x, f(x)=\log(x)$$

In [10]:
#ここにコードを書いてください。

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

以下のコードでは、数列Fのlog値を計算して、グラフを描いています。$\log(F)$のグラフはどちらの関数に近いでしょうか?

In [2]:
#演習2
import matplotlib.pyplot as plt
import math

N=10
F=[1,1]
n_list = [1,2]
for n in range(2,N):
    F.append(F[n-1]+F[n-2])
    n_list.append(n+1)

log_F=[]
for n in range(N):
    log_F.append( math.log( F[n] ) )
plt.plot(n_list, log_F, marker='o')
plt.grid()
plt.show()

演習 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 \times \exp(b) $$

  • $c:=exp(b), r:=exp(a)$とおく。$N$の値が大きくなる時、$c$と$r$の値を計算してください。$r$の値と黄金比との関係を検討しなさい。

$$ \mbox{黄金比}=\frac{\sqrt{5}+1}{2} $$

補足1:多くの点に当てはまる直線 y=ax+b の係数 a と bの求める方法

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

例えば、「a=0.45,b=-0.40」に対応する直線のグラフを描いてみます。 係数a,bの値を調整しながら、もっと良いa,bの値を推定してください。

In [3]:
#補足1
x_list=range(1,N+1)
a=0.45; b=-0.40;
y_list =[]
for x in x_list:
     y_list.append(a*x + b)
plt.scatter(n_list, log_F, marker='o')
plt.plot(x_list, y_list,"r-")
plt.grid()

補足1:2点に通っている直線 y=ax+b の係数 a と bの求める方法

$(x_1,y_1), (x_2,y_2)$を通っている直線を求めます。

$$ y = y_1 + (x-x_1)\frac{y_2-y_1}{x_2-x_1} = \left(\frac{y_2-y_1}{x_2-x_1} \right) x + \left( y_1 -x_1 \frac{y_2-y_1}{x_2-x_1} \right) $$

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

In [4]:
#補足1
x1=2; x2=7;
y1=math.log(F[1]); y2=math.log(F[6]);
a=(y2-y1)/(x2-x1)
b= y1 - x1*a

y_list =[]
x_list=[x1,x2]
for x in x_list:
     y_list.append(a*x + b)
     
plt.plot(n_list, log_F, 'o')
plt.plot(x_list, y_list, "r.-")
plt.grid()
plt.show()

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}) $$

演習 3 (オプション)

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)$のグラフに当てはまる直線を求めなさい。

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

In [6]:
#演習3
import matplotlib.pyplot as plt

import numpy as np

A=np.ones((4,2)) #全ての成分が1となる 4*2の行列を作成する。
y=np.zeros((4,1)) #全ての成分が0となる 4*1の行列を作成する。

A[:,0]=np.array([1,2,3,4])
y[:,0]=np.array([4.444, 5.079, 6.534, 8.572]);

x_list=A[:,0] #行列の第1列を取って、xにする。
y_list=y[:,0]
plt.plot(x_list,y_list,'o')

alpha=np.linalg.inv(A.T@A)@(A.T@y) # A.T@y : 行列Aの転置とベクトルyの積; @:行列の積計算の演算子

a=alpha[0,0]; b=alpha[1,0];
print("a=", a, " b=",b)

new_y_list=a*x_list+b #Numpyの配列は一体として四則演算などの数学関数の計算ができます。

plt.plot(x_list,new_y_list,'r-.')
plt.grid()
plt.show()
a= 1.3839  b= 2.6975

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. コーヒーブレーク

フィボナッチ数列によるゴールデン・スパイラルの図を描く。Nを他の値にしてからグラフを描いてみよう。 <!-- 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]) -->

In [25]:
import math
N=7
N_list=range(1,N+1)
F=[]
x_indx=[]
y_indx=[]
for k in range(0,N+1):
    if(k<2):
        F.append(1) 
    else:
        F.append(F[k-1] +F[k-2]) 
    y_indx.append(pow(-1,int(k/2)))
    x_indx.append(pow(-1,int((k+1)/2)))
    
import matplotlib.patches as pat
fig = plt.figure()
ax = fig.add_subplot(111)    

x=[0];y=[0];
#text( -1,0.3, '1');
for k in range(0,N):
    x.append(x[k]+ x_indx[k]*F[k+1])
    y.append(y[k]+ y_indx[k]*F[k+1])
    s1=90*(k-1)
    x0=x[k]-F[k+1]*math.cos(s1/180*math.pi)
    y0=y[k]-F[k+1]*math.sin(s1/180*math.pi)
    c=pat.Arc( xy= (x0, y0), width=2*F[k+1], height=2*F[k+1],theta1=s1, theta2=s1+90)
    ax.add_patch(c)
    x_list = [x[k],x[k+1],x[k+1],x[k]]; #x座標のリスト
    y_list = [y[k],y[k],y[k+1],y[k+1]]; #y座標のリスト
    plt.fill(x_list,y_list,alpha=0.5)
    plt.text((x[k]+x[k+1])/2, (y[k]+y[k+1])/2,'%d'%F[k+1])
#plt.plot(x,y)

plt.grid()
plt.show()

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

The original matlab code is revised to python language by Wada Kaoru. 2020/04

計算機演習A・Bのノート 授業担当:劉雪峰