北航数值分析第二次大作业概述.docx
《北航数值分析第二次大作业概述.docx》由会员分享，可在线阅读，更多相关《北航数值分析第二次大作业概述.docx（27页珍藏版）》请在冰豆网上搜索。
北航数值分析第二次大作业概述
数值分析第二次大作业
姓名：
李潇学号：
SY1303314
题目:
使用带双步位移的QR分解法求矩阵
的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知：
（i,j=1,2,……,10）
1、算法设计：
①矩阵的拟上三角化：
对实矩阵A进行相似变换化为拟上三角矩阵
，其变换矩阵采用Householder矩阵，变换过程如下：
若
，则
；
否则，
，
，
，
，
。
当
时，得
，令
又
是对称正交矩阵，于是
成立，因而
与
相似。
②矩阵的QR分解：
矩阵的QR分解过程与拟上三角化过程相似，在这里不再重复其原理。
③求全部特征值
矩阵拟上三角化后利用带双步位移的QR方法，采用书本Page63页具体算法实现。
为了使编程方便，并体会goto语句使用的灵活性，程序中的跳转均使用gotoLoop*实现。
④求A的相应于实特征值的特征向量
求实特征值对应的特征向量，即是求解线性方程组
，
。
因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n次（n为矩阵阶数）。
线性方程组的求解采用列主元素Gauss消去法。
2、程序源代码：
#include"Stdio.h"
#include"Conio.h"
#include"math.h"
//****************************************************************************//
//N：
待求矩阵的阶数//
//Epsilon：
精度水平//
//MaxL：
限制求解的最大迭代次数//
//****************************************************************************//
#defineN10
#defineEpsilon0.000000000001
#defineMaxL200
voidCreat_MatrixA（doublearray[N+1][N+1]）;
voidLoad_MatrixA（doubleMatrixA[N+1][N+1],doublearrayA[N+1][N+1]）;
voidClear_Elements（doublearray[N+1][N+1]）;
voidInitial_to_I（doublearray[N+1][N+1]）;
voidNishang（doublearrayA[N+1][N+1]）;
voidQR_Fenjie（doublearrayA[N+1][N+1],doubleQ[N+1][N+1],intm）;
voidSolution_namda（doublearrayA[N+1][N+1],doublenamda[N+1][2]）;
voidSolution_A_Next（doublearrayA[N+1][N+1],intm,doubles,doublet）;
voidSolution_vector（doublearray[N+1][N+1],doublenamda[N+1][2],doublex[N+1][N+1]）;
voidGaussElemitation_select（intm,intflag,doublearray[N+1][N+1],doublex[N+1][N+1]）;
voidPrint_test（doublearray[N+1][N+1]）;
voidPrint_testAB（doublearrayA[N+1][N+1],doublearrayB[N+1][N+1]）;
voidPrint_namda（doublenamda[N+1][2]）;
voidPrint_vector（doublex[N+1][N+1],doublenamda[N+1][2]）;
intmain（void）
{
//****************************************************************************//
//定义主程序中用到的变量//
//****************************************************************************//
//MatrixA：
用来存储源矩阵A//
//Q：
用来存储对A（n-1）进行QR分解得到的正交阵//
//namda：
用来存储矩阵A的全部特征值。
第1列存储实部，第2列存储虚部//
//****************************************************************************//
doubleMatrixA[N+1][N+1],Q[N+1][N+1],namda[N+1][2],x[N+1][N+1];
//****************************************************************************//
//对矩阵A进行拟上三角化，再进行QR分解见//
//****************************************************************************//
Creat_MatrixA（MatrixA）;
printf（"A经过拟上三角化得到的矩阵A（n-1）为：
\n"）;
/*对矩阵A进行拟上三角化，并输出矩阵A（n-1）*/
Nishang（MatrixA）;
Print_test（MatrixA）;
printf（"\nA（n-1）经过QR分解得到的矩阵Q、R和RQ分别为：
\n"）;
/*对矩阵A（n-1）进行QR分解*/
QR_Fenjie（MatrixA,Q,N）;
/*输出分解得到的Q、R和RQ*/
printf（"Q:
\n"）;Print_test（Q）;
printf（"R:
\n"）;Print_test（MatrixA）;
printf（"RQ:
\n"）;Print_testAB（MatrixA,Q）;
//****************************************************************************//
//用带双步位移的QR方法求矩阵A的全部特征值//
//****************************************************************************//
/*重新生成矩阵A*/
Creat_MatrixA（MatrixA）;
/*调用Solution_namda求矩阵A的全部特征值，并输出*/
Solution_namda（MatrixA,namda）;
Print_namda（namda）;
//****************************************************************************//
//求实特征值对应的特征向量//
//****************************************************************************//
/*重新生成矩阵A*/
Creat_MatrixA（MatrixA）;
printf（"\n实特征值对应的特征向量：
"）;
/*调用Solution_vector求矩阵A实特征值对应的特征向量，并输出*/
Solution_vector（MatrixA,namda,x）;
Print_vector（x,namda）;
getch（）;
return0;
}
/*创建源矩阵A*/
voidCreat_MatrixA（doublearray[N+1][N+1]）
{
inti,j;
/*矩阵第0行和第0列赋为0*/
for（i=0;i<=N;i++）array[i][0]=0;
for（j=0;j<=N;j++）array[0][j]=0;
for（i=1;i<=N;i++）
for（j=1;j<=N;j++）
{
if（i==j）array[i][j]=1.5*cos（i+1.2*j）;
elsearray[i][j]=sin（0.5*i+0.2*j）;
}
}
/*做矩阵A的拷贝*/
/*
求某个特征值对应的特征向量，实际上是求解一个n元线性方程组的根，
因此每求一个特征向量都要重新载入矩阵A。
*/
voidLoad_MatrixA（doubleMatrixA[N+1][N+1],doublearrayA[N+1][N+1]）
{
inti,j;
for（i=1;i<=N;i++）
for（j=1;j<=N;j++）arrayA[i][j]=MatrixA[i][j];
}
/*输出给定矩阵的所有元素，按2个元素一行输出*/
voidPrint_test（doublearray[N+1][N+1]）
{
inti,j;
for（i=1;i<=N;i++）
{
for（j=1;j<=N;j++）
{
printf（"a[%d][%d]=%.12e",i,j,array[i][j]）;
/*2个元素一行输出*/
if（j%2==0）printf（"\n"）;
}
}
}
/*输出给定矩阵A、B的积AB的所有元素，按2个元素一行输出*/
voidPrint_testAB（doublearrayA[N+1][N+1],doublearrayB[N+1][N+1]）
{
inti,j,k;
doubletemp;
for（i=1;i<=N;i++）
{
for（j=1;j<=N;j++）
{
temp=0;
for（k=1;k<=N;k++）temp+=arrayA[i][k]*arrayB[j][k];
printf（"%.12e",temp）;
/*2个元素一行输出*/
if（j%2==0）printf（"\n"）;
}
}
}
/*输出求得的所有特征值*/
voidPrint_namda（doublenamda[N+1][2]）
{
inti;
for（i=1;i<=N;i++）
{
if（namda[i][1]!
=0）
printf（"λ（%d）=%.12e+%.12ei\n",i,namda[i][0],namda[i][1]）;
else
printf（"λ（%d）=%.12e\n",i,namda[i][0]）;
}
}
/*输出矩阵A的实特征值对应的特征向量，按3个元素一行输出*/
voidPrint_vector（doublex[N+1][N+1],doublenamda[N+1][2]）
{
inti,j;
for（i=1;i<=N;i++）
{
if（namda[i][1]==0）
{
printf（"\nλ（%d）=%.12e对应的特征向量为：
\n",i,namda[i][0]）;
for（j=1;j<=N;j++）
{
printf（"%.12e",x[i][j]）;
/*3个元素一行输出*/
if（j%3==0）printf（"\n"）;
}
}
else;
}
}
/*将低于精度水平的矩阵元素清0*/
/*
拟上三角化过程后矩阵A的下三角元素的精确值应该为0，但由于计算机中存储的数与精确值都有一定误差，
因此，实际上得到的下三角元素不为0，而是一个绝对值很小的数。
为了避免这些数对后面的计算造成影响，
故要对拟上三角化后的矩阵A下三角元素进行清0。
*/
voidClear_Elements（doublearray[N+1][N+1]）
{
inti,j;
for（i=1;i<=N;i++）
for（j=1;j<=i-1;j++）{if（fabs（array[i][j]）}
/*将矩阵初Q始化为单位矩阵I*/
voidInitial_to_I（doublearray[N+1][N+1]）
{
inti,j;
for（i=0;i<=N;i++）
for（j=0;j<=N;j++）array[i][j]=0;
for（i=1;i<=N;i++）array[i][i]=1;
}
/*矩阵A的拟上三角化*/
voidNishang（doublearrayA[N+1][N+1]）
{
intr,i,j,flag;
doubled_r,c_r,h_r,t_r,temp;
doublew_r[N+1],p_r[N+1],q_r[N+1],u[N+1];
for（r=1;r<=N-2;r++）
{
/*判断（i=r+1,…，n），A[i][r]是否全为0*/
flag=0;
for（i=r+2;i<=N;i++）
{
if（arrayA[i][r]!
=0）{flag=1;break;}
else;
}
/*A[i][r]不全为0时，由A（r）求A（r+1）*/
if（flag==1）
{
/*计算d_r、c_r和h_r*/
d_r=0;
for（i=r+1;i<=N;i++）d_r+=arrayA[i][r]*arrayA[i][r];
d_r=sqrt（d_r）;
c_r=（arrayA[r+1][r]>0）?
（-d_r）:
（d_r）;
h_r=c_r*（c_r-arrayA[r+1][r]）;
/*给u_r各元素赋值*/
for（i=1;i<=N;i++）
{
if（i<=r）u[i]=0;
elseif（i==r+1）u[i]=arrayA[r+1][r]-c_r;
elseu[i]=arrayA[i][r];
}
/*计算p_r、q_r各元素*/
for（i=1;i<=N;i++）
{
temp=0;
for（j=1;j<=N;j++）temp+=arrayA[j][i]*u[j];
p_r[i]=temp/h_r;
}
for（i=1;i<=N;i++）
{
temp=0;
for（j=1;j<=N;j++）temp+=arrayA[i][j]*u[j];
q_r[i]=temp/h_r;
}
/*计算t_r的值和w_r各元素*/
temp=0;
for（i=1;i<=N;i++）
temp+=p_r[i]*u[i];
t_r=temp/h_r;
for（i=1;i<=N;i++）
w_r[i]=q_r[i]-t_r*u[i];
/*由A（r）计算A（r+1）*/
for（i=1;i<=N;i++）
for（j=1;j<=N;j++）arrayA[i][j]=arrayA[i][j]-w_r[i]*u[j]-u[i]*p_r[j];
}//endif
/*A[i][r]全为0时，A（r+1）=A（r）*/
else;
}
/*调用Clear_Elements将下三角清0*/
Clear_Elements（arrayA）;
}
/*矩阵A的QR分解*/
voidQR_Fenjie（doublearrayA[N+1][N+1],doubleQ[N+1][N+1],intm）
{
intr,i,j,flag;
doubled_r,c_r,h_r,temp;
doublew_r[N+1],p_r[N+1],u[N+1];
Initial_to_I（Q）;
for（r=1;r<=m-1;r++）
{
flag=0;
for（i=r+1;i<=m;i++）
{
if（arrayA[i][r]!
=0）{flag=1;break;}
else;
}
if（flag==1）
{
d_r=0;
for（i=r;i<=m;i++）d_r+=arrayA[i][r]*arrayA[i][r];
d_r=sqrt（d_r）;
c_r=（arrayA[r][r]>0）?
（-d_r）:
（d_r）;
h_r=c_r*（c_r-arrayA[r][r]）;
for（i=1;i<=m;i++）
{
if（i<=r-1）u[i]=0;
elseif（i==r）u[i]=arrayA[r][r]-c_r;
elseu[i]=arrayA[i][r];
}
for（i=1;i<=m;i++）
{
temp=0;
for（j=1;j<=m;j++）temp+=Q[i][j]*u[j];
w_r[i]=temp;
}
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）Q[i][j]=Q[i][j]-w_r[i]*u[j]/h_r;
for（i=1;i<=m;i++）
{
temp=0;
for（j=1;j<=m;j++）temp+=arrayA[j][i]*u[j];
p_r[i]=temp/h_r;
}
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）arrayA[i][j]=arrayA[i][j]-u[i]*p_r[j];
}//endif
else;
}
}
//****************************************************************************//
//带双步位移的QR分解求矩阵的全部特征值//
//****************************************************************************//
voidSolution_namda（doublearrayA[N+1][N+1],doublenamda[N+1][2]）
{
intk,m;
doubledet,s,t,namda1[2],namda2[2];
/*初始化存放特征值的数组元素*/
for（k=0;k<=N;k++）
for（m=0;m<=1;m++）namda[k][m]=0;
/*Loop1：
对矩阵A进行拟上三角化*/
Nishang（arrayA）;
/*Loop2：
迭代次数初始化为1，矩阵阶数初始化为N*/
k=1;
m=N;
Loop3:
if（fabs（arrayA[m][m-1]）<=Epsilon）
{
namda[m][0]=arrayA[m][m];
m--;
gotoLoop4;/*转4*/
}
else
gotoLoop5;/*转5*/
Loop4:
if（m>1）gotoLoop3;/*转3*/
elseif（m==1）
{
namda[m][0]=arrayA[m][m];
gotoLoop11;/*转11结束*/
}
elseif（m==0）gotoLoop11;/*转11结束*/
Loop5:
/*求二阶子阵的两个特征值*/
s=arrayA[m-1][m-1]+arrayA[m][m];
t=arrayA[m-1][m-1]*arrayA[m][m]-arrayA[m-1][m]*arrayA[m][m-1];
det=pow（s,2）-4*t;
if（det>=0）
{
namda1[0]=（s+sqrt（det））/2;
namda2[0]=（s-sqrt（det））/2;
}
elseif（det<0）
{
namda1[0]=s/2;
namda2[0]=namda1[0];
namda1[1]=sqrt（-1*det）/2;
namda2[1]=-1*namda1[1];
}
gotoLoop6;
Loop6:
if（m==2）
{
if（det>=0）
{
namda[m][0]=namda1[0];
namda[m-1][0]=namda2[0];
m=m-2;
}
else
{
namda[m][0]=namda1[0];
namda[m][1]=namda1[1];
namda[m-1][0]=namda2[0];
namda[m-1][1]=namda2[1];
m=m-2;
}
gotoLoop11;
}
else
gotoLoop7;
Loop7:
if（fabs（arrayA[m-1][m-2]）<=Epsilon）
{
if（det>=0）
{
namda[m][0]=namda1[0];
namda[m-1][0]=namda2[0];
m=m-2;
}
else
{
namda[m][0]=namda1[0];
namda[m][1]=namda1[1];
namda[m-1][0]=namda2[0];
namda[m-1][1]=namda2[1];
m=m-2;
}
gotoLoop4;
}
else
gotoLoop8;
Loop8:
if（k==MaxL）gotoLoop12;/*计算终止，未得到全部特征值*/
else
gotoLoop9;
Loop9:
Solution_A_Next（arrayA,m,s,t）;
gotoLoop10;
Loop10:
k++;
gotoLoop3;
Loop11:
printf（"\n求解成功:
总共迭代k=%d次\n",k）;
gotoLoop13;
Loop12:
printf（"迭代次数超出限制，算法失败！
"）;
gotoLoop13;
Loop13:
;
}
/*由A（k）计算A（k+1）*/
voidSolution_A_Next（doublearrayA[N+1][N+1],intm,doubles,doublet）
{
inti,j,k;
doubleQ[N+1][N+1],M[N+1][N+1];
doubletemp;
/*计算矩阵M_k：
M_k=AA-SA+TI*/
for（i=1;i<=m;i++）
for（j=1;j<=m;j++）
{
升级会员 