《数值分析》上机实习.docx
《《数值分析》上机实习.docx》由会员分享，可在线阅读，更多相关《《数值分析》上机实习.docx（16页珍藏版）》请在冰豆网上搜索。
《数值分析》上机实习
2006级研究生
《数值分析》上机试题
■院系电气工程学院
■专业测检技术与自动化装置
■姓名姜启
■学号2006020067
第一题
▲题目：
已知
12.3841202.115237-1.0610741.112336-0.1135840.7187191.7423823.067813-2.031743
2.11523719.141823-3.125432-1.0123452.1897361.563849-0.7841651.1123483.123124
-1.061074-3.12543215.5679143.1238482.0314541.836742-1.0567810.336993-1.010103
A=1.112336-1.0123453.12384827.1084374.101011-3.7418562.101023-0.718280-0.037585
-0.1135842.1897362.0314544.10101119.8979190.431637-3.1112232.1213141.784317
0.7187191.5638491.836742-3.7418560.4316379.789365-0.103458-1.1034560.238417
1.742382-0.784165-1.0567812.101023-3.111223-0.10345814.7138463.123789-2.213474
3.0678131.1123480.336993-0.7182802.121314-1.1034563.12378930.7193344.446782
-2.0317433.123124-1.010103-0.0375851.7843170.238417-2.2134744.44678240.00001
b=（2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392）;
（1）用Hoseholder变换，把A化为三对角矩阵B（并打印B）
（2）用超松驰法求解BX=b（取松驰因子W=1.4，X0=0，迭代九次）
（3）用列主元素消去法求解BX=b
▲解题方法的理论依据：
①Household变换：
任何实非奇异矩阵A，都可通过初等反射阵相似变化化为HOSSENBERG阵，因此若A是对称阵，则可化为相似的三对角阵，从而解决了实对称三对角化问题。
本题中，可依次对矩阵的列向量应用Household变换，完成矩阵的三对角化。
算法：
（1）令A0=A,aij
（1）=aij,已知Ar-1 即Ar-1=（aij（r））；
（2）Sr=（∑ni=r+1（air（r））2）1/2;
（3）ar=Sr2+|a（r）r+1,r|Sr;
ur=[0,…0,a（r）r+1,r+Sign（a（r）r+1,r）Sr,a（r）r+2,r,…,a（r）nr]T
（4）yr=Ar-1ur/ar;
（5）kr=0.5*urTyr/ar;
（6）qr=yr-krur;Ar=Ar-1-（qrurT+urqrT）;r=1,2….n-2.
②松驰法：
在GS方法法已求出x（m）基础上，经过重新组合而得的新序列，此新序列使收敛加快
算法如下：
③列列主消元法：
为改善Guass消去法的缺陷
▲主程序如下：
#include
#include
main（）
{intn=9,r,i,j,k;
floatsr=0,h=0,xr,kr=0;
floatw=1.4,c=0,c1=0,c2=0,c3=0,c4=0,c5=0;
floatz[3],d[9];
floatur[9][1];
floaturt[1][9];
floatyr[9][1];
floatqr[9][1];
floatqrt[1][9];
floatm[9][9];
floatt[9][9];
floatx[9][10]={0};
floatg[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
floata[9][9]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};
/*HOUSEHLODER变换*/
for（r=0;r<7;r++）/*总循环次数，r=n-2*/
{
sr=0;
for（i=r+1;i<9;i++）/*求SR*/
{sr=a[i][r]*a[i][r]+sr;}
sr=sqrt（sr）;
xr=sr*sr+fabs（a[r+1][r]）*sr;/*求公式中的ar*/
for（i=0;i<=r;i++）/*求向量ur，把ur设成一列九行的数组*/
{ur[i][0]=0;}
ur[r+1][0]=a[r+1][r]+sr*a[r+1][r]/fabs（a[r+1][r]）;/*符号函数用A/|A|的形式表式*/
for（i=r+2;i<9;i++）
{ur[i][0]=a[i][r];}
for（j=0;j<9;j++）/*求向量yr*/
{
h=0;
for（i=0;i<9;i++）/*求矩阵Ar-1与ur的乘积*/
{h=a[j][i]*ur[i][0]+h;}
h=h/xr;
yr[j][0]=h;
}
for（i=0;i<9;i++）/*求kr*/
{urt[0][i]=ur[i][0];}
kr=0;
for（i=0;i<9;i++）
{kr=urt[0][i]*yr[i][0]+kr;}
kr=kr/（2*xr）;
for（i=0;i<9;i++）/*求qr*/
{qr[i][0]=yr[i][0]-kr*ur[i][0];}
for（i=0;i<9;i++）
{qrt[0][i]=qr[i][0];}
for（j=0;j<9;j++）/*求qr乘ur的转置*/
{
for（i=0;i<9;i++）
{m[j][i]=qr[j][0]*urt[0][i];}
}
for（j=0;j<9;j++）/*求ur乘qr的转置*/
{
for（i=0;i<9;i++）
{t[j][i]=ur[j][0]*qrt[0][i];}
}
for（j=0;j<9;j++）/*把上面求得的两个矩阵相加*/
{
for（i=0;i<9;i++）
{t[j][i]=m[j][i]+t[j][i];}
}
for（j=0;j<9;j++）/*求Ar*/
{
for（i=0;i<9;i++）
{a[j][i]=a[j][i]-t[j][i];}
}
}
printf（"housholdis:
\n"）;
for（j=0;j<9;j++）
{
for（i=0;i<9;i++）
{printf（"%.6f",a[j][i]）;}
printf（"\n"）;
}
/*超松驰法求方解*/
for（i=1;i<10;i++）/*迭代九次*/
{
for（j=1;j<9;j++）c2=c2+a[0][i]*x[j][i-1]/a[0][0];/*求X0的一次迭代*/
c=w*（g[0]/a[0][0]-c2）;
x[0][i]=c-0.4*x[0][i-1];
c2=0;
for（j=1;j<8;j++）/*求X1到X7的一次迭代*/
{
for（k=0;kfor（k=j+1;k<9;k++）c2=c2+a[j][k]*x[k][i-1]/a[j][j];/*公式中的第二个求和式*/
c=w*（g[j]/a[j][j]-c1-c2）;
x[j][i]=c-0.4*x[j][i-1];
c1=0;
c2=0;
}
for（j=0;j<8;j++）c1=c1+a[8][j]*x[j][i]/a[8][8];/*求X8的一次迭代*/
c=w*（g[8]/a[8][8]-c1）;
x[8][i]=c-0.4*x[8][i-1];
c1=0;
}
printf（"sorresultare:
\n"）;
for（i=0;i<9;i++）printf（"%.5f,",x[i][8]）;
printf（"\n"）;
/*列主元素求解*/
for（i=0;i<8;i++）/*比较大小，如果同列中有大于该列中对角元素值的就把该两行交换*/
{
if（a[i+1][i]>a[i][i]）
{for（j=i,k=0;jc3=g[i];g[i]=g[i+1];g[i+1]=c3;}
c5=a[i+1][i];
for（j=i;j{a[i+1][j]=a[i+1][j]-a[i][j]*c5/a[i][i];}
g[i+1]=g
升级会员 