二阶常系数线性差分方程
齐次差分方程
二阶常系数线性差分方程的一般形式为:
yn+2+ayn+1+byn=f(n),n=0,1,2,....(1−1)
其中a,b为已知常数,且b≠0,f(n)为已知函数。
方程(1−1)的对应齐次方程为:
yn+2+ayn+1+byn=0................................(1−2)
设yn=λn为方程(1−2)的特解,其中λ为非零待定系数,代入方程后,有
λn(λ2+aλ+b)=0
因λ≠0,故函数yn=λn是方程(1−2)的特解的充分必要条件是λ满足方程
λ2+aλ+b=0............................................(1−3)
方程(1-3)成为方程(1-1)或(1-2)的特征方程,特征方程的解称为特征根或特征值。
根据二次代数方程(1-3)解的三种情况,可以仿照二阶常系数齐次线性微分方程,分别给出方程(1-1)的通解。
1.1特征方程有两个相异实根
即当判别式Δ=a2−4b>0时,方程(1−2)有两个相异实根
λ1=12(−a−Δ−−√),λ2=12(−a+Δ−−√)....(1−4)
于是方程(1-2)有两个特解
yn1=λn1,yn2=λn2
且由yn1yn2=(λ1λ2)n≠常数,知yn1与yn2线性无关,从而得到方程(1−2)的通解
yn=C1λn1+C2λn2,其中λ1,λ2由式(1−4)给出,C1,C2为任意常数。
1.2特征方程有二重根
即当判别式Δ=a2−4b=0时,方程(1−3)有二重根λ1=λ2=−12a,于是方程(1−2)有一个特解yn1=(−12a)n,可验证方程(1−2)有另一特解yn2=n(−12a)n
(注:系数n可参照二阶常系数线性齐次微分方程的通项证明求得)
且由yn1yn2=1n≠常数,知yn1与yn2线性无关,从而得到方程(10−2)的通解
yn=(C1+C2n)(−12a)n,其中C1,C2为任意常数。
1.3特征方程有两个共轭复根
即当判别式Δ=a2−4b<0时,方程(1−3)有两个共轭复根
λ1=12(−a+i−Δ−−−√),λ2=12(−a−i−Δ−−−√)
通过直接验证可知,方程(1-2)有两个特解
yn1=rncosβn,yn2=rnsinβn
其中r=(−a2)2+(−Δ√2)2−−−−−−−−−−−−−−√=b√........(详见其他参考)
斐波那契数列的递推式及通项公式
递推式:yn+2=⎧⎩⎨0,1,yn+1+yn,n=0n=1……..(1−5)n>1
由式(1−2)可知,yn+2−yn+1−yn=0............(1−6)
其中a=-1,b=-1
令yn=λn,式(1−6)化成λn(λ2−λ−1)=0,解得
yn=C1(1+5√2)n+C2(1−5√2)n,由y0=0,y1=1可知
C1=15√,C2=−15√
yn=15√(1+5√2)n−15√(1−5√2)n即为通解
可知,时间复杂度为指数形式。
改进方法:
1.时间复杂度为O(n)的方法:
int Fibonacci(int n) {
if (n<=0) {
return 0;
}
if (n==1) {
return 1;
}
int min=0;
int max=1;
int i=2;
int result=0;
while (i<=n) {
result=min+max;
min=max;
max=result;
++i;
}
return result;
2.时间复杂度为O(logn)的方法:
根据式(1-5),可以化成矩阵的形式:
(yn+2yn+1)=(1110)(yn+1yn)=…=(1110)n+1(y1y0)
因而计算yn就简化为计算矩阵的(n+1)次方,再对矩阵分解,计算矩阵(n−2)/2次方的平方,逐步分解,因而时间复杂度为O(logn)
#include <iostream>
using namespace std;
class Matrix
{
public:
int n;
int **m;
Matrix(int num)
{
m=new int*[num];
for (int i=0; i<num; i++) {
m[i]=new int[num];
}
n=num;
clear();
}
void clear()
{
for (int i=0; i<n; ++i) {
for (int j=0; j<n; ++j) {
m[i][j]=0;
}
}
}
void unit()
{
clear();
for (int i=0; i<n; ++i) {
m[i][i]=1;
}
}
Matrix operator=(const Matrix mtx)
{
Matrix(mtx.n);
for (int i=0; i<mtx.n; ++i) {
for (int j=0; j<mtx.n; ++j) {
m[i][j]=mtx.m[i][j];
}
}
return *this;
}
Matrix operator*(const Matrix &mtx)
{
Matrix result(mtx.n);
result.clear();
for (int i=0; i<mtx.n; ++i) {
for (int j=0; j<mtx.n; ++j) {
for (int k=0; k<mtx.n; ++k) {
result.m[i][j]+=m[i][k]*mtx.m[k][j];
}
}
}
return result;
}
};
int main(int argc, const char * argv[]) {
unsigned int num=2;
Matrix first(num);
first.m[0][0]=1;
first.m[0][1]=1;
first.m[1][0]=1;
first.m[1][1]=0;
int t;
cin>>t;
Matrix result(num);
result.unit();
int n=t-2;
while (n) {
if (n%2) {
result=result*first;
}
first=first*first;
n=n/2;
}
cout<<(result.m[0][0]+result.m[0][1])<<endl;
return 0;
}