北航数值分析实习第二题Word文档格式.docx
- 文档编号:3596753
- 上传时间:2023-05-02
- 格式:DOCX
- 页数:17
- 大小:267.83KB
北航数值分析实习第二题Word文档格式.docx
《北航数值分析实习第二题Word文档格式.docx》由会员分享,可在线阅读,更多相关《北航数值分析实习第二题Word文档格式.docx(17页珍藏版)》请在冰点文库上搜索。
voidQR();
//对矩阵进行QR分解
voidQRshuangbu();
//对矩阵进行带双步位移的QR分解
voidgauss();
//列主元的高斯消元法求解特征向量
inti,j,s,p,k,ik,nR,nC;
doubleA[n][n],q[n][n],r[n][n],rq[n][n],I[n][n];
doubleP[n],W[n],u[n],Q[n];
doubledr,cr,hr,ar,tr;
doubles1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n];
voidmain()
{for(i=1;
i<
n;
i++)
{for(j=1;
j<
j++)
{if(i==j)
{I[i][j]=1;
q[i][j]=1;
}
else{I[i][j]=0;
q[i][j]=0;
}}}
caculateA();
hessenberg();
QR();
QRshuangbu();
gauss();
printf("
运行时间为%d\n"
clock());
voidcaculateA()//计算矩阵A的系数
{
for(i=1;
{for(j=1;
{if(j!
=i)
A[i][j]=sin(0.5*i+0.2*j);
else
A[i][j]=1.52*cos(i+1.2*j);
}
voidhessenberg()//将A拟上三角化
for(s=1;
s<
n-2;
s++)
{for(ar=0.0,i=s+2;
ar+=A[i][s]*A[i][s];
if(ar<
=err)
continue;
else{
ar+=A[s+1][s]*A[s+1][s];
dr=sqrt(ar);
if(A[s+1][s]>
0)
cr=-dr;
else
cr=dr;
hr=cr*cr-cr*A[s+1][s];
for(i=1;
=s;
u[i]=0.0;
u[s+1]=A[s+1][s]-cr;
for(i=s+2;
u[i]=A[i][s];
for(j=1;
{for(P[j]=0.0,i=1;
P[j]+=A[i][j]*u[i]/hr;
for(tr=0.0,i=1;
{tr+=P[i]*u[i]/hr;
for(i=1;
{for(Q[i]=0.0,j=1;
Q[i]+=A[i][j]*u[j]/hr;
{W[i]=Q[i]-tr*u[i];
{for(j=1;
A[i][j]-=(W[i]*u[j]+u[i]*P[j]);
}
for(i=1;
{for(j=1;
if(fabs(A[i][j])<
err)
A[i][j]=0.0;
printf("
数值分析计算实习第二次作业\n"
);
A1班\n"
\n"
矩阵A经过拟上三角化后得到的矩阵为:
for(i=1;
{
for(j=1;
j++)
printf("
%1.12e,"
A[i][j]);
printf("
voidQR()//对矩阵进行QR分解
doubleu[n],w[n],F[n];
for(p=1;
p<
p++)
r[s][p]=A[s][p];
n-1;
{
for(dr=0.0,i=s;
dr+=r[i][s]*r[i][s];
dr=sqrt(dr);
if(dr<
continue;
else
{if(A[s][s]>
0)cr=-dr;
elsecr=dr;
hr=cr*cr-cr*r[s][s];
s;
u[i]=0;
u[s]=r[s][s]-cr;
for(i=s+1;
u[i]=r[i][s];
{for(F[i]=0.0,j=1;
F[i]+=r[j][i]*u[j]/hr;
{for(j=1;
r[i][j]=r[i][j]-u[i]*F[j];
for(i=1;
{for(w[i]=0.0,j=1;
w[i]+=q[i][j]*u[j];
i++)
{for(j=1;
q[i][j]-=w[i]*u[j]/hr;
}
}
for(i=1;
{for(j=1;
{for(rq[i][j]=0.0,s=1;
rq[i][j]+=r[i][s]*q[s][j];
}}
A经过QR分解后得到的正交矩阵Q为:
for(i=1;
if(fabs(q[i][j])<
q[i][j]=0.0;
{for(j=1;
q[i][j]);
A经过QR分解后得到的上三角矩阵R为:
if(fabs(r[i][j])<
r[i][j]=0.0;
r[i][j]);
正交矩阵Q乘以上三角矩阵R得到的RQ为:
if(fabs(rq[i][j])<
rq[i][j]=0.0;
rq[i][j]);
voidQRshuangbu()//对矩阵进行带双步位移的QR分解
intK=1,m=10;
nR=0,nC=0;
loop3:
if(fabs(A[m][m-1])<
{nR++;
tzR[nR]=A[m][m];
m--;
gotoloop4;
elsegotoloop5;
loop4:
if(m==1)
{
nR++;
tzR[nR]=A[1][1];
gotoloop9;
elseif(m==2)
{
s1=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];
x=s1*s1-4*t;
if(x>
=0)
nR++;
tzR[nR]=(s1+sqrt(x))/2;
tzR[nR]=(s1-sqrt(x))/2;
nC++;
tzC[0][nC]=s1/2;
tzC[1][nC]=sqrt(-x)/2;
tzC[1][nC]=-sqrt(-x)/2;
elsegotoloop3;
loop5:
if(fabs(A[m-1][m-2])<
s1=A[m-1][m-1]+A[m][m];
m--;
gotoloop4;
elsegotoloop6;
loop6:
{if(K>
=L)
printf("
计算终止,未能得到全部特征值\n"
elsegotoloop7;
loop7:
s1=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];
for(i=1;
=m;
i++)//生成M
for(j=1;
for(sum=0.0,s=1;
sum+=A[i][s]*A[s][j];
M[i][j]=sum-s1*A[i][j]+t*I[i][j];
for(s=1;
for(ar=0.0,i=s+1;
ar+=M[i][s]*M[i][s];
if(ar==0)
continue;
else
{
ar+=M[s][s]*M[s][s];
dr=sqrt(ar);
if(M[s][s]>
0)
cr=-dr;
elsecr=dr;
hr=cr*cr-cr*M[s][s];
for(i=1;
u[i]=0.0;
u[s]=M[s][s]-cr;
for(i=s+1;
u[i]=M[i][s];
for(j=1;
j++)//v
for(v[j]=0.0,i=1;
v[j]+=M[i][j]*u[i]/hr;
i++)//Mr+1
for(j=1;
M[i][j]-=u[i]*v[j];
j++)//P
for(P[j]=0.0,i=1;
P[j]+=A[i][j]*u[i]/hr;
i++)//Q
for(Q[i]=0.0,j=1;
Q[i]+=A[i][j]*u[j]/hr;
for(tr=0.0,i=1;
i++)//tr
tr+=P[i]*u[i]/hr;
i++)//W
W[i]=Q[i]-tr*u[i];
i++)//Ar+1
A[i][j]-=(W[i]*u[j]+u[i]*P[j]);
}
gotoloop8;
loop8:
{k++;
gotoloop3;
loop9:
;
}
矩阵的全部特征值为:
=nR;
%1.12e\n"
tzR[i]);
=nC;
%1.12e"
tzC[0][i]);
if(tzC[1][i]>
printf("
+%1.12e*i\n"
tzC[1][i]);
elseprintf("
%1.12e*i\n"
voidgauss()//列主元的高斯消元法求解特征向量
doublech,m[n],x[n][n];
for(p=1;
caculateA();
for(i=1;
A[i][j]-=tzR[p]*I[i][j];
for(k=1;
k<
k++)
{
for(ik=k,i=k;
if(A[ik][k]<
A[i][k])
ik=i;
for(j=k;
{
ch=A[ik][j];
A[ik][j]=A[k][j];
A[k][j]=ch;
}
for(i=k+1;
m[i]=A[i][k]/A[k][k];
for(j=k+1;
A[i][j]-=m[i]*A[k][j];
}
x[p][n-1]=1.0;
for(k=n-2;
k>
0;
k--)
for(ar=0.0,j=k+1;
=n;
ar-=A[k][j]*x[p][j]/A[k][k];
x[p][k]=ar;
}
printf("
对应特征值%1.12e的一个特征向量α为(全部特征向量为kα,k取不为零的所有实数):
tzR[p]);
printf("
x[p][i]);
3运行结果
运行结果如下:
-8.945216982281e-001,-9.933136491826e-002,-1.099831758877e+000,-7.665038709077e-
001,1.707601141456e-001,-1.934882558889e+000,-8.390208705246e-002,9.132565113143
e-001,-6.407977009188e-001,1.946733678685e-001,
-2.347878362416e+000,2.372057921598e+000,1.827998552316e+000,3.266556884714e-001
2.082360583635e-001,2.088987009941e+000,1.847861910289e-001,-1.263015266080e+00
0,6.790694668499e-001,-4.672150886500e-001,
0.000000000000e+000,1.735954469946e+000,-1.165023367477e+000,-1.246744443518e+00
0,-6.298225489084e-001,-1.984820180992e+000,2.975750060800e-001,6.339300596595e-
001,-1.308518928772e-001,3.040301036095e-001,
0.000000000000e+000,0.000000000000e+000,-1.292937563924e+000,-1.126239225902e+00
0,1.190782911924e+000,-1.308772983895e+000,1.860151662666e-001,4.236733936881e-0
01,-1.019600826545e-001,1.943660914505e-001,
0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,1.577711153032e+000,
8.169358328160e-001,4.461531723828e-001,-4.365092541609e-002,-4.665979167188e-00
1,2.941231566184e-001,-1.034421113665e-001,
0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,
-7.728975134989e-001,-1.601028244046e+000,-2.912685474827e-001,-2.434337858321e-
001,6.736286084510e-001,2.624772904937e-001,
0.000000000000e+000,-7.296773946362e-001,-7.965456279818e-003,9.710739102007e-00
1,-1.298967368574e-001,2.780242081241e-002,
0.000000000000e+000,0.000000000000e+000,7.945539612976e-001,-4.525143454606e-001
5.048901527575e-001,-1.211210193512e-001,
0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,7.039911373514e-001,
1.267535523498e-001,-3.714696735513e-001,
-4.919586872214e-001,4.081509766399e-001,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 实习 第二