北航数值分析A第二题.docx
《北航数值分析A第二题.docx》由会员分享，可在线阅读，更多相关《北航数值分析A第二题.docx（22页珍藏版）》请在冰豆网上搜索。
北航数值分析A第二题
《数值分析A》计算实习第二题
一、设计方案
1.定义矩阵
。
由函数A（double（*a）[N]）定义矩阵
中的元素，方便初始化。
2.对矩阵
进行拟上三角化，得到
。
由函数AN（double（*a）[N]）对
进行拟上三角化，减少计算量。
3.对
采用带双步位移的
分解法求解特征值。
在每次
分解前，如果可以直接得到矩阵
的一个特征值或者一对特征值，则对矩阵
降阶。
函数SQR（double（*a）[N],double（*lc）[2]）代表带双步位移的
分解，其中lc[N][0]存储特征值的实部，lc[N][1]存储特征值的虚部。
函数HAH（double（*a）[N],double*u,doubleh,intn）通过参数u[N]、h的传递对矩阵
左乘、右乘Householder矩阵。
4.采用列主元Gauss消去法求实特征值对应的特征向量。
对
初始化后，采用列主元素Gauss消去法求
的解，消元过程中，矩阵
的最后一行元素为0，令向量
的最后一个元素为1，回代可求得对应于特征值
的特征向量
。
所用函数为G（double*x,doublelc）。
5.对
初始化，研究拟上三角化对一般
分解法求得收敛的分块上三角阵所需
分解次数的影响，对比分块上三角阵一阶、二阶对角块给出的特征值与带双步位移的
分解法求解结果的一致性。
所用函数为QR（double（*a）[N],double（*q）[N],double（*r）[N]），采用一阶对角块的元素、二阶对角块的迹组成的向量
，通过判断
，确定
分解的收敛性。
二、源代码
#include
#include
#defineN10//矩阵阶数
#defineK50000//最大迭代次数
#defineE1e-12//精度
voidA（double（*a）[N]）//初始化矩阵A
{
inti,j;
for（i=0;ifor（j=0;ja[i][j]=i!
=j?
sin（0.5*i+0.2*j+0.7）:
1.5*cos（i+1.2*j+2.2）;
}
voidHAH（double（*a）[N],double*u,doubleh,intn）//对矩阵A左乘右乘Householder矩阵
{
doublep[N],q[N],t,w[N],s1,s2;
inti,j;
for（i=0;i{s1=0.0;
s2=0.0;
for（j=0;j{
s1+=a[j][i]*u[j]/h;
s2+=a[i][j]*u[j]/h;
}
p[i]=s1;
q[i]=s2;
}
t=0.0;
for（i=0;it+=p[i]*u[i];
t/=h;
for（i=0;iw[i]=q[i]-t*u[i];
for（i=0;ifor（j=0;ja[i][j]-=（w[i]*u[j]+u[i]*p[j]）;
}
voidAN（double（*a）[N]）//拟上三角化
{
doubleu[N],d,c,h;
intr,i,t;
for（r=1;r<=N-2;r++）
{
t=0;
for（i=r+2;i<=N;i++）
if（fabs（a[i-1][r-1]）t+=1;
if（t{
d=0.0;
for（i=r+1;i<=N;i++）
d+=a[i-1][r-1]*a[i-1][r-1];
d=sqrt（d）;
c=a[r][r-1]<=0?
d:
-d;
h=c*（c-a[r][r-1]）;
for（i=1;i<=N;i++）
u[i-1]=i0:
a[i-1][r-1];
u[r]=a[r][r-1]-c;
HAH（a,u,h,N）;
}
}
}
voidSQR（double（*a）[N],double（*lc）[2]）
{
doubles,t,s1r,s1i,s2r,s2i,b,c,d,m[N][N],sum,h,u[N],v[N];
intk=0,n=N,i,j,r;
while（n>1）
{
if（fabs（a[n-1][n-2]）{
lc[N-n][0]=a[n-1][n-1];
lc[N-n][1]=0.0;
printf（"QR分解次数为%d，矩阵阶数为%d时，由最后1行倒数第2个元素为0，求得特征值\nλ%d=%.12e\n",k,n,N-n+1,lc[N-n][0]）;
printf（"此时矩阵A的%d阶顺序主子式为\n",n）;
for（i=0;i{
for（j=0;j<（nn:
N/2）;j++）
printf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
if（n>N/2）
{
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
}
n-=1;
}
else
{
//解特征值s1，s2
b=-a[n-2][n-2]-a[n-1][n-1];
c=a[n-2][n-2]*a[n-1][n-1]-a[n-1][n-2]*a[n-2][n-1];
d=b*b-4*c;
if（d<0）
{
s1r=-b/2.0;
s2r=-b/2.0;
s1i=sqrt（-d）/2.0;
s2i=-sqrt（-d）/2.0;
}
else
{
s1r=（-b+sqrt（d））/2.0;
s2r=（-b-sqrt（d））/2.0;
s1i=0.0;
s2i=0.0;
}
if（n==2||fabs（a[n-2][n-3]）{
lc[N-n][0]=s1r;
lc[N-n][1]=s1i;
lc[N-n+1][0]=s2r;
lc[N-n+1][1]=s2i;
if（n==2）
{
if（fabs（s1i）printf（"QR分解次数为%d，矩阵阶数为%d时，直接求得特征值\nλ%d=%.12e\nλ%d=%.12e\n",k,n,N-n+1,lc[N-n][0],N-n+2,lc[N-n+1][0]）;
else
printf（"QR分解次数为%d，矩阵阶数为%d时，直接求得特征值\nλ%d=%.12e+%.12ei\nλ%d=%.12e+%.12ei\n",k,n,N-n+1,lc[N-n][0],lc[N-n][1],N-n+2,lc[N-n+1][0],lc[N-n+1][1]）;
}
else
{
if（fabs（s1i）printf（"QR分解次数为%d，矩阵阶数为%d时，由倒数第2行倒数第3个元素为0，求得特征值\nλ%d=%.12e\nλ%d=%.12e\n",k,n,N-n+1,lc[N-n][0],N-n+2,lc[N-n+1][0]）;
else
printf（"QR分解次数为%d，矩阵阶数为%d时，由倒数第2行倒数第3个元素为0，求得特征值\nλ%d=%.12e+%.12ei\nλ%d=%.12e+%.12ei\n",k,n,N-n+1,lc[N-n][0],lc[N-n][1],N-n+2,lc[N-n+1][0],lc[N-n+1][1]）;
}
printf（"此时矩阵A的%d阶顺序主子式为\n",n）;
for（i=0;i{
for（j=0;j<（nn:
N/2）;j++）
printf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
if（n>N/2）
{
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
}
n-=2;
}
elseif（k{
s=a[n-2][n-2]+a[n-1][n-1];
t=a[n-2][n-2]*a[n-1][n-1]-a[n-1][n-2]*a[n-2][n-1];
for（i=0;ifor（j=0;j{
sum=0.0;
for（r=0;rsum+=a[i][r]*a[r][j];
m[i][j]=sum-s*a[i][j]+t*（i==j?
1:
0）;
}
for（r=1;r{
sum=0.0;
for（i=r+1;i<=n;i++）
if（fabs（m[i-1][r-1]）sum+=1;
if（sum{
d=0.0;
for（i=r;i<=n;i++）
d+=m[i-1][r-1]*m[i-1][r-1];
d=sqrt（d）;
c=m[r-1][r-1]<=0?
d:
-d;
h=c*（c-m[r-1][r-1]）;
for（i=1;i<=n;i++）
u[i-1]=i0:
m[i-1][r-1];
u[r-1]=m[r-1][r-1]-c;
for（i=0;i{
sum=0.0;
for（j=0;jsum+=m[j][i]*u[j];
v[i]=sum/h;
}
for（i=0;ifor（j=0;jm[i][j]=m[i][j]-u[i]*v[j];
HAH（a,u,h,n）;
}
}
k+=1;
}
else
{
printf（"未求出全部特征值\n"）;
break;
}
}
}
if（n==1）
{
lc[N-n][0]=a[0][0];
lc[N-n][1]=0.0;
printf（"QR分解次数为%d，矩阵阶数为%d时，直接求得特征值\nλ%d=%.12e\n",k,n,N-n+1,lc[N-n][0]）;
printf（"此时矩阵A的1阶顺序主子式为\n%20.12e\n\n",a[0][0]）;
}
}
voidG（double*x,doublelc）//列主元素gauss消去法求特征向量
{
doublea[N][N],t,m;
inti,j,k,ik;
A（a）;
for（i=0;ifor（j=0;ja[i][j]-=i==j?
lc:
0;
for（k=1;k{
t=0.0;
for（i=k;i<=N;i++）
if（fabs（a[i-1][k-1]）>=fabs（t））
{
t=a[i-1][k-1];
ik=i;
}
for（j=0;j{
t=a[k-1][j];
a[k-1][j]=a[ik-1][j];
a[ik-1][j]=t;
}
for（i=k+1;i<=N;i++）
{
m=a[i-1][k-1]/a[k-1][k-1];
for（j=k+1;j<=N;j++）
a[i-1][j-1]-=m*a[k-1][j-1];
}
}
x[N-1]=1.0;
for（k=N-1;k>0;k--）
{
t=0.0;
for（j=k+1;j<=N;j++）
t+=a[k-1][j-1]*x[j-1];
x[k-1]=-t/a[k-1][k-1];
}
}
voidQR（double（*a）[N],double（*q）[N],double（*r）[N]）
{
doubled,c,h,t,s,u[N],w[N],p[N];
inti,j,k;
for（i=0;ifor（j=0;j{
q[i][j]=i==j?
1:
0;
r[i][j]=a[i][j];
}
for（k=1;kif（fabs（a[k][k-1]）>E）
{
d=0.0;
for（i=k;i<=N;i++）
d+=a[i-1][k-1]*a[i-1][k-1];
d=sqrt（d）;
c=a[k-1][k-1]<=0?
d:
-d;
h=c*（c-a[k-1][k-1]）;
for（i=1;i<=N;i++）
u[i-1]=i0:
a[i-1][k-1];
u[k-1]=a[k-1][k-1]-c;
for（i=0;i{
t=0.0;
s=0.0;
for（j=0;j{
t+=q[i][j]*u[j];
s+=r[j][i]*u[j];
}
w[i]=t;
p[i]=s/h;
}
for（i=0;ifor（j=0;j{
q[i][j]-=w[i]*u[j]/h;
r[i][j]-=u[i]*p[j];
}
HAH（a,u,h,N）;
}
}
voidmain（）
{
doublea[N][N],lc[N][2],x[N],t,q[N][N],r[N][N],b[N][N];//lc存储复特征值
inti,j,k,m=0;
A（a）;//定义矩阵A
AN（a）;//对矩阵A进行拟上三角化
printf（"对矩阵A进行拟上三角化\n"）;
for（i=0;i{
for（j=0;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
SQR（a,lc）;//对拟上三角矩阵A进行带双步位移的QR分解，求解特征值
for（i=0;iif（fabs（lc[i][1]）m+=1;
printf（"所有%d个实特征值及其对应的特征向量\n",m）;
for（i=0,j=1;iif（fabs（lc[i][1]）{
printf（"λ%d=%.12e\n",j,lc[i][0]）;
printf（"X%d=",j）;
G（x,lc[i][0]）;//gauss消去法求特征向量
for（k=0;k{
printf（"\t%20.12e",x[k]）;
if（k==N/2-1||k==N-1）
printf（"\n"）;
}
j++;
}
printf（"所有%d对复特征值\n",（N-m）/2）;
for（i=0,j=m+1;iif（fabs（lc[i][1]）>E）
{
printf（"λ%d,%d=%.12e±%.12ei\n",j,j+1,lc[i][0],lc[i][1]）;
i++;
j+=2;
}
A（a）;//初始化矩阵A
AN（a）;//对矩阵A进行拟上三角化
for（k=0,t=1;k{
if（t>E）
{
for（i=0;ifor（j=0;jb[i][j]=a[i][j];
QR（a,q,r）;//对拟上三角阵A进行第k次QR分解，求Q、R及乘积矩阵RQ
}
else
{
printf（"\nQR分解次数为%d\n",k）;
break;
}
t=0.0;
for（i=0;iif（fabs（a[i+1][i]）t+=（b[i][i]-a[i][i]）*（b[i][i]-a[i][i]）;
else
{
t+=（b[i][i]+b[i+1][i+1]-a[i][i]-a[i+1][i+1]）*（b[i][i]+b[i+1][i+1]-a[i][i]-a[i+1][i+1]）;
i+=1;
}
t=sqrt（t）;
//以矩阵B-A的一阶对角块和二阶对角块的迹构成的向量的2-范数<ε判断QR分解法收敛
}
printf（"第%d次QR分解中，R=\n",k）;//输出第k次分解时的R
for（i=0;i{
for（j=0;jprintf（"%20.12e\t",fabs（r[i][j]）0:
r[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（r[i][j]）0:
r[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
printf（"第%d次QR分解中，Q=\n",k）;//输出第k次分解时的Q
for（i=0;i{
for（j=0;jprintf（"%20.12e\t",fabs（q[i][j]）0:
q[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（q[i][j]）0:
q[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
printf（"第%d次QR分解中，RQ=\n",k）;//输出第k次分解时的RQ
for（i=0;i{
for（j=0;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
for（i=0;i{
for（j=N/2;jprintf（"%20.12e\t",fabs（a[i][j]）0:
a[i][j]）;
printf（"\n"）;
}
printf（"\n"）;
}
三、计算结果
拟上三角矩阵
带双步位移的
分解法求特征值过程
所有特征值及实特征值对应的特征向量
对
采用
分解法求
、
及
，下面只给出最后一次迭代的结果
四、计算结果分析
1.对矩阵
进行拟上三角化可以减少
分解的迭代次数。
对矩阵
直接进行
分解，得到收敛的分块上三角阵所需次数为866次；对拟上三角阵
进行
分解，得到收敛的分块上三角阵所需次数为489次。
2.
分解法得到收敛的分块上三角阵，由一阶、二阶对角块求得的特征值与带双步位移的
分解法求得的特征值一致。
3.带双步位移的
分解法通过降阶大大减少了QR分解的次数。
采用带双步位移的
分解法求
的特征值，只需14次QR分解。
升级会员 