抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

求斐波那契数列通项-线性代数大作业

发布于：2023年1月10日

设斐波那契数列第 kk 项为 FkF_k，则递推公式：

Fk=Fk−1+Fk−2F_k=F_{k-1}+F_{k-2}

则：

[FkFk−1]=[1110][Fk−1Fk−2]\begin{bmatrix} F_k\\ F_{k-1} \end{bmatrix} = \begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix} \begin{bmatrix} F_{k-1}\\ F_{k-2} \end{bmatrix}

则：

[FkFk−1]=[1110]k−1[F1F0]\begin{bmatrix} F_k\\ F_{k-1} \end{bmatrix} = \begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix}^{k-1} \begin{bmatrix} F_1\\ F_0 \end{bmatrix}

那么：

Fk=[1110]1,1k−1F_k=\begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix}^{k-1}_{1,1}

问题转化为求一个矩阵的 nn 次幂。

由于矩阵运算满足结合律，我们可以采用矩阵快速幂。

1
2
3
4
5
6
7
8
9
Matrix operator ^ (Matrix A,int k){
    Matrix base=A;
    while (k){
        if (k&1) base=base*A;
        A=A*A;
        k>>=1;
    }
    return base;
}

编写代码如下

1
2
3
4
5
6
7
8
9
10
11
#include <fraction.h>
#define Num frac
#include <matrix.h>
using namespace std;
int main(){
    Matrix A(2,2);
    A[1][1]=1,A[1][2]=1,A[2][1]=1;
    for (int i=1;i<=10;++i){
        cout<<"Fib"<<i<<"="<<(A^(i-1))[1][1]<<endl;
    }
}

如果要求解析解，我们可以采用相似对角化，求出 PP，使得 P−1AP=ΛP^{-1}AP=\Lambda，则：

An=(PΛP−1)n=PΛnP−1A^n=(P\Lambda P^{-1})^n=P\Lambda^nP^{-1}

首先，计算特征值：

∣A−λE∣=0|A-\lambda E|=0

∣1−λ11−λ∣=−λ+λ2−1\begin{vmatrix} 1-\lambda &1\\ 1 &-\lambda\\ \end{vmatrix} =-\lambda+\lambda^2-1

有两个解：

λ1=12+125,λ2=12−125\lambda_1=\frac{1}{2}+\frac{1}{2}\sqrt{5},\lambda_2=\frac{1}{2}-\frac{1}{2}\sqrt{5}

对应 C++ 代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Matrix A(2,2);
A.Message="A";
A[1][1]=1,A[1][2]=1,A[2][1]=1;
out<<begin_latex()<<endl<<endl;
out<<to_latex(A)<<endl<<endl;
Matrix B=A;
for (int i=1;i<=A.row;++i){
    B[i][i]=B[i][i]-Num("l");
}
cout<<"Eigen Poly"<<endl;
_poly x=Determinant(B).x;//计算行列式
cout<<x<<endl;
out<<to_latex(B,"vmatrix")<<"="<<to_latex(x)<<endl<<endl;
upoly _x;
_x.init_from_poly(x);
cpoly v=Factorization(_x);//因式分解为一次和二次因式
v.sort();
cout<<"Factorization"<<endl<<v<<endl;
out<<"Factorize it then we have "<<to_latex(v)<<endl<<endl;

对于每一个特征值，用基础解系求得一个特征向量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
for (int i=0;i<v.v.size();++i){
    if (v.v[i].first.deg()==1){
        Num lambda=poly(poly_ele((frac)(0)-v.v[i].first[0]));
        out<<"For eigen value $\\lambda="<<to_latex(lambda,0)<<"$ ,we have vectors,"<<endl<<endl;
        cout<<"lambda="<<lambda<<endl;
        cout<<"n="<<v.v[i].second<<endl;//代数重数
        Matrix B=A-lambda*Matrix(A.row,A.col,1);
        vector<Matrix>baseS=baseSolution(B);
        // baseS=Schmidt(baseS); 如果需要正交矩阵的话需要施密特正交化
        cout<<"m="<<baseS.size()<<endl;//几何重数，几何重数不超过代数重数
        cout<<B<<endl;
        for (int i=0;i<baseS.size();++i){
            cout<<baseS[i]<<endl;
            s.push_back(baseS[i]);
            out<<to_latex(baseS[i])<<endl<<endl;
        }
    }
    else if (v.v[i].first.deg()==2){
        Num lambda_1,lambda_2;
        frac a=v.v[i].first[2],b=v.v[i].first[1],c=v.v[i].first[1];
        frac delta=b*b-4*a*c;
        lambda_1=poly_ele(sqrtNum(-b/(2*a),1/(2*a),delta));
        lambda_2=poly_ele(sqrtNum(-b/(2*a),-1/(2*a),delta));
        out<<"For eigen value $\\lambda="<<to_latex(lambda_1,0)<<"$ ,we have vectors,"<<endl<<endl;
        cout<<"lambda="<<lambda_1<<endl;
        cout<<"n="<<v.v[i].second<<endl;//代数重数
        Matrix B=A-lambda_1*Matrix(A.row,A.col,1);
        vector<Matrix>baseS=baseSolution(B);
        cout<<"m="<<baseS.size()<<endl;//几何重数，几何重数不超过代数重数
        cout<<B<<endl;
        for (int i=0;i<baseS.size();++i){
            cout<<baseS[i]<<endl;
            s.push_back(baseS[i]);
            out<<to_latex(baseS[i])<<endl<<endl;
        }
        out<<"For eigen value $\\lambda="<<to_latex(lambda_2,0)<<"$ ,we have vectors,"<<endl<<endl;
        cout<<"lambda="<<lambda_2<<endl;
        cout<<"n="<<v.v[i].second<<endl;//代数重数
        B=A-lambda_2*Matrix(A.row,A.col,1);
        baseS=baseSolution(B);
        cout<<"m="<<baseS.size()<<endl;//几何重数，几何重数不超过代数重数
        cout<<B<<endl;
        for (int i=0;i<baseS.size();++i){
            cout<<baseS[i]<<endl;
            s.push_back(baseS[i]);
            out<<to_latex(baseS[i])<<endl<<endl;
        }
    }
}

