北航数值分析大作业2Word下载.docx
《北航数值分析大作业2Word下载.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业2Word下载.docx（23页珍藏版）》请在冰豆网上搜索。
=Epsilon|,则得到A的一个特征值Lamda[m][0]=A[m][m]，置m=m-1,转（3）;
否则转（4）;
（3）ifm=1,则得到A的一个特征值Lamda[m][0]=A[1][1]，转（10）;
ifm=0，则直接转（10）；
ifm>
1,则转
（2）;
（4）求二阶子阵的两个特征值，即计算二次方程的两个根，
x*x-s1[0]*x+s2[0]=0（*）
其中s1[0]=A[m-1][m-1]+A[m][m];
s2[0]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]
（该过程可用求解二次方程子程序实现，注意讨论根的情况）；
（5）ifm=2,则得到A的两个特征值，转（10）；
否则转（7）；
（6）if|A[m-1][m-2]）|<
=Epsilon,则得到A的两个特征根，置m=m-2，转（3）;
（7）ifk=L,printf“FailtogetalltheeigenvaluesofmatrixA!
”；
否则转（8）；
（8）对i,j=1,…,m，计算
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
Mk[i][j]=Mk[i][j]-s*A[i][j]+t*I[i][j]
Mk=Qk*Rk（对Mk作QR分解，可调用QR分解子程序）
Ak+1=Qk’*Ak*Qk
（9）置k=k+1，转
（2）；
（10）printf“GetalltheeigenvaluesofmatrixA!
Congratulations!
”，停止计算。
Mk的QR分解算法：
记B[][N]=Mk[][N],C[][N]=A[][N]，置flag=0（B[i][r]=0）。
对r=1到m-1执行
（1）if|B[i][r]|>
=Epsilon,flag=1，转
（2）;
（3）对i=r+1到m执行ur[i]=A[i][r];
（4）计算
vr=Br’*ur/hr
Br+1=Br-ur*vr’
pr=Crur/hr
qr=Crr/hr
Cr+1=Cr-wr*ur-ur*pr’
（5）continue.
求二元一次方程算法：
置根Lamda[][3]={0}。
对方程（*），令t1=s1[0];
t2=s2[0];
计算
delt=t1*t1-4*t2;
sqdelt=sqrt（fabs（delt））;
（2）ifdelt>
=0.0，则s1[0]=（t1+sqdelt）/2.0，s2[0]=（t1-sqdelt）/2.0;
否则转（3）;
（3）置Lamda[m][2]=1;
Lamda[m-1][2]=1;
s1[0]=t1/2.0;
s1[1]=sqdelt/2.0;
s2[0]=t1/2.0;
s2[1]=-sqdelt/2.0;
4.用列主元Gauss消去法求解实特征根对应的特征向量
（1）ifLamda[m][2]=0，转
（2）;
（2）调用GaussOptimal子程序；
（3）printf“Whentheeigenvalueis%4.12e,thecorrespondingeigenvectoris\n"
Lamda[i][0]”
GaussOptimal算法
记a1=aij，令b[N]={0}。
对i=1到10执行A1[i][i]=A1[i][i]-Lamda[p][0]
消元过程：
对于k=1到n-1执行
（1）选行号ik,使得|aik,k|=max|aik|;
（2）交换akj与aik,j（j=1到n）；
（3）对于i=k+1到n计算
mik=A1[i][k]/A1[k][k]
A1[i][j]=A1[i][j]-mik*A1[k][j]（j=k+1到n）
回代过程
（1）iffabs（A1[n][n]）<
=1e-9则iffabs（A1[n-1][n-1]）<
=1e-9，printf（"
两重根\n"
），xm=（0,0,…,0,1）,xm-1=（*,*,…,*,1,0）;
否则转（3）；
（2）对k=n-2到1执行
1）置temp=0.0;
2）对j=k+1到n执行：
temp=temp+A1[k][j]*x[p][j];
3）x[p-1][k]=-temp/A1[k][k];
（3）printf（"
\n单根"
）;
x[p][n]=1.0;
对k=n-1到1执行
1）temp=0.0;
2）对j=k+1到n执行temp=temp+A1[k][j]*x[p][j];
3）x[p][k]=-temp/A1[k][k];
（4）特征向量归一化.
二、源程序
#include<
stdio.h>
/*定义头文件*/
stdafx.h>
math.h>
#defineEpsilon1e-20/*定义精度*/
#definen10/*定义矩阵和向量的长度*/
#defineN11/*定义数组用，11时和日常表达数组相符合*/
#defineL1000/*定义最大迭代次数*/
voidHessenberg（doubleA[][N],doubleP[][N]）;
/*拟上三角化程序*/
voidQRreduce（doubleC[][N],doubleB[][N],intm）;
/*QR分解法*/
intDoublestep（doubleA[][N],doubleLamda[][3],FILE*fp）;
/*双步位移的QR分解法*/
voidroot（doubles1[],doubles2[],doubleLamda[][3],intm）;
/*求二元一次方程的根*/
voidGaussOptimal（doublea1[][N],doublex[][N],doubleLamda[][3],intp）;
/*高斯消元法*/
voidEigenvector（doubleA[][N],doublex[][N],doubleLamda[][3],FILE*fp）;
/*求解特征向量对应的特征根*/
intmain（）
{
inti=0,j=0,flag=0;
doublea[N][N]={0.0},P[N][N]={0.0},Lamda[N][3]={0.0},x[N][N]={0.0},a1[N][N]={0.0};
FILE*fp;
fp=fopen（"
QRdata.txt"
"
wt"
/*以只写方式打开QR-data文本文档*/
printf（"
*****************************************************\n"
用带双步位移的QR分解方法求解矩阵（10*10）的全部特征值\n"
/*输出标题*/
/*计算矩阵A*/
for（i=1;
i<
=n;
i++）
{
for（j=1;
j<
j++）
if（i==j）
a[i][j]=1.5*cos（i+1.2*j）;
P[i][j]=1.0;
}
else
a[i][j]=sin（0.5*i+0.2*j）;
/*输出矩阵A*/
fprintf（fp,"
\n原始矩阵A\n"
%4.12e"
a[i][j]）;
\n"
\n********************************************\n"
Hessenberg（a,P）;
/*拟上三角化*/
HessenbergisOK!
\n拟上三角化矩阵A（n-1）\n"
a1[i][j]=a[i][j];
\n拟上三角化变换矩阵P\n"
P[i][j]）;
flag=Doublestep（a,Lamda,fp）;
/*带双步位移QR分解*/
\nDoublestepisOK!
\n矩阵A（n-1）进行QR分解后所得矩阵\n"
=N;
if（flag）
\nTheeigenvaluesofmatrixA\n"
%4.12e%4.12e%4.12e\n"
Lamda[i][0],Lamda[i][1],Lamda[i][2]）;
\nFailtogettheeigenvaluesofmatrixA.Pleasecheckit!
return（0）;
Eigenvector（a1,x,Lamda,fp）;
/*求实特征值对应特征向量*/
\nEigenvectorisOK!
Getalleigenvectorsoftherealeigenvalues!
fclose（fp）;
return
（1）;
}
/*****************拟上三角化子程序*********************/
voidHessenberg（doubleA[][N],doubleP[][N]）
inti,j,r;
intflag=0;
/*flag=0全为零，flag=1不全为零*/
doublecr,dr,hr,tr,temp;
doublepr[N]={0.0},qr[N]={0.0},ur[N]={0.0},wr[N]={0.0},p[N]={0.0};
for（r=1;
r<
=n-2;
r++）
cr=0.0;
dr=0.0;
hr=0.0;
tr=0.0;
temp=0.0;
pr[i]=0;
qr[i]=0;
ur[i]=0;
wr[i]=0;
p[i]=0;
for（i=r+2;
i++）/*判断第r列是否全为，若有一个不是则flag=1*/
if（fabs（A[i][r]）>
=Epsilon）
flag=1;
break;
if（flag==1）/*计算dr，cr，hr；
给ur赋值*/
for（i=r+1;
temp=temp+A[i][r]*A[i][r];
dr=sqrt（temp）;
if（A[r+1][r]>
Epsilon）
cr=-dr;
cr=dr;
hr=cr*（cr-A[r+1][r]）;
ur[i]=A[i][r];
ur[r+1]=ur[r+1]-cr;
i++）/*计算pr，qr，tr，wr和A*/
pr[i]=pr[i]+A[j][i]*ur[j];
qr[i]=qr[i]+A[i][j]*ur[j];
p[i]=p[i]+P[i][j]*ur[j];
pr[i]=pr[i]/hr;
qr[i]=qr[i]/hr;
p[i]=p[i]/hr;
tr=tr+pr[i]*ur[i];
tr=tr/hr;
wr[i]=qr[i]-tr*ur[i];
if（（j==r）&
&
（i>
j+1））/*ur前r个元素为零*/
A[i][j]=0.0;
A[i][j]=A[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
P[i][j]=P[i][j]-p[i]*ur[j];
/*****************QR分解子程序*********************/
voidQRreduce（doubleC[][N],doubleB[][N],intm）
doublepr[N],qr[N],ur[N],vr[N],wr[N];
=m-1;
=m;
pr[i]=0.0;
qr[i]=0.0;
ur[i]=0.0;
vr[i]=0.0;
wr[i]=0.0;
if（fabs（B[i][r]）>
if（flag==1）
for（i=r;
temp=temp+B[i][r]*B[i][r];
if（B[r][r]>
hr=cr*（cr-B[r][r]）;
ur[i]=B[i][r];
ur[r]=ur[r]-cr;
vr[i]=vr[i]+B[j][i]*ur[j];
pr[i]=pr[i]+C[j][i]*ur[j];
qr[i]=qr[i]+C[i][j]*ur[j];
vr[i]=vr[i]/hr;
C[i][j]=C[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
B[i][j]=B[i][j]-ur[i]*vr[j];
/****************带双步位移QR分解子程序****************/
intDoublestep（doubleA[][N],doubleLamda[][3],FILE*fp）
inti,j,k,m,l;
intflag=1;
doubles1[2]={0.0},s2[2]={0.0},s,t;
doubleMk[N][N]={0.0},I[N][N]={0.0};
k=0;
m=n;
I[i][i]=1.0;
while
（1）
Loop1:
if（fabs（A[m][m-1]）<
Lamda[m][0]=A[m][m];
m=m-1;
Loop2:
if（m==1）
\n1Getalltheeigenvalues!
Congratulations!
\n最大迭代次数k：
%4d\n"
k）;
return（flag）;
elseif（m==0）
\n0Getalltheeigenvalues!
gotoLoop1;
s1[0]=A[m-1][m-1]+A[m][m];
s1[1]=0;
s2[0]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
s2[1]=0;
root（s1,s2,Lamda,m）;
if（m==2）
for（i=0;
2;
Lamda[m][i]=s1[i];
Lamda[m-1][i]=s2[i];
\n2Getalltheeigenvalues!
if（fabs（A[m-1][m-2]）<
m=m-2;
gotoLoop2;
if（k==L）
\nLFailtogetalltheeigenvalues!
Pleasecheck!
flag=0;
/*计算中间过度矩阵Mk*/
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m-1][m]*A
升级会员 