数值分析——带双步位移的QR分解求特征值算法

数 值 分 析（B ） 大 作 业（二）

1、算法设计：

①矩阵的拟上三角化：

对实矩阵A 进行相似变换化为拟上三角矩阵(1)A n -，其变换矩阵采用Householder 矩阵，变换过程如下：

若()

a 0(2,,)r ir i r n ==+ ，则r H I =；

否则，(r)(r)T

r r+1,r n,r s =(0,,0,a ,a ) ，

(r)

r r+1,r r 2c = -sgn(a )||s ||，

()()()

r r r r+11,2,u =s -c e (0,,0,,,,)r r r r r r r r nr a c a a ++=- ，

2T

r r r 2

H =I-2u u /r u ，

(1)

()

r r r r A

H A

H +=。

当2r n =-时，得(1)(2)(1)

222112n n n n n n A

H A H H H A H H ------=== ，令12n-2P=H H H 又r H 是对称正交矩阵，于是n-1T A =P AP 成立，因而n-1A 与 A 相似。

②矩阵的QR 分解：

矩阵的QR 分解过程与拟上三角化过程相似，在这里不再重复其原理。 ③求全部特征值

矩阵拟上三角化后利用带双步位移的QR 方法，采用书本Page 63页具体算法实现。为了使编程方便，并体会goto 语句使用的灵活性，程序中的跳转均使用goto Loop*实现。

④求A 的相应于实特征值的特征向量

求实特征值对应的特征向量，即是求解线性方程组(λI-A)=0i i x ，

(1,,)i n = 。因此，为得到全部实特征值对应的特征向量，解线性方程组的过

程要循环n 次（n 为矩阵阶数）。线性方程组的求解采用列主元素Gauss 消去法。

#include

#include

#define ERR 1.0e-12 //误差限

#define N10 //矩阵行列数

#define L 1.0e5 //最大迭代次数

double A[N][N]={0};

void Init_A() //初始化矩阵

{

int i,j;

for(i=0;i

for(j=0;j

{

if(i==j)

A[i][j]=1.5*cos((i+1)+1.2*(j+1));

else

A[i][j]=sin(0.5*(i+1)+0.2*(j+1));

}

}

/*

void Display_A()

{

int i,j;

for(i=0;i

{

for(j=0;j

printf("%8.4f",A[i][j]);

printf("\n");

}

}

*/

int Sgn(double a)

{

if(a>=0)

return 1;

else

return -1;

}

void On_To_The_Triangle() //矩阵拟上三角化{

int i,j,r,flag=0;

double cr,dr,hr,tr,temp;

double ur[N],pr[N],qr[N],wr[N];

for(r=1;r<=N-2;r++)

{

flag=0;

for(i=r+2;i<=N;i++)

if(A[i-1][r-1]!=0)

{

flag=1;

break;

}

if(0==flag)

continue;

dr=0;

for(i=r+1;i<=N;i++)

dr+=A[i-1][r-1]*A[i-1][r-1];

dr=sqrt(dr);

if(0==A[r][r-1])

cr=dr;

else cr=-Sgn(A[r][r-1])*dr;

hr=cr*cr-cr*A[r][r-1];

for(i=1;i<=r;i++)

ur[i-1]=0;

ur[r]=A[r][r-1]-cr;

for(i=r+2;i<=N;i++)

ur[i-1]=A[i-1][r-1];

for(i=1;i<=N;i++)

{

temp=0;

for(j=1;j<=N;j++)

temp+=A[j-1][i-1]*ur[j-1];

pr[i-1]=temp/hr;

}

for(i=1;i<=N;i++)

{

temp=0;

for(j=1;j<=N;j++)

temp+=A[i-1][j-1]*ur[j-1];

qr[i-1]=temp/hr;

}

temp=0;

for(i=1;i<=N;i++)

{

temp+=pr[i-1]*ur[i-1];

tr=temp/hr;

}

for(i=1;i<=N;i++)

{

wr[i-1]=qr[i-1]-tr*ur[i-1];

}

for(i=1;i<=N;i++)

for(j=1;j<=N;j++)

A[i-1][j-1]=A[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];

}

}

void Get_Roots(double eigenvalue[][2],int m,double ss,double tt) //求一元二次方程的根{

double discriminant=ss*ss-4*tt; //

if(discriminant<0)

{

*(*(eigenvalue+m-2))=0.5*ss;

*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);

*(*(eigenvalue+m-1))=0.5*ss;

*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);

}