得到：

η1=[12+1251],η2=[12−1251]\eta_1=\begin{bmatrix} \frac{1}{2}+\frac{1}{2}\sqrt{5}\\ 1\\ \end{bmatrix},\eta_2=\begin{bmatrix} \frac{1}{2}-\frac{1}{2}\sqrt{5}\\ 1\\ \end{bmatrix}

组合形成

P=[12+12512−12511]P=\begin{bmatrix} \frac{1}{2}+\frac{1}{2}\sqrt{5} &\frac{1}{2}-\frac{1}{2}\sqrt{5}\\ 1 &1\\ \end{bmatrix}

且：

Λ=PAP−1=[12+1250012−125]\Lambda=PAP^{-1}=\begin{bmatrix} \frac{1}{2}+\frac{1}{2}\sqrt{5} &0\\ 0 &\frac{1}{2}-\frac{1}{2}\sqrt{5}\\ \end{bmatrix}

解得：

Fn=(12+1105)(12+125)n−1+(12−1105)(12−125)n−1F_n=(\frac{1}{2}+\frac{1}{10}\sqrt{5})(\frac{1}{2}+\frac{1}{2}\sqrt{5})^{n-1}+(\frac{1}{2}-\frac{1}{10}\sqrt{5})(\frac{1}{2}-\frac{1}{2}\sqrt{5})^{n-1}

对应代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Matrix P=s[0];
out<<"Combine the vectors together then we get, "<<endl<<endl;
for (int i=1;i<s.size();++i){
    P=addH(P,s[i]);
}
cout<<P.message("P")<<endl;
cout<<(P^(-1)).message("P^{-1}")<<endl;
cout<<((P^-1)*A*P).message("P^{-1}AP")<<endl;
Matrix Lambda=(P^-1)*A*P;
out<<to_latex(P.message("P"))<<endl<<endl;
out<<to_latex(((P^-1)*A*P).message("\\Lambda=PAP^{-1}"))<<endl<<endl;
out<<"So $A^n=P^{-1}\\Lambda^n P$"<<endl<<endl;
out<<"For calculating $A^n$, assume $\\Lambda^n$=[a,0;0,b], then we have:"<<endl<<endl;
Matrix _A(2,2);
_A[1][1]=poly("a"),_A[2][2]=poly("b");
cout<<P*_A*(P^-1)<<endl;
out<<to_latex((P*_A*(P^-1))[1][1])<<endl<<endl;
out<<"Then, substitute $a=("<<to_latex(Lambda[1][1],0)<<")^{n-1}$ and $b=("<<to_latex(Lambda[2][2],0)<<")^{n-1}$, we have the final answer"<<endl<<endl;

一些小小的推广

求 Fk=Fk−1+2Fk−2F_k=F_{k-1}+2F_{k-2} 的通项。

则：

A=[1210]A=\begin{bmatrix}1&2\\1 &0\end{bmatrix}

利用程序自动生成结果：

求 Fk=−2Fk−1+Fk−2+2Fk−3F_k=-2F_{k-1}+F_{k-2}+2F_{k-3} 的通项。

则：

A=[−212100010]A=\begin{bmatrix} -2 &1 &2\\ 1 &0 &0\\ 0 &1 &0\\ \end{bmatrix}

image-20221215155720027

注意到，程序算出的 AA 的特征多项式就是我们常说的数列的特征多项式，在后面严格证明。

更进一步推广：常系数齐次线性递推

一个 kk 阶常系数齐次递推数列满足：

fn=∑i=1kaifn−if_n=\sum_{i=1}^ka_if_{n-i}

