求解线性递推方程第n项的一般方法

概述

系数为常数，递推项系数均为一次的，形如下面形式的递推式，称为线性递推方程。

\[f[n]=\begin{cases} C &n\in Value\\ a_1 f[n-1]+a_2 f[n-2]+⋯a_t f[n-t]&n∉Value \end{cases}\]

\((a_1,a_2,…,a_t,C∈\mathbb{R},0<t<n)\)
其中\(Value\)为终止条件的集合。

例如：斐波那契\((Fibonacci)\)数列则通过下面这个线性递推方程定义

\[f[n]=\begin{cases} 1 &n=0\ \mathrm{or}\ n=1\\ f[n-1]+f[n-2] &n\geq2 \end{cases}\]

其中\(n\in \mathbb{N}\).

直接求解法

从小到大从终止条件正推至\(x=n\)项，或利用递归从第\(n\)项开始反推到递推终止条件。
时间复杂度\(O(n)\).

待定系数法（常用于递推项为两项的）

一般对于递推项仅有两项的，有递推关系式

\[f[n]=\begin{cases} C &n∈Value\\ a_1 f[n-1]+a_2 f[n-2] &n∉Value \end{cases}\]

假设系数\(x,y\)满足\(f[n]+yf[n-1]=x(f[n-1]+yf[n-2])\)，整理得\(\begin{cases}x-y=a_1\\xy=a_2\end{cases}\)，可以解出系数\(x,y\)。
根据等比数列性质

\[a_n=a_1\cdot q^{n-1} \]

得到\(f[n]\)和\(f[t](t∈Value)\)的关系式，得到时间复杂度为\(O(1)\)的通解。

例如：对于斐波那契\((Fibonacci)\)数列，我们得到方程组

\[\begin{cases}x-y=1\\xy=1\end{cases} \]

解得

\[\begin{cases}x_1=\frac{\sqrt{5}+1}{2}\\y_1=\frac{\sqrt{5}-1}{2}\end{cases},\begin{cases}x_2=\frac{1-\sqrt{5}}{2}\\y_2=-\frac{1+\sqrt{5}}{2}\end{cases} \]

我们得到

\[\begin{cases}f[n]+y_1f[n-1]=x_1^{n-1}(f[1]+y_1f[0])=x_1^n\\f[n]+y_2f[n-1]=x_2^{n-1}(f[1]+y_2f[0])=x_2^n\end{cases} \]

联立两式，将\(f[n-1]\)消去，可以得到

\[f[n]=\cfrac{x_2^ny_1-x_1^ny_2}{y_1-y_2}=\cfrac{\left(\frac{1+\sqrt5}{2}\right)^{n+1}-\left(\frac{1-\sqrt5}{2}\right)^{n+1}}{\sqrt5} \]

即为\(f[n]\)的通解。

矩阵乘法

对于递推项数较多，规模较大的递推式，我们一般采用矩阵运算的方式。
先将递推式写成一个如下式所呈现的\(t\)行方程组：

\[\begin{cases} f[n]&=a_1&f[n-1]&+a_2&f[n-2]&+⋯+a_{t-1}&f[n-t+1]&+a_t&f[n-t]\\ f[n-1]&=1&f[n-1]&+0&f[n-2]&+⋯+0&f[n-t+1]&+0&f[n-t]\\ f[n-2]&=0&f[n-1]&+1&f[n-2]&+⋯+0&f[n-t+1]&+0&f[n-t]\\ ⋮&&⋮\\ f[n-t+1]&=0&f[n-1]&+0&f[n-2]&+⋯+1&f[n-t+1]&+0&f[n-t]\\\end{cases}\]

写作\(t×t\)的矩阵与列向量的乘积形式

\[\left[\begin{matrix}f[n]\\f[n-1]\\f[n-2]\\⋮\\f[n-t+1]\end{matrix}\right]=\left[\begin{matrix}a_1&a_2&⋯&a_{t-1}&a_t\\1&0&⋯&0&0\\0&1&⋯&0&0\\⋮&⋮&⋱&⋮&⋮\\0&0&⋯&1&0\end{matrix}\right]\left[\begin{matrix}f[n-1]\\f[n-2]\\f[n-3]\\⋮\\f[n-t]\end{matrix}\right]\]

