Skip to content
Lex Blog

用矩阵法与差分方程法推导斐波那契通项公式

本文给出两种推导斐波那契数列的方法：矩阵法与差分方程法。

在斐波那契的著作《计算之书》中，斐波那契数列定义如下：

Fn={0若 n=01若 n=1Fn−1+Fn−2若 n≥2F_n = \begin{cases} 0 & \text{若 } n = 0 \\ 1 & \text{若 } n = 1 \\ F_{n-1} + F_{n-2} & \text{若 } n \geq 2 \end{cases}

可以证明，其闭式公式为：

Fn=15((1+52)n−(1−52)n)F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)

以我目前的知识，只能理解下面两种证明方法：

矩阵法

首先讨论 n≥2n \geq 2 的情形。现在的目标是把斐波那契数列的递推公式化为矩阵形式。怎么做呢？我们可以从线性方程组的角度入手。

首先，已知：

Fn−1+Fn−2=FnF_{n-1} + F_{n-2} = F_{n}

我们可以加上方程 Fn−1+0⋅Fn−2=Fn−1F_{n-1} + 0 \cdot F_{n-2} = F_{n-1}，构成一个线性方程组：

{Fn−1+Fn−2=FnFn−1+0⋅Fn−2=Fn−1\begin{cases} F_{n-1} + F_{n-2} = F_{n} \\ F_{n-1} + 0 \cdot F_{n-2} = F_{n-1} \end{cases}

这可以化为如下的矩阵形式：

[1110][Fn−1Fn−2]=[FnFn−1]\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix} = \begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix}

现在，我们可以反复迭代这一过程：

[FnFn−1]=[1110][Fn−1Fn−2]=[1110]2[Fn−2Fn−3]=⋯=[1110]n−1[F1F0]\begin{align} \begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix} &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix} \\ &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^2 \begin{bmatrix} F_{n-2} \\ F_{n-3} \end{bmatrix} \\ &=\cdots \\ &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-1} \begin{bmatrix} F_{1} \\ F_{0} \end{bmatrix} \end{align}

记矩阵 A=[1110]\boldsymbol{A} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}。于是问题转化为求 An−1\boldsymbol{A}^{n-1}；之后我们还可以计算 An\boldsymbol{A}^{n}，并把所有 nn 换成 n−1n-1。

注意到矩阵 A\boldsymbol{A} 是方阵，我们可以利用矩阵的特征值与特征向量。

特征向量可以理解为：右乘该矩阵后，得到与自身平行的向量。特征值则是右乘矩阵时特征向量被缩放的比例因子。换言之：

Ax=λx(1)\boldsymbol{A}\boldsymbol{x} = \lambda\boldsymbol{x}\tag{1}

这里 λ\lambda 是特征值，非零向量 x∈Rn\boldsymbol{x} \in \mathbb{R}^n（其中 nn 为方阵的阶数）是矩阵 A\boldsymbol{A} 对应于特征值 λ\lambda 的特征向量。

因此具体做法是：先求出特征值 λ1\lambda_1 与 λ2\lambda_2（一般 nn 阶矩阵有 nn 个特征值），得到对角矩阵 diag{λ1,λ2}\rm{diag}\{\lambda_1, \lambda_2\}；然后寻找可逆矩阵 P\boldsymbol{P}，使得：

P−1AP=diag{λ1,λ2}\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P} = \rm{diag}\{\lambda_1, \lambda_2\}

利用矩阵乘法的性质：

(P−1AP)n=P−1A(PP−1)A(P⋯P−1)AP=P−1AnP(2)(\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P})^n = \boldsymbol{P}^{-1}\boldsymbol{A}(\boldsymbol{P}\boldsymbol{P}^{-1})\boldsymbol{A}(\boldsymbol{P}\cdots\boldsymbol{P}^{-1})\boldsymbol{A}\boldsymbol{P} = \boldsymbol{P}^{-1}\boldsymbol{A}^n\boldsymbol{P}\tag{2}

这样就可以算出 An\boldsymbol{A}^{n}。

先求特征值。我们可以把方程 (1)(1) 改写成：

(A−λE)x=0(3)(\boldsymbol{A}-\lambda\boldsymbol{E})\boldsymbol{x} = \boldsymbol{0}\tag{3}

这里 E\boldsymbol{E} 是单位矩阵，可以算出 A−λE=[1−λ11−λ]\boldsymbol{A}-\lambda\boldsymbol{E} = \begin{bmatrix} 1-\lambda & 1 \\ 1 & -\lambda \end{bmatrix}。为保证存在非零解，我们求解：

∣A−λE∣=∣1−λ11−λ∣=λ2−λ−1=0\left| \boldsymbol{A}-\lambda\boldsymbol{E} \right| = \begin{vmatrix} 1-\lambda & 1 \\ 1 & -\lambda \end{vmatrix} = \lambda^2 - \lambda - 1 = 0

