《数值分析》上机实习Word文件下载.docx

《数值分析》上机实习Word文件下载.docx

《《数值分析》上机实习Word文件下载.docx》由会员分享，可在线阅读，更多相关《《数值分析》上机实习Word文件下载.docx（16页珍藏版）》请在冰豆网上搜索。

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<

stdio.h>

math.h>

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;

=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++）

{ur[i][0]=a[i][r];

for（j=0;

j<

j++）/*求向量yr*/

h=0;

for（i=0;

i++）/*求矩阵Ar-1与ur的乘积*/

{h=a[j][i]*ur[i][0]+h;

h=h/xr;

yr[j][0]=h;

i++）/*求kr*/

{urt[0][i]=ur[i][0];

kr=0;

{kr=urt[0][i]*yr[i][0]+kr;

kr=kr/（2*xr）;

i++）/*求qr*/

{qr[i][0]=yr[i][0]-kr*ur[i][0];

{qrt[0][i]=qr[i][0];

j++）/*求qr乘ur的转置*/

{m[j][i]=qr[j][0]*urt[0][i];

j++）/*求ur乘qr的转置*/

{t[j][i]=ur[j][0]*qrt[0][i];

j++）/*把上面求得的两个矩阵相加*/

{t[j][i]=m[j][i]+t[j][i];

j++）/*求Ar*/

{a[j][i]=a[j][i]-t[j][i];

printf（"

housholdis:

\n"

）;

j++）

{printf（"

%.6f"

a[j][i]）;

printf（"

\n"

/*超松驰法求方解*/

for（i=1;

10;

i++）/*迭代九次*/

{

for（j=1;

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;

8;

j++）/*求X1到X7的一次迭代*/

for（k=0;

k<

j;

k++）c1=c1+a[j][k]*x[k][i]/a[j][j];

/*公式中的第一个求和式*/

for（k=j+1;

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;

}

for（j=0;

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];

sorresultare:

i++）printf（"

%.5f,"

x[i][8]）;

/*列主元素求解*/

i++）/*比较大小，如果同列中有大于该列中对角元素值的就把该两行交换*/

if（a[i+1][i]>

a[i][i]）

{for（j=i,k=0;

i+3;

j++,k++）{z[k]=a[i][j];

a[i][j]=a[i+1][j];

a[i+1][j]=z[k];

c3=g[i];

g[i]=g[i+1];

g[i+1]=c3;

c5=a[i+1][i];

for（j=i;

{a[i+1][j]=a[i+1][j]-a[i][j]*c5/a[i][i];

g[i+1]=g