令\(A\)为中间这个\(t\)阶方阵，则有

\[\left[\begin{matrix}f[n]\\f[n-1]\\f[n-2]\\⋮\\f[n-t+1]\end{matrix}\right]=A^{k}\left[\begin{matrix}f[n-k]\\f[n-k-1]\\f[n-k-2]\\⋮\\f[n-k-t+1]\end{matrix}\right] \]

其中\(n-k-t+1>0\)且\(n-k∈Value\)。

例如：对于斐波那契\((Fibonacci)\)数列，我们可以将其写作

\[\left[\begin{matrix}f[n]\\f[n-1]\end{matrix}\right]=\left[\begin{matrix}1&1\\1&0\end{matrix}\right]^{n-1}\left[\begin{matrix}f[1]\\f[0]\end{matrix}\right] \]

对于如何求解\(A^n\)（矩阵的幂），一般有如下方法：

朴素的矩阵连乘法

时间复杂度为\(O(nt^3)\).
常规的写法为

void multiply()
{
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            for (int k = 0; k < n; ++k)
                t[i][j] += x[i][k] * y[k][j];
}

稍加优化的写法为

void multiply()
{
    for (int i = 0; i < n; ++i)
        for (int k = 0; k < n; ++k)
            if (x[i][k])
            for (int j = 0; j < n; ++j)
                t[i][j] += x[i][k] * y[k][j];
}

我们研究的问题主要是对阶数不大的矩阵进行乘法运算，因此朴素的矩阵乘法效果更好。
对于规模较大的矩阵（例如阶数在512阶及以上），用\(Strassen\)矩阵乘法可以将矩阵相乘的时间复杂度从\(O(t^3)\)降低到\(O(t^{2.81})\)，但是过程较为繁琐，一般作为固定的模板使用。
具体过程可以参考下面这些网页的介绍：
矩阵乘法Strassen算法-知乎
快速矩阵乘法的研究——上-csdn

矩阵快速幂(Fast Power)

求矩阵的\(n\)次幂可以通过迭代的方法加速进程，具体操作如下：
求\(A^n\)，其中\(A\)是\(t\)阶方阵。
我们知道，任意的正整数\(n\)都可以写成

\[n=2^{k_1}+2^{k_2}+\cdots+2^{k_t} \]

其中\(k_1>k_2>\cdots>k_t\geq0\)
通过迭代，可以快速求得\(A^n\):

1. \(A^{2^k}=\left(A^{2^{k-1}}\right)^2\)
2. \(A^{2^a+2^b}=A^{2^a}×A^{2^b}\)

这样的办法叫做矩阵快速幂。

时间复杂度为\(O(t^3 \log⁡n)\).

相似对角矩阵法

若方阵\(A\)可相似对角化，当且仅当\(A\)有\(t\)个线性无关特征向量。
存在可逆矩阵\(P\)满足\(D=P^{-1}AP\)
其中\(D=\left[\begin{matrix}λ_1&⋯&0\\⋮&⋱&⋮\\0&⋯&λ_n\end{matrix}\right]\),\(λ_1,⋯,λ_t\)是\(A\)的\(t\)个线性无关特征值。
那么则有

\[A=PDP^{-1}⇒A^n=(PDP^{-1})^n=PD^nP^{-1} \]

其中\(D^n=\left[\begin{matrix}λ_1^n&⋯&0\\⋮&⋱&⋮\\0&⋯&λ_n^n\end{matrix}\right]\)

求解特征值的方法

\(A\)的特征多项式定义为\(p(λ)=|λE_t-A|\).
求解方程\(p(λ)=0\)可以解出\(t\)个特征值\(λ_1,⋯,λ_t\)。
以斐波那契\((Fibonacci)\)数列为例，

\[p(\lambda)=\left|\begin{matrix}λ-1&-1\\-1&λ\end{matrix}\right|=\lambda^2-\lambda-1=0 \]