对于斐波那契数列：

a1=a2=1a_1=a_2=1

转移矩阵：

A=(a1a2a3…ak−2ak−1ak100…000010⋯000⋮⋮⋮⋱⋮⋮⋮000⋯100000⋯010)A=\begin{pmatrix}a_1&a_2&a_3&\dots&a_{k-2}&a_{k-1}&a_k\\1&0&0&\dots&0&0&0\\0&1&0&\cdots&0&0&0\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\0&0&0&\cdots&1&0&0\\0&0&0&\cdots&0&1&0\end{pmatrix}

初值：

F=(fk−1fk−2⋮f0)F=\begin{pmatrix}f_{k-1}\\f_{k-2}\\\vdots\\f_0\end{pmatrix}

则答案变成 (AnF)k,1(A^nF)_{k,1}。还是转化为计算 AnA^n。

不妨设多项式

g(λ)=∣λE−A∣g(\lambda)=|\lambda E-A|

在第一行第一列展开，我们得到 λk\lambda^k，在第一行第二列展开，我们得到 −λk−1a1-\lambda^{k-1}a_1，在第一行第三列展开，我们得到 −λk−2a2-\lambda^{k-2}a_2，之后一直到 −λ0ak-\lambda^0a_k。

于是

g(λ)=λk−a1λk−1−⋯−akg(\lambda)=\lambda^k-a_1\lambda^{k-1}-\cdots-a_k

由 Cayley-Hamilton 定理，我们也知道 g(A)=0g(A)=0，那么我们希望凑出 g(A)g(A) 的形式，就可以进行抵消。

不妨设 An=f(A)g(A)+r(A)A^n=f(A)g(A)+r(A)，其中 f,rf,r 都是关于 AA 的多项式。那么就可以转化为多项式除法和取余的问题。答案是 r(A)r(A)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <poly.h>
#define Num frac
#include <matrix.h>
#define MAXN 1005
using namespace std;
frac a[MAXN];
int main(){
    int n,k;
    cin>>n>>k;
    upoly g;
    g.v.resize(k+1);
    g[k]=frac(1);
    for (int i=1;i<=k;++i){
        cin>>a[i];
        g[k-i]=-a[i];
    }
    Matrix A(k,k);
    for (int i=1;i<=k;++i) A[1][i]=a[i];
    for (int i=1;i<=k-1;++i) A[i+1][i]=1;
    upoly M;
    M.v.resize(n+1);
    M[n]=frac(1);
    upoly r=M%g;
    Matrix An=Matrix(k,k);
    for (int i=r.v.size()-1;i>=0;--i) An=An*A+r.v[i]*Matrix(k,k,1);
    Matrix F(k,1);
    for (int i=1;i<=k;++i) cin>>F[k-i+1][1];
    cout<<(An*F)[k][1]<<endl;
}

例如，计算 fib8fib_{8}，则输入

1
2
3
8 2
1 1
0 1

不用线性代数的另解：生成函数

构造f(x)=∑j=1fibjxjf(x)=\sum _{j=1} fib_j x^j

所以f(x)×x=∑j=2fibj−1xjf(x) \times x=\sum_{j=2} fib_{j-1}x^j

f(x)−f(x)×x=fib1x+∑j=3fibj−2xj=x+x2×f(x)f(x)-f(x)\times x=fib_1x+\sum_{j=3}^{}fib_{j-2}x^j=x+x^2\times f(x)

所以f(x)=x1−x−x2f(x)=\frac{x}{1-x-x^2}

令ϕ=1+52,ϕ′=1−52\phi =\frac{1+\sqrt{5}}{2},\phi'=\frac{1-\sqrt{5}}{2}。

分解：x1−x−x2=x(1−ϕ′x)(1−ϕx)=15(11−ϕx−11−ϕ′x)\frac{x}{1-x-x^2}=\frac{x}{(1-\phi'x)(1-\phi x)}=\frac{1}{\sqrt{5}}(\frac{1}{1-\phi x} - \frac{1}{1-\phi' x})

根据公式11−x=∑j=0xj\frac{1}{1-x}=\sum _{j=0} x^{j}，前面一项化为(ϕ)n(\phi)^n后面一项化成(ϕ′)n(\phi')^n。

于是我们得到fibn=15(ϕn−ϕ′n)fib_n=\frac{1}{\sqrt{5}}(\phi ^n-\phi'^n)。

即fibn=−15(1−52)n+15(1+52)nfib_n=-\frac{1}{\sqrt 5}(\frac{1-\sqrt 5}{2})^n+\frac{1}{\sqrt 5}(\frac{1+\sqrt 5}{2})^n。

相关文章
P4783 【模板】矩阵求逆P7112 【模板】行列式求值QR分解n维向量与线性方程组解的结构基本概念

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：https://gaisaiyuno.github.io/archives/2e3b6d72.html

更新于：2023年1月10日

线性代数

评论