数值分析实习第二题.docx
《数值分析实习第二题.docx》由会员分享，可在线阅读，更多相关《数值分析实习第二题.docx（31页珍藏版）》请在冰豆网上搜索。
数值分析实习第二题
数值分析实习第二题
一、算法设计方案
1.采用C语言进行算法设计输出。
2.由于矩阵A的为10×10的矩阵，所以采用双步位移的QR分解法求矩阵A的全部特征值。
3.先对矩阵A进行初始化赋值。
4.先采用Householder矩阵将矩阵A拟上三角化，得到矩阵A（n-1）。
5.为了完成题目要求，减少运算次数，此时运行独立函数QR_UP_HE（）另外声明3个10×10局部的二维数组变量Q，R，RQ对A（n-1）进行一次QR分解，得到相对应的矩阵Q，R和乘积矩阵RQ，用于输出。
6.之后开始用双步位移法求解矩阵A的特征值。
在求解过程中，遇到解一元二次方程组，采用数组S_Re[2]来存根的实数部分，用数组S_Im[2]来存根的虚数部分，用公式
来求解方程
。
QR分解时候选用精度ε=10-12，不限制最大迭代次数。
7.将双步位移的判断逻辑进行合理设定如下：
1）变量m从9（采用10维数组）开始循环运算，当m小于零的时候即认为QR分解完成，结束循环。
2）先判断m是否为零，若为零，那么A[0][0]即为最后一个特征值，结束循环。
否则继续。
3）之后判断A[m][m-1]的绝对值是否小于精度要求，若小于，则认为A[m][m]即为其中一个特征值，运算维数降为m–1，跳过这次循环。
否则继续。
4）此时求二阶子阵
的两个特征值s1和s2。
5）此时如果m=1，那么这就是A剩下的最后两个特征值，结束循环。
否则继续。
6）如果A[m-1][m-2]的绝对值小于精度要求，那么即可判断s1和s2是A的其中两个特征值，运算维数降为m–2，跳过这次循环。
否则继续。
7）上述判断完之后，即可开始用双步位移的QR分解法循环求A得全部特征值了。
8.全部特征值求得后，将其打印出来，之后开始求A的相应于实特征值的特征向量。
特征向量采用解方程组（λI-A）X=0的方法进行，先用Gauss列主元的方法，将A上三角化。
由于b为0，所以不用记录A对角化过程中的行交换的次序。
9.在反解X的过程中，需先判断A的行是否全为0（即小于一个精度值，程序中精度值定为ε=10-9，取太小会导致程序判断出错，从而得出错误的特征值）。
若最后一行不为零，则对特征向量X的最后一个值赋1，即X=[0,0,0,…,0,1]，若倒数第二行也不为零，则对特征向量X的倒数第二行的值也赋1，以此类推。
10.之后回代X，得到其中的一个特征向量X。
然后将其归一化，即可输出。
二、运算结果
1.矩阵经过拟上三角化后得到的矩阵
2.进行QR分解后所得到的矩阵Q、R和乘积矩阵RQ
3.矩阵的全部特征值
4.
相应于特征值的特征向量
三、源程序
#include
#include
#include
#defineEPS1.0e-12
#defineN10
doubleA[N][N];
doubleLambda_Re[N];
doubleLambda_Im[N];
voidINI_A（void）
//矩阵A进行赋值
{
chari,j;
for（i=0;i{
for（j=0;j{
if（j==i）
{
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））;
}
}
}
}
voidVECT_MULTI（charMax_N,doubleVector[],doubleVect_Ret[],doublek）
//向量数乘
{
charn;
for（n=0;n{
Vect_Ret[n]=Vector[n]*k;
}
}
doubleVECT_NORM（charRow,charStr,charEnd,doubleMat[][N]）
//求矩阵A中列向量的模
{
charn;
doublenorm;
for（n=Str,norm=0;n{
norm+=Mat[n][Row]*Mat[n][Row];
}
return（sqrt（norm））;
}
doubleVECT_NORM2（charMax_N,doubleVector_2[]）
//向量2范数
{
intn;
doublenorm_2;
for（n=0,norm_2=0;n{
norm_2+=（Vector_2[n]*Vector_2[n]）;
}
return（sqrt（norm_2））;
}
charSP_SGN（doubles）
//取符号
{
if（s>0）
{
return
（1）;
}
else
{
return（-1）;
}
}
doubleVECT_T_VEC（charMax_N,doubleVector_T[],doubleVector[]）
//向量转置乘向量
{
charn;
doublepower_norm;
for（n=0,power_norm=0;n{
power_norm+=Vector_T[n]*Vector[n];
}
return（power_norm）;
}
voidVECT_PLUS（charMax_N,doubleVector_1[],doubleVector_2[],doubleVector_Ret[]）
//向量加法
{
chari;
for（i=0;i{
Vector_Ret[i]=Vector_1[i]+Vector_2[i];
}
}
voidMAT_VEC（charMax_N,doubleMatrix[][N],doubleVect[],doubleVect_Ret[]）
//矩阵乘向量
{
chari,j;
doublecache[N]={0};
for（i=0;i{
for（j=0;j{
cache[i]+=Matrix[i][j]*Vect[j];
}
}
for（i=0;i{
Vect_Ret[i]=cache[i];
}
}
voidMAT_T_VEC（charMax_N,doubleMatrix[][N],doubleVect[],doubleVect_Ret[]）
//矩阵转置乘向量
{
chari,j;
doublecache[N]={0};
for（i=0;i{
for（j=0;j{
cache[i]+=Matrix[j][i]*Vect[j];
}
}
for（i=0;i{
Vect_Ret[i]=cache[i];
}
}
voidVECT_EVA（charMax_N,doubleVect_new[],doubleVect_bef[]）
//向量赋值
{
chari;
for（i=0;i{
Vect_new[i]=Vect_bef[i];
}
}
voidMAT_MULTI（charMax_N,doubleMat_A[][N],doubleMat_B[][N],doubleRet_Mat[][N]）
//矩阵乘法
{
chari,j,k;
for（i=0;i{
for（j=0;j{
Ret_Mat[i][j]=0;
for（k=0;k{
Ret_Mat[i][j]+=Mat_A[i][k]*Mat_B[k][j];
}
}
}
}
voidSOLVE_EQU（doubleCoe[2],doubleSol_Re[2],doubleSol_Im[2]）
//解一元二次方程
{
doubleDelta;
Delta=Coe[0]*Coe[0]/4-Coe[1];
if（Delta>=0）
{
Sol_Re[0]=（-1）*Coe[0]/2+sqrt（Delta）;
Sol_Re[1]=（-1）*Coe[0]/2-sqrt（Delta）;
Sol_Im[0]=0;
Sol_Im[1]=0;
}
else
{
Sol_Re[0]=（-1）*Coe[0]/2;
Sol_Re[1]=（-1）*Coe[0]/2;
Sol_Im[0]=sqrt（（-1）*Delta）;
Sol_Im[1]=（-1）*sqrt（（-1）*Delta）;
}
}
voidEVA_MAT（charMax_N,doubleBef[][N],doubleRet[][N]）
//矩阵相等赋值
{
chari,j;
for（i=0;i{
for（j=0;j{
Ret[i][j]=Bef[i][j];
}
}
}
voidEXCHANGE（double*No_1,double*No_2）
//数据交换
{
doubletemp;
temp=*No_1;
*No_1=*No_2;
*No_2=temp;
}
voidUPPER_HESSENBERG（void）
//矩阵拟上三角化
{
charr,i,k,j;
doubled,c,h;
doublep[N],q[N],w[N],u[N];
for（r=0;r{
for（i=r+2;i{
if（A[i][r]!
=0）
{
d=VECT_NORM（r,r+1,N,A）;
c=（-1）*SP_SGN（A[r+1][r]）*d;
h=c*c-c*A[r+1][r];
for（k=0;k{
if（k>r+1）
{
u[k]=A[k][r];
}
elseif（k==r+1）
{
u[k]=A[k][r]-c;
}
else
{
u[k]=0;
}
}
VECT_MULTI（N,u,w,1/h）;
MAT_T_VEC（N,A,w,p）;
MAT_VEC（N,A,w,q）;
VECT_MULTI（N,u,w,-VECT_T_VEC（N,p,w））;
VECT_PLUS（N,q,w,w）;
for（k=0;k{
for（j=0;j{
A[k][j]-=（w[k]*u[j]+u[k]*p[j]）;
}
}
break;
}
}
}
}
voidQR_DISSOLUTION（charMax_N,doubleMat[][N]）
//矩阵双步位移QR分解
{
chari,j,r,k;
doubled,c,h;
doubleu[Max_N],w[Max_N],p[Max_N],q[Max_N];
for（r=0;r{
for（i=r+1;i{
if（Mat[i][r]!
=0）
{
d=VECT_NORM（r,r,Max_N,Mat）;
c=（-1）*SP_SGN（Mat[r][r]）*d;
h=c*c-c*Mat[r][r];
for（k=0;k{
if（k>r）
{
u[k]=Mat[k][r];
}
elseif（k==r）
{
u[k]=Mat[k][r]-c;
}
else
{
u[k]=0;
}
}
VECT_MULTI（Max_N,u,w,1/h）;
MAT_T_VEC（Max_N,Mat,w,p）;
for（k=0;k{
for（j=0;j{
Mat[k][j]-=（u[k]*p[j]）;
}
}
MAT_T_VEC（Max_N,A,w,p）;
MAT_VEC（Max_N,A,w,q）;
VECT_MULTI（Max_N,u,w,-VECT_T_VEC（Max_N,p,w））;
VECT_PLUS（Max_N,q,w,w）;
for（k=0;k{
for（j=0;j{
A[k][j]-=（w[k]*u[j]+u[k]*p[j]）;
}
}
break;
}
}
}
}
voidDOUBLE_DISPLACE（void）
//双步位移法求特征值
{
chari,j,m;
doubleS_Re[2],S_Im[2],D_Coe[2];
doubleMk[N][N];
for（m=N-1;m>=0;）
{
if（m==0）
{
Lambda_Re[0]=A[0][0];
Lambda_Im[m]=0;
break;
}
if（fabs（A[m][m-1]）<=EPS）
{
Lambda_Re[m]=A[m][m];
Lambda_Im[m]=0;
m--;
continue;
}
D_Coe[0]=（-1）*（A[m-1][m-1]+A[m][m]）;
D_Coe[1]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
SOLVE_EQU（D_Coe,S_Re,S_Im）;
if（m==1）
{
Lambda_Re[0]=S_Re[0];
Lambda_Re[1]=S_Re[1];
Lambda_Im[0]=S_Im[0];
Lambda_Im[1]=S_Im[1];
break;
}
if（fabs（A[m-1][m-2]）<=EPS）
{
Lambda_Re[m-1]=S_Re[0];
Lambda_Re[m]=S_Re[1];
Lambda_Im[m-1]=S_Im[0];
Lambda_Im[m]=S_Im[1];
m=m-2;
continue;
}
MAT_MULTI（（m+1）,A,A,Mk）;
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{
Mk[i][j]=Mk[i][j]+D_Coe[0]*A[i][j];
if（i==j）
{
Mk[i][j]=Mk[i][j]+D_Coe[1];
}
}
}
QR_DISSOLUTION（（m+1）,Mk）;
}
}
voidQR_UP_HE（void）
//拟上三角矩阵QR分解
{
doubleR[N][N]={0},Q[N][N]={0},RQ[N][N];
chari,j,r,k;
doublep[N],w[N],u[N];
doubled,c,h;
FILE*fp;
EVA_MAT（N,A,R）;
for（i=0;i{
Q[i][i]=1;
}
for（r=0;r{
for（i=r+1;i{
if（R[i][r]!
=0）
{
d=VECT_NORM（r,r,N,R）;
c=（-1）*SP_SGN（R[r][r]）*d;
h=c*c-c*R[r][r];
for（k=0;k{
if（k>r）
{
u[k]=R[k][r];
}
elseif（k==r）
{
u[k]=R[k][r]-c;
}
else
{
u[k]=0;
}
}
VECT_MULTI（N,u,p,1/h）;
MAT_VEC（N,Q,u,w）;
for（k=0;k{
for（j=0;j{
Q[k][j]-=（w[k]*p[j]）;
}
}
MAT_T_VEC（N,R,p,w）;
for（k=0;k{
for（j=0;j{
R[k][j]-=（u[k]*w[j]）;
}
}
break;
}
}
}
if（（fp=fopen（"A（n-1）-QR","w"））==NULL）
//输出QR分解结果
{
printf（"cannotopenfile\n"）;
exit（0）;
}
fprintf（fp,"SY1103104胡延勍数值分析第二题\n"）;
fprintf（fp,"拟上三角化得到矩阵A（n-1）进行QR分解后得到：
\n"）;
fprintf（fp,"Q矩阵为：
\n"）;
for（i=0;i{
for（j=0;j{
fprintf（fp,"%.13e",Q[i][j]）;
}
fprintf（fp,"\n\n"）;
}
fprintf（fp,"R矩阵为：
\n"）;
for（i=0;i{
for（j=0;j{
fprintf（fp,"%.13e",R[i][j]）;
}
fprintf（fp,"\n\n"）;
}
MAT_MULTI（N,R,Q,RQ）;
fprintf（fp,"乘积矩阵RQ：
\n"）;
for（i=0;i{
for（j=0;j{
fprintf（fp,"%.13e",RQ[i][j]）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
}
voidGAUSS_PIVOT（doubleLamb,doublex[]）
//Gauss选主元求特征向量
{
chari,k,j,line;
doublem;
for（i=0;i{
for（k=0;k{
if（i==k）
{
A[i][k]=Lamb-A[i][k];
}
else
{
A[i][k]=（-1）*A[i][k];
}
}
}
for（k=0;k{
for（i=k,line=k;i
升级会员 