北航数值分析大作业第二次的Word下载.docx
《北航数值分析大作业第二次的Word下载.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业第二次的Word下载.docx（22页珍藏版）》请在冰点文库上搜索。
源程序：
#include<
stdio.h>
math.h>
conio.h>
#definee0.000000000001/*设置精度水平*/
#defineN10/*设置矩阵阶数*/
#defineL2500/*设置迭代次数*/
/*子程序*/
/*sgn函数*//*保证ci与aii符号相反*/
doublesgn（doublea）
{
doublex;
if（a>
e）x=1;
elsex=-1;
returnx;
}
/*求根函数*//*计算双步位移法中的两个位移量，即右下角二阶子式的两个特征值*/
doublepig（doubleb,doublec）
doublem;
m=b*b-4*c;
returnm;
/*矩阵的拟上三角化*//*只要对A1进行拟上三角化就能保证每个矩阵Ak都为拟上三角矩阵，大大简化计算量*/
voidhessenberg（doubleA[N][N]）
inti,j,k;
intm=0;
doubled,c,h,t;
doubleu[N],p[N],q[N],w[N];
for（i=0;
i<
N-2;
i++）
{
for（j=i+2;
j<
N;
j++）if（A[j][i]<
=e）m=m+1;
if（m==（N-2-i））continue;
for（j=i+1,d=0;
j++）d=d+A[j][i]*A[j][i];
d=sqrt（d）;
c=-1*sgn（A[i+1][i]）*d;
h=c*c-c*A[i+1][i];
for（j=i+2;
j++）u[j]=A[j][i];
for（j=0;
i+2;
j++）u[j]=0;
u[i+1]=A[i+1][i]-c;
j++）
{for（k=i+1,p[j]=0;
k<
k++）p[j]=A[k][j]*u[k]+p[j];
p[j]=p[j]/h;
}
{for（k=i+1,q[j]=0;
k++）q[j]=A[j][k]*u[k]+q[j];
q[j]=q[j]/h;
for（j=0,t=0;
j++）t=t+p[j]*u[j];
t=t/h;
j++）w[j]=q[j]-t*u[j];
j++）
{
for（k=0;
k++）A[j][k]=A[j][k]-w[j]*u[k]-u[j]*p[k];
}
}
/*矩阵的QR分解*/
voidQRdiv（doubleA[N][N],doubleQ[N][N],doubleR[N][N]）
//intm=0;
doubled,c,h;
doubleu[N],w[N],p[N];
{for（j=0;
j++）{if（i==j）Q[i][j]=1;
elseQ[i][j]=0;
}}
j++）R[i][j]=A[i][j];
N-1;
//for（j=i+1;
j++）if（R[j][i]<
=E）m=m+1;
//if（m==（N-1-i））continue;
for（j=i,d=0;
j++）d=d+R[j][i]*R[j][i];
c=-1*sgn（R[i][i]）*d;
h=c*c-c*R[i][i];
for（j=i+1;
j++）u[j]=R[j][i];
i;
u[i]=R[i][i]-c;
j++）{for（k=0,w[j]=0;
k++）w[j]=Q[j][k]*u[k]+w[j];
j++）{for（k=0;
k++）Q[j][k]=Q[j][k]-w[j]*u[k]/h;
{for（k=i,p[j]=0;
k++）p[j]=R[k][j]*u[k]+p[j];
k++）R[j][k]=R[j][k]-u[j]*p[k];
/*求解矩阵的所有特征值*/
voidcharacteristic（doubleA[N][N],doublechaR[N],doublechaI[N]）
intk=0,m=N-1;
inti,j;
intl;
doubles,t,x;
doubleM[N][N],B[N][N];
intf=0;
doubleQ[N][N],R[N][N];
for（l=0;
l<
L;
l++）
next:
if（m==0）{chaR[0]=A[0][0];
chaI[0]=0;
break;
if（fabs（A[m][m-1]）<
=e）
chaR[m]=A[m][m];
chaI[m]=0;
m--;
gotonext;
}
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];
if（m==1）
x=pig（s,t）;
if（x>
=e）{x=sqrt（x）;
chaR[m]=s/2+x/2;
chaR[m-1]=s/2-x/2;
chaI[m-1]=0;
}
else{x=sqrt（fabs（x））;
chaR[m]=s/2;
chaR[m-1]=s/2;
chaI[m]=x/2;
chaI[m-1]=-x/2;
break;
if（fabs（A[m-1][m-2]）<
m=m-2;
gotonext;
for（i=0;
=m;
for（j=0;
{
if（i==j）
{
for（k=0,M[i][j]=0;
k++）M[i][j]=A[i][k]*A[k][j]+M[i][j];
M[i][j]=M[i][j]-s*A[i][j]+t;
}
else
M[i][j]=M[i][j]-s*A[i][j];
}
/*带双步位移的QR分解中M的分解*/
for（i=0;
{for（j=0;
j++）R[i][j]=M[i][j];
m;
for（j=i+1;
=e）f=f+1;
if（f==（m-i））continue;
for（j=i,d=0;
d=sqrt（d）;
c=-1*sgn（R[i][i]）*d;
h=c*c-c*R[i][i];
for（j=0;
u[i]=R[i][i]-c;
{for（k=0,w[j]=0;
{for（k=0;
for（k=0;
}
for（k=0;
k++）M[j][k]=Q[j][k];
}
for（i=0;
{for（k=0,B[i][j]=0;
k++）B[i][j]=M[k][i]*A[k][j]+B[i][j];
{for（k=0,A[i][j]=0;
k++）A[i][j]=B[i][k]*M[k][j]+A[i][j];
/*求矩阵的所有特征值的特征向量*/
voideigenvector（doubleV[N][N],doubleT[N]）
doubleA[N][N],baoz[N][N],guod[N];
doublec;
inti,j,k,m,t;
intW=0;
i++）for（j=0;
j++）baoz[i][j]=V[i][j];
for（t=0;
t<
6;
t++）
j++）A[i][j]=baoz[i][j];
i++）A[i][i]=A[i][i]-T[t];
i++）
for（j=i;
j++）if（fabs（A[j][i]）>
e）{k=j;
for（j=i;
j++）{guod[j]=A[i][j];
A[i][j]=A[k][j];
A[k][j]=guod[j];
c=A[j][i];
if（fabs（c）>
e）for（m=i;
m<
m++）A[j][m]=A[j][m]/c;
}
c=A[j][i];
if（j!
=i）{for（m=i;
m++）A[j][m]=A[j][m]-A[i][m]*c;
V[t][N-1]=-1;
for（i=N-2;
i>
=0;
i--）
V[t][i]=A[i][N-1];
/*主函数*/
voidmain（）
doublea[N][N],b[N][N],chaR[N],chaI[N];
doubleq[N][N],r[N][N],qr[N][N];
doubleDOG[N];
doublef,g;
inti,j,k;
if（i!
=j）a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
elsea[i][j]=1.5*cos（（i+1）+1.2*（j+1））;
hessenberg（a）;
printf（"
拟上三角化后A（n-1）:
\n"
）;
N-5;
if（fabs（a[i][j]）<
e）a[i][j]=0;
printf（"
%22.11e"
a[i][j]）;
}
printf（"
for（j=N-5;
QRdiv（a,q,r）;
生成的上三角矩阵R：
j++）{if（fabs（r[i][j]）<
e）r[i][j]=0;
printf（"
%22.12e"
r[i][j]）;
}printf（"
生成的正交矩阵Q：
j++）{if（fabs（q[i][j]）<
e）q[i][j]=0;
q[i][j]）;
j++）for（k=0,qr[i][j]=0;
k++）qr[i][j]=qr[i][j]+r[i][k]*q[k][j];
生成的R*Q阵：
j++）{if（fabs（qr[i][j]）<
e）qr[i][j]=0;
qr[i][j]）;
characteristic（a,chaR,chaI）;
for（i=1;
if（i<
3）{f=chaR[i];
g=chaI[i];
chaR[i]=chaR[7+i];
chaI[i]=chaI[7+i];
chaR[7+i]=f;
chaI[7+i]=g;
if（i==5）{f=chaR[i];
chaR[i]=chaR[7];
chaI[i]=chaI[7];
chaR[7]=f;
chaI[7]=g;
矩阵的全部特征值：
for（j=0;
if（fabs（chaI[j]）<
=e）printf（"
λ%2d=%19.11e\n"
j+1,chaR[j]）;
elseif（chaI[j]>
e）printf（"
λ%2d=%18.11e+%18.11ei\n"
j+1,chaR[j],chaI[j]）;
elseprintf（"
λ%2d=%19.11e%19.11ei\n"
}//重新输入矩阵A
λ%d对应的特征向量为：
\n"
（i+1））;
j++）printf（"
for（j=N-5;
getch（）;
计算结果：
拟上三角化后A（n-1）
-0.894521698228-0.0993313649183-1.09983175888-0.7665038709080.170760
114146-1.93488255889-0.08390208705250.913256511314-0.640797700919
0.194673367868
-2.347878362422.37205792161.827998552320.3266556884710.208236058364
2.088987009940.184********9-1.263015266080.67906946685-0.46721508865
01.73595446995-1.16502336748-1.24674444352-0.629822548908-1.98482
0180990.297575006080.633930059659-0.130********70.30403010361
0-3.42119688154e-017-1.29293756392-1.12623922591.19078291192
-1.30877298390.186********70.423673393688-0.1019600826550.194366091451
05.94010094099e-017-3.28119028842e-0171.577711153030.816935
8328160.446153172383-0.0436509254161-0.4665979167190.294123156618
-0.103442111366
0-3.49494913888e-0178.59590318918e-0179.17131065135e-018
-0.772897513499-1.60102824405-0.291268547483-0.2434337858320.673628608451
0.262477290494
0-2.13633528681e-016-3.94299506639e-017-8.80991580387e-017
0-0.729677394636-0.007965456279830.971073910201-0.129896736857
.027*********
0-6.48100629069e-0178.50487510592e-017-5.02328897109e-018
0-5.3832322921e-0170.794553961298-0.4525143454610.504890152758
-0.121121019351
07.034391
升级会员 