else

{

*(*(eigenvalue+m-2))=0.5*(ss+sqrt(discriminant));

*(*(eigenvalue+m-2)+1)=0;

*(*(eigenvalue+m-1))=0.5*(ss-sqrt(discriminant));

*(*(eigenvalue+m-1)+1)=0;

}

}

void Get_Mk(double mk[][N],int m,double ss,double tt) //获取Mk,用于带双步位移的QR分解

{

int i,j,k;

for(i=0;i

for(j=0;j

*(*(mk+i)+j)=0;

for(i=0;i

for(j=0;j

{

for(k=0;k

*(*(mk+i)+j)+=A[i][k]*A[k][j];

*(*(mk+i)+j)-=ss*A[i][j];

if(j==i)

*(*(mk+i)+j)+=tt;

}

}

void QR_Reslove(double mk[][N],int m) //QR分解{

int i,j,r,flag=0;

double cr,dr,hr,tr,temp;

double ur[N],vr[N],pr[N],qr[N],wr[N];

double B[N][N],C[N][N];

for(i=0;i

for(j=0;j

{

B[i][j]=*(*(mk+i)+j);

C[i][j]=A[i][j];

}

for(r=1;r<=m-1;r++)

{

flag=0;

for(i=r+1;i<=m;i++)

if(B[i-1][r-1]!=0)

{

flag=1;

break;

}

if(0==flag)

continue;

dr=0;

for(i=r;i<=m;i++)

dr+=B[i-1][r-1]*B[i-1][r-1];

dr=sqrt(dr);

if(0==B[r-1][r-1])

cr=dr;

else cr=-Sgn(B[r-1][r-1])*dr;

hr=cr*cr-cr*B[r-1][r-1];

for(i=1;i

ur[i-1]=0;

ur[r-1]=B[r-1][r-1]-cr;

for(i=r+1;i<=m;i++)

ur[i-1]=B[i-1][r-1];

for(i=1;i<=m;i++)

{

temp=0;

for(j=1;j<=m;j++)

temp+=B[j-1][i-1]*ur[j-1];

vr[i-1]=temp/hr;

}

for(i=0;i

for(j=0;j

B[i][j]-=ur[i]*vr[j];

for(i=1;i<=m;i++)

{

temp=0;

for(j=1;j<=m;j++)

temp+=C[j-1][i-1]*ur[j-1];

pr[i-1]=temp/hr;

}

for(i=1;i<=m;i++)

{

temp=0;

for(j=1;j<=m;j++)

temp+=C[i-1][j-1]*ur[j-1];

qr[i-1]=temp/hr;

}

temp=0;

for(i=1;i<=m;i++)

{

temp+=pr[i-1]*ur[i-1];

tr=temp/hr;

}

for(i=1;i<=m;i++)

{

wr[i-1]=qr[i-1]-tr*ur[i-1];

}

for(i=1;i<=m;i++)

for(j=1;j<=m;j++)

C[i-1][j-1]=C[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];

}

for(i=0;i

for(j=0;j

A[i][j]=C[i][j];

}

void Display_Eigenvalue(double value[][2]) //显示特征值

{

int i;

for(i=0;i

{

printf("λ%d=%8.4f",i+1,*(*(value+i)));

if(*(*(value+i)+1)>0)

printf("+%8.4f",*(*(value+i)+1));

else if(*(*(value+i)+1)<0)

printf("%8.4f",*(*(value+i)+1));

printf("\n");

}

printf("\n");

}

int QR_With_Double_Step_Displacement(double eigenvalue[][2]) //带双步位移QR分解求特征值

{

int i,j,k=1,m=N;

double s,t;

double Mk[N][N];

for(i=0;i

for(j=0;j<2;j++)

eigenvalue[i][j]=0;

do

{

k++;

if(m==1)

{

eigenvalue[m-1][0]=A[m-1][m-1];

m--;

continue;

}

else if(m==2)

{

s=A[m-2][m-2]+A[m-1][m-1];

t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-1][m-2];

Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根

m=0;

continue;

}

else if(m==0)

return 0;

else if(fabs(A[m-1][m-2])<=ERR)

{

eigenvalue[m-1][0]=A[m-1][m-1];

m--;

continue;

}

else

{

s=A[m-2][m-2]+A[m-1][m-1];

t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];

if(fabs(A[m-2][m-3])<=ERR)

{

Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根

m-=2;

continue;

}

else

{

Get_Mk(Mk,m,s,t); //获取Mk,用于带双步位移的QR分解

QR_Reslove(Mk,m); //QR分解

}

}

}

while(k

printf("Errer\n"); //迭代过程不成功，报错

}

void Select_Principal_Element(int k) //Gauss消元法求解方程组之选主元

{

int i,max=k;

double trans[N];

for(i=k+1;i

if(fabs(A[i][k])>fabs(A[max][k]))

max=i;

if(k==max)

return;

else

{

for(i=k;i

{

trans[i]=A[k][i];

A[k][i]=A[max][i];

A[max][i]=trans[i];

}

}

}

void Eliminant(int k) //Gauss消元法求解方程组之消元

{

int i,j;

double rate;

for(i=k+1;i

{

rate=A[i][k]/A[k][k];

for(j=k;j

A[i][j]=A[i][j]-rate*A[k][j];

}

}

void Back_Substitution(double Eigenvector[]) //Gauss消元法求解方程组之回代

{

int i,j,k;

double Temp_M[N][N+1];

for(i=0;i

{

for(j=i;j

Temp_M[i][j]=A[i][j];

Temp_M[i][N]=0;

}

Temp_M[N-1][N-1]=1;

Temp_M[N-1][N]=1;

for(k=N-1;k>=0;k--)

{

*(Eigenvector+k)=Temp_M[k][N];

for(i=k;i

*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_M[k][i+1];

*(Eigenvector+k)=*(Eigenvector+k)/Temp_M[k][k];

}

}

void Display_Eigenvector(double Eigenvector[],double eigenvalue) //显示特征值对应特征向量

{

int i;

printf("When λ=%8.4f\nEigenvector=\n",eigenvalue);

for(i=0;i

printf("%8.4f",*(Eigenvector+i));

printf("\n");

}

void Get_Eigenvector(double value[][2]) //利用Gauss消元法求解特征向量{

int i,j;

double eigenvalue;

double Eigenvector[N];

for(i=0;i

if(*(*(value+i)+1)==0)

{

Init_A(); //初始化矩阵

eigenvalue=*(*(value+i));

for(j=0;j

A[j][j]-=eigenvalue;

for(j=0;j

{

Select_Principal_Element(j); //Gauss消元法求解方程组之选主元

Eliminant(j); //Gauss消元法求解方程组之消元}

Back_Substitution(Eigenvector); //Gauss消元法求解方程组之回代

Display_Eigenvector(Eigenvector,eigenvalue); //显示特征值对应特征向量}

}

main()

{

double eigenvalue[N][2];

Init_A(); //初始化矩阵

On_To_The_Triangle(); //矩阵上三角化

QR_With_Double_Step_Displacement(eigenvalue); //带双步位移QR分解求特征值printf("Contact me : 731363227@https://www.360docs.net/doc/de14185695.html, \n\n");

Display_Eigenvalue(eigenvalue); //显示特征值

Get_Eigenvector(eigenvalue); //利用Gauss消元法求解特征向量

return 0;

}

PS: For further details, please contact me :731363227@https://www.360docs.net/doc/de14185695.html,