数值分析——带双步位移的QR分解求特征值算法.docx
文本预览下载声明
数 值 分 析（B） 大 作 业（二）
1、算法设计：
①矩阵的拟上三角化：
对实矩阵A进行相似变换化为拟上三角矩阵，其变换矩阵采用Householder矩阵，变换过程如下：
若，则；
否则，，
，
，
，
。
当时，得，令又是对称正交矩阵，于是成立，因而与 相似。
②矩阵的QR分解：
矩阵的QR分解过程与拟上三角化过程相似，在这里不再重复其原理。
③求全部特征值
矩阵拟上三角化后利用带双步位移的QR方法，采用书本Page 63页具体算法实现。为了使编程方便，并体会goto语句使用的灵活性，程序中的跳转均使用goto Loop*实现。
④求A的相应于实特征值的特征向量
求实特征值对应的特征向量，即是求解线性方程组，。因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n次（n为矩阵阶数）。线性方程组的求解采用列主元素Gauss消去法。
#include stdio.h
#include math.h
#define ERR 1.0e-12 //误差限
#define N 10 //矩阵行列数
#define L 1.0e5 //最大迭代次数
double A[N][N]={0};
void Init_A() //初始化矩阵
{
int i,j;
for(i=0;iN;i++)
for(j=0;jN;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;iN;i++)
{
for(j=0;jN;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(discriminant0)
{
*(*(eigenvalue+m-2)
显示全部