解此方程，得到 λ1=1+52\lambda_1 = \frac{1+\sqrt{5}}{2} 与 λ2=1−52\lambda_2 = \frac{1-\sqrt{5}}{2}。

因此，对角矩阵为：

diag{λ1,λ2}=[1+52001−52]\rm{diag}\{\lambda_1, \lambda_2\} = \begin{bmatrix} \frac{1+\sqrt{5}}{2} & 0 \\ 0 & \frac{1-\sqrt{5}}{2} \end{bmatrix}

假设特征向量为 x=[xy]T\boldsymbol{x} = \begin{bmatrix} x & y \end{bmatrix}^T，将 λ1\lambda_1 与 λ2\lambda_2 分别代入方程 (2)(2)，得到两个方程组：

{(1−λ1)x+y=0x−λ1y=0,{(1−λ2)x+y=0x−λ2y=0\begin{cases} (1-\lambda_1)x + y = 0 \\ x - \lambda_1y = 0 \end{cases}, \begin{cases} (1-\lambda_2)x + y = 0 \\ x - \lambda_2y = 0 \end{cases}

在两个方程组中都令 y=1y = 1，即可得到矩阵的两个特征向量：

x1=[1+521]T,x2=[1−521]T\boldsymbol{x}_1 = \begin{bmatrix} \frac{1+\sqrt{5}}{2} & 1 \end{bmatrix}^T, \boldsymbol{x}_2 = \begin{bmatrix} \frac{1-\sqrt{5}}{2} & 1 \end{bmatrix}^T

于是，可逆矩阵 P\boldsymbol{P} 由这两个特征向量构成：

P=[1+521−5211]\boldsymbol{P} = \begin{bmatrix} \frac{1+\sqrt{5}}{2} & \frac{1-\sqrt{5}}{2} \\ 1 & 1 \end{bmatrix}

为什么是这样呢？设 P=[x1x2y1y2]\boldsymbol{P} = \begin{bmatrix} x_1 & x_2 \\ y_1 & y_2 \end{bmatrix}（其中 x1,y1,x2,y2x_1, y_1, x_2, y_2 分别是特征向量 x1,x2\boldsymbol{x}_1, \boldsymbol{x}_2 的分量）。

现在，若计算 Pdiag{λ1,λ2}\boldsymbol{P}\rm{diag}\{\lambda_1, \lambda_2\}，它恰好等于 [λ1x1λ2x2λ1y1λ2y2]\begin{bmatrix} \lambda_1x_1 & \lambda_2x_2 \\ \lambda_1y_1 & \lambda_2y_2 \end{bmatrix}。这说明 AP=Pdiag{λ1,λ2}\boldsymbol{A}\boldsymbol{P} = \boldsymbol{P}\rm{diag}\{\lambda_1, \lambda_2\} 是必然成立的。两边左乘 P−1\boldsymbol{P}^{-1}，得到 P−1AP=diag{λ1,λ2}\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P} = \rm{diag}\{\lambda_1, \lambda_2\}。因此，P=[x1x2y1y2]\boldsymbol{P} = \begin{bmatrix} x_1 & x_2 \\ y_1 & y_2 \end{bmatrix} 是合理的。

其逆矩阵很容易计算：

P−1=P∗∣P∣=15[11−5211+52]\boldsymbol{P}^{-1} = \frac{\boldsymbol{P}^*}{\left|\boldsymbol{P}\right|} = \frac{1}{\sqrt{5}}\begin{bmatrix} 1 & \frac{1-\sqrt{5}}{2} \\ 1 & \frac{1+\sqrt{5}}{2} \end{bmatrix}

代入方程 (2)(2)：

An=P(P−1AP)nP−1=Pdiagn{λ1,λ2}P−1\boldsymbol{A}^n = \boldsymbol{P}(\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P})^n\boldsymbol{P}^{-1} = \boldsymbol{P}\rm{diag}^n\{\lambda_1, \lambda_2\}\boldsymbol{P}^{-1}

其中对角矩阵为：

diagn{λ1,λ2}=[(1+52)n00(1−52)n]\rm{diag}^n\{\lambda_1, \lambda_2\} = \begin{bmatrix} \left(\frac{1+\sqrt{5}}{2}\right)^n & 0 \\ 0 & \left(\frac{1-\sqrt{5}}{2}\right)^n \end{bmatrix}

因此，

An=15[(1+52)n+1+(1−52)n+1(1+52)n+1(1−52)+(1−52)n+1(1+52)(1+52)n+(1−52)n(1+52)n(1−52)+(1−52)n(1+52)]\boldsymbol{A}^n = \frac{1}{\sqrt{5}}\begin{bmatrix} \left(\frac{1+\sqrt{5}}{2}\right)^{n+1} + \left(\frac{1-\sqrt{5}}{2}\right)^{n+1} & \left(\frac{1+\sqrt{5}}{2}\right)^{n+1}\left(\frac{1-\sqrt{5}}{2}\right) + \left(\frac{1-\sqrt{5}}{2}\right)^{n+1}\left(\frac{1+\sqrt{5}}{2}\right) \\ \left(\frac{1+\sqrt{5}}{2}\right)^n + \left(\frac{1-\sqrt{5}}{2}\right)^n & \left(\frac{1+\sqrt{5}}{2}\right)^{n}\left(\frac{1-\sqrt{5}}{2}\right) + \left(\frac{1-\sqrt{5}}{2}\right)^{n}\left(\frac{1+\sqrt{5}}{2}\right) \end{bmatrix}

