斐波拉契数列的特征值解法

2020-12-31 本文已影响0人 MatrixYe

斐波拉契数列是是指这样一个数列：0，1，1，2，3，5，8，13，每一项都等于前两项之和，其表达式为：
$a_{n}=\left\{\begin{array}{ll} n & n=0,1 \\ a_{n-1}+a_{n-2} & n \geqslant 2 \end{array}\right.$
直接递归时间复杂度很高，对于高次项计算耗时严重。如下列代码:

#直接迭代求法
def fb(n:int)->int:
    if n==0 or n==1:
        return n
    else:
        return fb(n-1)+fb(n-2)

# 求1～9 项
print([fb(i) for i in range(1,10)])

因此一般采用线性代数特征值解法来求解斐波拉契的通项公式。将相邻项作为二维向量， $x_1=(1,0),x_2=(1,1),x_3=(2,1)$ 可以推导出斐波拉契数列的线性代数版。
$\left[\begin{array}{l}a_2 \\a_1\end{array}\right]=\left[\begin{array}{l}1 &1\\1 & 0\end{array}\right]\left[\begin{array}{l}1 \\0\end{array}\right] \\\left[\begin{array}{l}a_3\\a_2\end{array}\right]=\left[\begin{array}{l}1 &1\\1 & 0\end{array}\right]\left[\begin{array}{l}1 &1\\1 & 0\end{array}\right]\left[\begin{array}{l}1 \\0\end{array}\right]\\ ......\\ \left[\begin{array}{l}a_n \\a_{n-1}\end{array}\right]=\left[\begin{array}{l}1 &1\\1 & 0\end{array}\right]^{n-1}\left[\begin{array}{l}1 \\0\end{array}\right]$
将上述通项简化即：
$x_n=A^{n-1}x_1$
其中 $x_1$ 为已知量 $\left[\begin{array}{l}1 \\0\end{array}\right]$

因此只要求解矩阵 $A$ 的n-1次幂即可获得通项式。但矩阵的计算时间复杂度更高，因此还需要简化，此时需要使用到矩阵的特征向量。矩阵的特征向量是指经过矩阵线性变化后，方向无变化，长度进行了压缩或拉升的向量。向量长度拉升的比例即为特征值 $\lambda$ 。
$A\vec{v}=\lambda \vec{v} \\$
求解矩阵A的特征值得到，特征值 $\lambda_1=1.618033988749895,\lambda_2=-0.6180339887498948$ ,对应的特征向量 $v_1=(0.85065081, -0.52573111),v_2=(0.52573111,0.85065081)$ ,注意这里的特征值是确定的，但特征向量实际上可以取任意一组。

因为向量 $x_1$ 与特征向量线性 $v_1,v_2$ 不相关，因此可以将 $x_1$ 替换为 $v_1$ 和 $v_2$ 的线性组合：
$x_1=c_1v_1+c_2v_2$
已知 $x_1,v_1,v_2$ 求线性方程组得到 $c_1=0.8506508083520401,c_2=-0.5257311121191336$ 。

此时，公式(3)可以表述为:
$x_n=A^{n-1}x_1=x_1=A^{n-1}(c_1v_1+c_2v_2)=c_1A^{n-1}v_1+c_2A^{n-1}v_2$
此时，方程的右侧 $v_1,v_2$ 是矩阵A的特征向量，经过n-1次A的线性映射后，其方向依旧不变，大小被连续缩放为 $\lambda$ 次。因此公式6可以写成:
$x_n=c_1A^{n-1}v_1+c_2A^{n-1}v_2=c_1\lambda_1^{n-1}v_1+c_2\lambda_2^{n-1}v_2$
此时， $c_1,c_2,v_1,v_2,\lambda_1,\lambda_2$ 都已经在上面求得，带入公式获得 $x_n$ 的通项公式。最后取 $x_n$ 的第一个分量，即为斐波拉契数列第n项 $a_n$

的表达式。由表达式可以看出，两个特征值的绝对值一个是大于1的，另一个是小于1的，因此随着n的递增， $\lambda_2$ 的影响会越来越弱，最终表达式会近似等于:
$x_n=c_1\lambda_1^{n-1}v_1$
而上面的表达式中， $\lambda_1=1.618033988749895$ ，即黄精分割。

最终代码实现如下：

"""
斐波拉契数列 的通项式
"""
import numpy as np
# X_n=A`n-1 x_1
# 向量 x_1
x1=np.array([1,0])
# 矩阵A
A=np.array([[1,1],[1,0]])
# 求矩阵a的特征值和特征向量
ks,vs=np.linalg.eig(A)
k_1=ks[0] # 特征值1
k_2=ks[1] # 特征值2
v_1=vs[0] #特征向量1
v_2=vs[1] # 特征向量2
print(f"求矩阵A的特征值为k_1={k_1},k_2={k_2}")
print(f"求矩阵A的特征向量为v_1={v_1},v_2={v_2}")
# 求解c1 和 c2,用v1和v2 的线性组合替换x_1, c1*v_1+c2*v_2=x1
c1,c2=np.linalg.solve(vs,x1)
print(f"求解得c1={c1},c2={c2}")
# 带入元公式 xn=A`(n-1)x1=A~(n-1) (c1v_1+c2v_2)=c1A`v_1+c2Av_2=c1k_1`v_1+c2k_2`v_2
# x_n的第一个分项就是a_n
fblq=lambda n: c1*k_1**(n-1)*v_1[0]+c2*k_2**(n-1)+v_2[0]
# 测试求前50项
print([fblq(i) for i in range(1,50)])

运行结果:

求矩阵A的特征值为k_1=1.618033988749895,k_2=-0.6180339887498948
求矩阵A的特征向量为v_1=[ 0.85065081 -0.52573111],v_2=[0.52573111 0.85065081]
求解得c1=0.8506508083520401,c2=-0.5257311121191336
[0.7236067977499789, 2.021471201601977, 2.2193468872328226, 3.7150869767156656, 5.408702751829356, 8.598058616425888, 13.481030256136112, 21.553357760442864, 34.50865690445984, 55.53628355278358, 89.51920934512428, 144.52976178578874, 233.52324001879393, 377.52727069246356, 610.5247795991384, 987.526319179483, 1597.525367666502, 2584.5259557338663, 4181.52559228825, 6765.525816909996, 10946.525678086125, 17711.525763884005, 28657.525710858015, 46368.5257436299, 75025.5257233758, 121393.52573589359, 196418.52572815728, 317811.5257329387, 514229.52572998387, 832040.5257318105, 1346269.5257306825, 2178309.5257313815, 3524578.5257309517, 5702887.525731222, 9227465.525731063, 14930352.525731174, 24157817.525731124, 39088169.525731176, 63245986.5257312, 102334155.52573125, 165580141.52573135, 267914296.52573153, 433494437.52573174, 701408733.5257323, 1134903170.5257328, 1836311903.525734, 2971215073.525736, 4807526976.52574, 7778742049.5257435]

斐波拉契数列的特征值解法

猜你喜欢

热点阅读