数值分析计算实习二Word格式.docx
《数值分析计算实习二Word格式.docx》由会员分享，可在线阅读，更多相关《数值分析计算实习二Word格式.docx（19页珍藏版）》请在冰豆网上搜索。
先判断迭代矩阵右下角是否为二阶分块矩阵，若是，再求解其一对特征值。
d）根据求得的实特征值计算对应的特征向量。
因为特征值已知，只需要求解以下线性方程组的一非零解即可：
为了方便起见，直接利用QR分解函数分解矩阵B，则原方程组等价于
然后利用直接法解出方程组即可。
其中矩阵R右下角元素必为0，则使对应行未知数为1，然后回代即可。
当矩阵R对应特征值的重数为2时，在矩阵R中除了右下角元素为0外，必有另一个对角元素为0，则令最后行对应未知数为0，该行对应未知数为1.当重数大于2时依次类推，即可求出对应于特定重数特征值的特征向量。
结果：
1、
拟上三角矩阵
如下。
2、迭代最终矩阵如下，迭代次数为14次。
3、矩阵所有特征值如下。
从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。
4、对应于实特征值的特征向量如下
问题发现与讨论
1、双步位移QR迭代过程中迭代矩阵
始终为拟上三角矩阵。
2、双步位移QR方法与基本QR方法比较
基本QR方法迭代过程中迭代矩阵
下图为基本QR方法得到的分块上三角矩阵：
下图为基本QR方法得到的特征值：
在基本QR迭代中发现，当矩阵收敛到分块上三角矩阵后，再次迭代它仍旧为分块上三角矩阵，然而矩阵中的元素改变了，不过块的阶数和位置不会改变。
基本QR迭代收敛速度为545次。
可以看出虽然两者收敛到的分块上三角矩阵不同，但得到的最终特征值是一样的。
也从而可初步验证所求得得特征值的正确性。
收敛速度上基本QR方法明显慢于双步位移QR方法。
源程序：
#include"
stdafx.h"
stdio.h"
math.h"
stdlib.h"
#defineN10
#defineEp1e-12
int_tmain（intargc,_TCHAR*argv[]）
{
intar0（doublea[N][N],ints,intn）;
//判断矩阵多少位后是0
intsgn（doublex）;
//0为1的符号函数
voidatra（doublea[N][N],intn）;
//矩阵转置
voidaadd（doublea[N],doubleb[N],doublec[N],intn,intco）;
//向量加法
voidaaadd（doublea[N][N],doubleb[N][N],doublec[N][N],intn,intco）;
//矩阵加法
voidnamul（doublea[N],doublex,doubleb[N],intn）;
//向量数乘
voidnaamul（doublea[N][N],doublex,doubleb[N][N],intn）;
//矩阵数乘
voidamul（doublea[N],doubleb[N],doublec[N][N],intn）;
//向量乘法
doubleaneiji（doublea[N],doubleb[N],intn）;
//向量内积
voidaamul（doublea[N],doubleb[N][N],doublec[N],intn）;
//矩阵向量乘法
voidaaamul（doublea[N][N],doubleb[N][N],doublec[N][N],intn）;
//矩阵矩阵乘法
voidsol2（doublea11,doublea12,doublea21,doublea22,double（*p1）[2],double（*p2）[2]）;
//计算二阶矩阵的特征值
voidqrs（doublea[N][N],doubleq[N][N]）;
voidGauss（doublea[N][N],doubley[N][N],int*p）;
doubleA[N][N];
doublenmd[N][2]={0};
for（inti=0;
i<
N;
i++）//对原矩阵赋值
{
for（intj=0;
j<
j++）
{
if（i==j）A[i][j]=1.52*cos（i+1+1.2*（j+1））;
elseA[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
}
}
for（intr=0;
r<
N-2;
r++）//将原矩阵拟上三角化
if（ar0（A,r+2,r）==0）continue;
else
doubled=0;
for（inti=r+1;
i++）d=d+pow（A[i][r],2）;
doubledr=sqrt（d）;
doublecr=-sgn（A[r+1][r]）*dr;
doublehr=pow（cr,2）-cr*A[r+1][r];
doubleur[N]={0};
ur[r+1]=A[r+1][r]-cr;
for（inti=r+2;
i++）ur[i]=A[i][r];
doublepr[N];
atra（A,N）;
aamul（ur,A,pr,N）;
namul（pr,1/hr,pr,N）;
doubleqr[N];
aamul（ur,A,qr,N）;
namul（qr,1/hr,qr,N）;
doubletr;
tr=aneiji（pr,ur,N）/hr;
doublewr[N];
namul（ur,-tr,wr,N）;
aadd（qr,wr,wr,N,1）;
doubleB[N][N];
amul（wr,ur,B,N）;
aaadd（A,B,A,N,-1）;
amul（ur,pr,B,N）;
//求一步操作后的矩阵Ar
i++）
printf（"
%.11lf"
A[i][j]）;
if（j==N-1）printf（"
\n"
）;
}
printf（"
//打印拟上三角矩阵Ar
doubleAqr[N][N]={0};
aaadd（Aqr,A,Aqr,N,1）;
doubleqqr[N][N]={0};
10000;
doubleAqr1[N][N]={0};
aaadd（Aqr,Aqr1,Aqr1,N,1）;
qrs（Aqr,qqr）;
doubleB[N][N]={0};
aaamul（Aqr,qqr,B,N）;
for（ints=0;
s<
s++）
for（intk=0;
k<
k++）
Aqr[s][k]=B[s][k];
ints=0;
for（intt=0;
t<
t++）
if（fabs（Aqr[t+1][t]）<
Ep）s=s+1;
if（s==N-2）break;
Aqr[i][j]）;
//打印A（n-1）进行qr迭代后的最终结果
intm=N-1;
intk=0;
while（m>
-1）
if（（fabs（A[m][m-1]）<
Ep）||（m==0））
nmd[m][0]=A[m][m];
nmd[m][1]=0;
m=m-1;
continue;
elseif（（fabs（A[m-1][m-2]）<
Ep）||（m==1））
{
sol2（A[m-1][m-1],A[m-1][m],A[m][m-1],A[m][m],（nmd+m）,（nmd+（m-1）））;
//求解二阶矩阵的特征值
m=m-2;
continue;
}
doubles=A[m-1][m-1]+A[m][m];
doublet=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];
doubleM[N][N]={0};
doubleD[N][N]={0};
aaamul（A,A,M,m+1）;
naamul（A,s,D,m+1）;
aaadd（M,D,M,m+1,-1）;
for（inti=0;
m+1;
M[i][i]=M[i][i]+t;
for（intr=0;
m;
r++）//M的QR分解和Ar的计算
if（ar0（A,r+1,r）==0）continue;
else
{
doubled=0;
for（inti=r;
m+1;
i++）
d=d+pow（M[i][r],2）;
doubledr=sqrt（d）;
doublecr=-sgn（M[r][r]）*dr;
doublehr=pow（cr,2）-cr*M[r][r];
doubleur[N]={0};
ur[r]=M[r][r]-cr;
for（inti=r+1;
i++）ur[i]=M[i][r];
doublevr[N];
atra（M,m+1）;
aamul（ur,M,vr,m+1）;
namul（vr,1/hr,vr,m+1）;
amul（ur,vr,D,m+1）;
aaadd（M,D,M,m+1,-1）;
doublepr[N];
atra（A,m+1）;
aamul（ur,A,pr,m+1）;
namul（pr,1/hr,pr,m+1）;
doubleqr[N];
aamul（ur,A,qr,m+1）;
namul（qr,1/hr,qr,m+1）;
doubletr;
tr=aneiji（pr,ur,m+1）/hr;
doublewr[N];
namul（ur,-tr,wr,m+1）;
aadd（qr,wr,wr,m+1,1）;
amul（wr,ur,D,m+1）;
aaadd（A,D,A,m+1,-1）;
amul（ur,pr,D,m+1）;
}
k++;
if（k>
10000）//算法失败判定
算法失败\n"
break;
printf（"
%1.11e"
nmd[i][0]）;
%1.11e\n"
nmd[i][1]）;
if（fabs（nmd[i][1]）>
Ep）continue;
for（intj=i+1;
if（fabs（nmd[j][1]）>
elseif（fabs（nmd[j][0]-nmd[i][0]）<
Ep）nmd[j][1]=1;
for（intk=0;
if（fabs（nmd[k][1]）>
for（inti=0;
i++）//求矩阵A-I*nmd
for（intj=0;
if（i==j）A[i][j]=1.52*cos（i+1+1.2*（j+1））-nmd[k][0];
elseA[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
doubleq[N][N];
doubley[N][N];
intp[1]={0};
qrs（A,q）;
//对A-I*nmd进行QR分解得到新矩阵R
Gauss（A,y,p）;
//齐次不满秩的高斯消去法,求非0根
对应于λ%d实特征值的特征向量为\n"
k+1）;
p[0];
printf（"
%1.12e\n"
y[i][j]）;
return0;
}
intar0（doublea[N][N],ints,intn）
for（inti=s;
if（fabs（a[i][n]）>
Ep）
return1;
}
intsgn（doublex）
if（x>
0）return1;
elsereturn-1;
voidatra（doublea[N][N],intn）
doublex;
n;
for（intj=i+1;
x=a[i][j];
a[i][j]=a[j][i];
a[j][i]=x;
voidaadd（doublea[N],doubleb[N],doublec[N],intn,intco）
if（co==1）
i++）c[i]=a[i]+b[i];
else
if（co==-1）
i++）c[i]=a[i]-b[i];
向量加减选择错误\n"
voidaaadd（doublea[N][N],doubleb[N][N],doublec[N][N],intn,intco）
{
if（co==1）
c[i][j]=a[i][j]+b[i][j];
for（intj=0;
c[i][j]=a[i][j]-b[i][j];
矩阵加减选择错误\n"
voidnamul（doublea[N],doublex,doubleb[N],intn）
b[i]=x*a[i];
voidnaamul（doublea[N][N],doublex,doubleb[N][N],intn）
b[i][j]=x*a[i][j];
voidamul（doublea[N],doubleb[N],doublec[N][N],intn）
c[i][j]=a[i]*b[j];
doubleaneiji（doublea[N],doubleb[N],intn）
doublet=0;
t=t+a[i]*b[i];
returnt;
voidaamul（doublea[N],doubleb[N][N],doublec[N],intn）
c[i]=0;
c[i]=c[i]+b[i][j]*a[j];
voidaaamul（doublea[N][N],doubleb[N][N],doublec[N][N],intn）
c[i][j]=0;
for（ints=0;
c[i][j]=c[i][j]+a[i][s]*b[s][j];
}
voidsol2（doublea11,doublea12,doublea21,doublea22,double（*p1）[2],double（*p2）[2]）
doubledeta;
deta=a11*a11+a22*a22-2*a11*a22+4*a12*a21;
if（deta<
0）
（*p1）[0]=（a11+a22）/2;
（*p2）[0]=（a11+a22）/2;
（*p1）[1]=sqrt（0-deta）/2;
（*p2）[1]=0-sqrt（0-deta）/2;
（*p1）[0]=（a11+a22）/2+sqrt（deta）/2;
（*p1）[1]=0;
（*p2）[0]=（a11+a22）/2-sqrt（deta）/2;
（*p2）[1]=0;
voidqrs（doublea[N][N],doubleq[N][N]）
if（i==j）q[i][j]=1;
elseq[i][j]=0;
N-1;
r++）
if（ar0（a,r+1,r）==0）continue;
doubled=0,dr;
doublecr;
doublehr;
for（inti=r
升级会员 