\[\Rightarrow\lambda_1=\frac{1+\sqrt5}{2},\lambda_2=\frac{1-\sqrt5}{2} \]

求解特征向量的方法

将\(t\)个特征值\(λ_1,⋯,λ_t\)中不相同的特征值分别代入\(B=λE_t-A\)得到\(t'(0<t'\leq t)\)个矩阵。
对于每一个矩阵\(B_i(0<i\leq t')\)，进行初等行变换得到简单方阵\(B'_i\)，此时可以计算出向量空间的维数\(dim(B'_i)=k_i=t-R(B'_i)\)，其中\(R(B'_i)\)表示矩阵\(B'_i\)的秩。
将矩阵\(B'_i\)写作\(t\)元线性方程（组）

\[B'_i\left[\begin{matrix}x_{1}\\⋮\\x_{t}\end{matrix}\right]=0 \]

即

\[\begin{cases}b_{11} x_1+b_{12} x_2+⋯+b_{1t} x_t=0\\⋮\\b_{k_i1} x_1+b_{k_i2} x_2+⋯+b_{k_it} x_t=0\end{cases} \]

可以计算出\(k_i\)个线性无关的特征向量，任意代入其中的\(R(B')\)个数据即可求得\(k_i\)个线性无关的列向量解。
对\(t'\)个矩阵\(B_i(0<i\leq t')\)进行相同操作，求得\(t\)个线性无关的列向量解。

\[\left[\begin{matrix}x_{11}\\⋮\\x_{t1}\end{matrix}\right],\cdots,\left[\begin{matrix}x_{1t}\\⋮\\x_{tt}\end{matrix}\right] \]

将这些得到的\(t\)个列向量依次结合，可以得到可逆矩阵\(P\)，将\(t\)个特征值依次排在主对角线上，可以得到\(D\)。
以斐波那契\((Fibonacci)\)数列为例，
将\(\lambda_1\)代入得到

\[B_1=λ_1E_2-A=\left[\begin{matrix}-λ_2&-1\\-1&λ_1\end{matrix}\right] \]

将\(\lambda_2\)代入得到

\[B_2=λ_2E_2-A=\left[\begin{matrix}-λ_1&-1\\-1&λ_2\end{matrix}\right] \]

将\(B_1\)化为简单矩阵\(B'_1\)

\[B'_1=\left[\begin{matrix}λ_2&1\\0&0\end{matrix}\right] \]

解方程\(B'_1T_1=0\)，得到特征向量

\[T_1=\left[\begin{matrix}1\\-λ_2\end{matrix}\right] \]

同样地，将\(B_2\)化为简单矩阵\(B'_2\)，并求出特征向量\(T_2\)

\[B'_2=\left[\begin{matrix}λ_1&1\\0&0\end{matrix}\right],T_2=\left[\begin{matrix}λ_2\\1\end{matrix}\right] \]

将\(T_1,T_2\)依次结合，即可得到可逆矩阵

\[P=\left[\begin{matrix}1&λ_2\\-λ_2&1\end{matrix}\right] \]

特征值依次排列在主对角线上，得到

\[D=\left[\begin{matrix}λ_1\\&λ_2\end{matrix}\right] \]

求解逆矩阵的方法

要求可逆矩阵\(P\)的逆矩阵\(P^{-1}\)，方法有很多种，如伴随矩阵法，初等行变换法，哈密尔顿-凯莱定理法等。
以斐波那契\((Fibonacci)\)数列为例，我们求得

\[P^{-1}=\frac{1}{\sqrt5}\left[\begin{matrix}λ_1&1\\-1&λ_1\end{matrix}\right] \]

因此我们得到

\[\begin{aligned} A^{n-1}&=PD^{n-1}P^{-1} \\&=\frac{1}{\sqrt5}\left[\begin{matrix}1&λ_2\\-λ_2&1\end{matrix}\right]\left[\begin{matrix}λ_1^{n-1}\\&λ_2^{n-1}\end{matrix}\right]\left[\begin{matrix}λ_1&1\\-1&λ_1\end{matrix}\right] \\&=\frac{1}{\sqrt5}\left[\begin{matrix}λ_1^{n}-λ_2^{n}&λ_1^{n-1}-λ_2^{n-1}\\λ_1^{n-1}-λ_2^{n-1}&λ_1^{n-2}-λ_2^{n-2}\end{matrix}\right] \end{aligned}\]

可以得到

\[\begin{aligned} f[n]&=\frac{1}{\sqrt5}\left[\left(λ_1^{n}-λ_2^{n}\right)f[1]+\left(λ_1^{n-1}-λ_2^{n-1}\right)f[0]\right]\\ &=\frac{1}{\sqrt5}\left(λ_1^{n+1}-λ_2^{n+1}\right)\\ &=\cfrac{\left(\frac{1+\sqrt5}{2}\right)^{n+1}-\left(\frac{1-\sqrt5}{2}\right)^{n+1}}{\sqrt5} \end{aligned}\]

时间复杂度

上述方法矩阵运算需要自行计算方可得到通解，根据通解求出\(f[n]\)时间复杂度为\(O(1)\)。
这种方法十分繁琐，一般更好适用于在有计算软件辅助下进行较大规模的矩阵运算。

哈密尔顿-凯莱定理(Cayley–Hamilton theorem)法

\(A\)的特征多项式定义为\(p(λ)=|λE_t-A|\).
哈密尔顿-凯莱定理断言：\(p(A)=\mathbf{0}\)
展开得\(A^t+c_{t-1} A^{t-1}+⋯+c_1 A+c_0 E_t=\mathbf{0}\)
通过这种方式可以将\(A^n (n≥t)\)降次为\(c'_{t-1} A^{t-1}+⋯+c'_1 A+c'_0 E_t\)的形式。
可以假设\(A^n=c'_{t-1} A^{t-1}+⋯+c'_1 A+c'_0 E_t\)
结合方程

\[\begin{cases}λ_1^n=c'_{t-1} λ_1^{t-1}+⋯+c'_1 λ_1+c'_0\\⋮\\λ_t^n=c'_{t-1} λ_t^{t-1}+⋯+c'_1λ_t+c'_0 \end{cases} \]

可以分别解得

\[\left[\begin{matrix}c'_{t-1}\\⋮\\c'_1\\c'_0\end{matrix}\right] \]

即可得到\(A^n\)并求出通项。
以斐波那契\((Fibonacci)\)数列为例，
欲求出\(A^{n-1}\)，我们列方程组

\[\begin{cases}λ_1^{n-1}=c'_1 λ_1+c'_0\\λ_2^{n-1}=c'_1λ_2+c'_0 \end{cases} \]

解得

\[\begin{cases}c'_1=\frac{1}{\sqrt5}\left(λ_1^{n-1}-λ_2^{n-1}\right)\\c'_0=\frac{1}{\sqrt5}\left(λ_1^{n-2}-λ_2^{n-2}\right)\end{cases} \]

所以

\[\begin{aligned} A^{n-1}&=c'_1\left[\begin{matrix}1&1\\1&0\end{matrix}\right]+c'_0\left[\begin{matrix}1&0\\0&1\end{matrix}\right]\\ &=\left[\begin{matrix}c'_1+c'_0&c'_1\\c'_1&c'_0\end{matrix}\right] \end{aligned}\]

解得通项

\[\begin{aligned} f[n]&=(c'_1+c'_0)f[1]+c'_1f[0]\\ &=\cfrac{\left(\frac{1+\sqrt5}{2}\right)^{n+1}-\left(\frac{1-\sqrt5}{2}\right)^{n+1}}{\sqrt5} \end{aligned}\]

用该方法进行矩阵运算比相似对角矩阵法要简便，同样需要自行计算方可得到通解，根据通解求出\(f[n]\)时间复杂度为\(O(1)\)。

参考资料

矩阵快速幂（练习）
https://www.luogu.com.cn/problem/P3390
矩阵快速幂
https://zhuanlan.zhihu.com/p/95902286
矩阵n次方计算法
矩阵相似对角化
https://zhuanlan.zhihu.com/p/138285148
哈密尔顿-凯莱定理
https://baike.baidu.com/item/哈密尔顿－凯莱定理

感谢支持！