北航数值分析第二次大作业--QR分解
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《数值分析A》计算实习题目二
姓名
学号
数值分析第二次大作业
一、算法设计方案
首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），算法如课本P61。
使用QR分解法对矩阵A（n-1）进行QR分解，算法如课本P59,
进而求出所得矩阵的Q、R、RQ矩阵。
然后对A（n-1）进行带双步位移的QR分解求矩阵的全部特征值，采用以下几步进行：
第一步：判断是否a m,m-1(k)<=ε ,若不是,则进入第四步。若是，则a m,m-1(k)为特征值，m=m-1，若m=1，则进入第二步，若m=2进入第三步，否则转第四步。
第二步：m=1，则a11(k)为特征值，转向结束步。
第三步：m=2,则可以求出A的两个特征值s1和s2，转向结束步。
第四步：判断是否a m-1,m-2(k)<=ε，若不是，进入第五步。若是，则得到A的两个特征值s1和s2，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。
第五步：判断是否达到循环上限，若达到，则结束，否则进入第六步。
第六步：对A进行双步位移QR分解，这里的算法如课本P64，分解后转向计数步。
计数步：对循环次数进行计数，并转向第一步。
结束步：显示所求得的特征值。
最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17.
二、源程序代码
#include
#include
#include
int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s;
double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10];
double X[10][10],Y[10][10],Qt[10][10],M[10][10];
double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0};
double epsilon=1e-12;
void main()
{
void Quasiuppertriangular(double A[][10]);
void QRdecomposition(double A[][10]);
void DoublestepsQR(double A[][10]);
int i,j;
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
A[i][j]=sin(0.5*(i+1)+0.2*(j+1));
Q[i][j]=0;
AA[i][j]=A[i][j];
}
A[i][i]=1.5*cos(2.2*(i+1));
AA[i][i]=A[i][i];
}
Quasiuppertriangular(A); //调用拟上三角化函数printf( "\n A经过拟上三角化矩阵为：\n\n");
for(i=0;i<10;i++) //输出拟上三角化矩阵
{
for(j=0;j<10;j++)
{
printf("%.12e ",A[i][j]); //输出拟上三角化矩阵
}
printf( "\n\n");
}
QRdecomposition(A); //调用QR分解函数printf( " 进行QR分解后，R矩阵为：\n\n"); //输出R矩阵
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
printf("%.12e ",R[i][j]);
}
printf( "\n\n");
}
printf( " Q矩阵为:\n\n"); //输出Q矩阵
for(i=0;i<10;i++)
{ for(j=0;j<10;j++)
{ printf("%.12e ",Q[i][j]);
}
printf( "\n\n");
}
printf( " RQ矩阵为：\n\n"); //输出RQ矩阵
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{ printf("%.12e ",RQ[i][j]);
}
printf( "\n\n");
}
DoublestepsQR(A); //调用双步位移函数printf( "\n\n 特征值实部依次为：\n\n"); //输出特征值实部
for(j=0;j<10;j++)
{
printf("%.12e ",Re[j]);
printf("\n\n 特征值虚部依次为：\n\n "); //输出特征值虚部
for(j=0;j<10;j++)
{
printf("%.12e ",Im[j]);
}
//按行输出特征向量
printf( "\n\n 按行输出实特征根相应特征向量为：\n\n");
for(i=0;i<10;i++)
{
if(i==1||i==2||i==5||i==6)
{
continue;
}
for(j=0;j<10;j++)
{
printf("%.12e ",X[i][j]);
}
printf( "\n\n");
}
getchar();
}
void Quasiuppertriangular(double A[][10])
{
for(j=0;j<8;j++)
{
for(i=0;i<10;i++)
{
U[i]=0;
P[i]=0;
T[i]=0;
W[i]=0;
}
m=0;
for(i=j+2;i<10;i++)