于是，

[FnFn−1]=An−1[F1F0]=15[(1+52)n+(1−52)n(1+52)n−1+(1−52)n−1]\begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix} = \boldsymbol{A}^{n-1} \begin{bmatrix} F_1 \\ F_0 \end{bmatrix} = \frac{1}{\sqrt{5}}\begin{bmatrix} \left(\frac{1+\sqrt{5}}{2}\right)^{n} + \left(\frac{1-\sqrt{5}}{2}\right)^{n} \\ \left(\frac{1+\sqrt{5}}{2}\right)^{n-1} + \left(\frac{1-\sqrt{5}}{2}\right)^{n-1} \end{bmatrix}

只考虑方程两边矩阵的第一行，即可得到斐波那契数列的闭式公式：

Fn=15((1+52)n−(1−52)n)F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)

代入 n=0n=0 与 n=1n=1 验证，发现两者都满足该方程。

差分方程法

把数列 {an}\{a_n\} 的一阶差分定义为 Δan=an+1−an\Delta a_n = a_{n+1} - a_n（这里采用向后差分），则可以定义二阶差分：

Δ2an=Δan+1−Δan=an+2−2an+1+an\Delta^2 a_n = \Delta a_{n+1} - \Delta a_n = a_{n+2} - 2a_{n+1} + a_n

进一步，可以定义 mm 阶差分：

Δman=Δm−1an+1−Δm−1an=∑i=0m(−1)iCmian+m−i\Delta^m a_n = \Delta^{m-1} a_{n+1} - \Delta^{m-1} a_n = \sum_{i=0}^{m} (-1)^i C_m^i a_{n+m-i}

以及

F(an,Δan,Δan,Δ2an,Δ3an,…)F(a_n, \Delta a_n, \Delta a_n, \Delta^2 a_n, \Delta^3 a_n, \ldots)

定义如 Fn=Fn−1+Fn−2F_n = F_{n-1} + F_{n-2} 的斐波那契数列，其二阶差分是常数：

Δ2Fn=Fn+2−2Fn+1+Fn=Fn+1+Fn+1−2Fn+1+Fn=Fn+1−Fn=Fn\Delta^2 F_n = F_{n+2} - 2F_{n+1} + F_n = F_{n+1} + F_{n+1} - 2F_{n+1} + F_n = F_{n+1} - F_n = F_n

现在，我们可以写出 FnF_n 的二阶线性齐次差分方程：

Δ2Fn−Fn=0\Delta^2 F_n - F_n = 0

这是一个特征方程，我们可以假设 Fn=rnF_n = r^n 来求解：

r2−1=0r^2 - 1 = 0

解此方程得到两个解：r=1r = 1 与 r=−1r = -1。因此，齐次差分方程的通解为：

Fn=c1⋅1n+c2⋅(−1)nF_n = c_1 \cdot 1^n + c_2 \cdot (-1)^n

现在，我们需要初始条件来求特解。已知 F0=0F_0 = 0、F1=1F_1 = 1。将它们代入通解：

F0=c1⋅10+c2⋅(−1)0=c1+c2=0F1=c1⋅11+c2⋅(−1)1=c1−c2=1\begin{align*} F_0 &= c_1 \cdot 1^0 + c_2 \cdot (-1)^0 = c_1 + c_2 = 0 \\ F_1 &= c_1 \cdot 1^1 + c_2 \cdot (-1)^1 = c_1 - c_2 = 1 \end{align*}

解这个方程组，得 c1=12c_1 = \frac{1}{2}、c2=−12c_2 = -\frac{1}{2}。因此，斐波那契数列的特解为：

Fn=12⋅1n−12⋅(−1)nF_n = \frac{1}{2} \cdot 1^n - \frac{1}{2} \cdot (-1)^n

利用 (−1)n=(−1)n−1⋅(−1)=−(−1)n−1(-1)^n = (-1)^{n-1} \cdot (-1) = -(-1)^{n-1}，可以进一步化简：

Fn=12−12⋅(−1)n−1F_n = \frac{1}{2} - \frac{1}{2} \cdot (-1)^{n-1}

这确实是斐波那契数列的闭式公式：

Fn=15((1+52)n−(1−52)n)F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)

这样，我们便用差分方程法成功地推导出了相同的结果。

总之，矩阵法与差分方程法都得出了斐波那契数列相同的闭式表达式，展现了数学以多种途径殊途同归之美。