上机实习题一一、题目:已知A 与b12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.0317432.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 A= -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784137 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.00001b={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392} 1.用Household 变换,把A 化为三对角阵B (并打印B )。
2.用超松弛法求解BX=b (取松弛因子ω=1.4,X (0)=0,迭代9次)。
3.用列主元素消去法求解BX=b 。
二、解题方法的理论依据:1 、用Householder 变换的理论依据﹝1﹞ 令A0=A,a(ij)1=a(ij),已知Ar_1即Ar_1=a(ij)r ﹝2﹞ Sr=sqrt(pow(a,2))﹝3﹞ a(r)=Sr*Sr+abs(a(r+1,r))*Sr ﹝4﹞ y(r)=A(r_1)*u®/a®﹝5﹞ Kr=(/2)*Ur 的转置*Yr/a® ﹝6﹞ Qr=Yr-Kr*Ur﹝7﹞ Ar=A(r-1)-(Qr*Ur 的转置+Ur*Qr 的转置) r=1,2,,……,n -2 2 、用超松弛法求解其基本思想:在高斯方法已求出x (m),x (m-1)的基础上,组合新的序列,从而加快收敛速度。
其算式:⎪⎪⎩⎪⎪⎨⎧=+⋅=+-+----=][][0][0][~][]][[][]1[0]][[]1][[]1[]][[]1][[][~i X i X i X i X w i X i i B i b i X i i B i i B i X i i B i i B i X 其中ω是超松弛因子,当ω>1时,可以加快收敛速度 3 、用消去法求解用追赶消去法求Bx=b 的方法:][]1[1i b i d =+ , ]][1[]2[1i i B i a +=+ ,][]1[]1][[][]][[]1[]1][[i b i X i i B i X i i B i X i i B =+⋅++⋅+-⋅-]][[]1[1i i B i b =+ , ]1][[]1[1+=+i i B i c , q1[0]=0 , u1[0]=0 ,8,,2,1]),[1][1][1(][1][1⋅⋅⋅=⋅+-=i i q i a i b i c i q9,,2,1]),1[1][1][1(])[1][1][1(}[1⋅⋅⋅=-⋅+⋅-=i i q i a i b i u i a i d i u x[9]=u1[9]1,,7,8],[1]1[][1][⋅⋅⋅=++⋅=i i u i x i q i x三、1.计算程序:#include "math.h" #include "stdio.h" #define ge 8 void main() {int sign(double x); double a[][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.713846,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}}; double k,h,s,w; int i,j,n,m,g; double u[9],x1[9],y[9],q[9],b1[9][10],x[9]; double b[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}; for(j=0;j<7;++j) /*Household 变换 */ { s=0.0; for(i=j+1;i<9;++i) s=s+a[i][j]*a[i][j]; s=sqrt(s); h=(a[j+1][j]>0)?(s*s+s*a[j+1][j]):(s*s-s*a[j+1][j]); for(g=0;g<9;++g) {if (g<=j)u[g]=0;else if (g==j+1)u[g]=a[j+1][j]+s*sign(a[j+1][j]);else u[g]=a[g][j];}for(m=0;m<9;++m){y[m]=0;for(n=0;n<9;++n)y[m]=y[m]+a[m][n]*u[n];y[m]=y[m]/h;}k=0;for(i=0;i<9;++i)k=k+u[i]*y[i];k=0.5*k/h;for(i=0;i<9;++i)q[i]=y[i]-k*u[i];for(n=0;n<9;++n)for(m=0;m<9;++m)a[m][n]=a[m][n]-(q[m]*u[n]+u[m]*q[n]);}printf("Household:\n");for(i=0;i<9;++i)for(j=0;j<9;++j){if (j%9==0)printf("\n");printf("%-9.5f",a[i][j]);}printf("\n");w=1.4; /*超松弛法*/ for(i=0;i<9;i++)x1[i]=0;for(i=0;i<9;i++)for(j=0;j<9;j++){if(i==j)b1[i][j]=0;else b1[i][j]=-a[i][j]/a[i][i];}for(i=0;i<9;i++)b1[i][9]=b[i]/a[i][i];for(n=0;n<9;n++)for(i=0;i<9;i++){s=0;for(j=0;j<9;j++)s=s+b1[i][j]*x1[j];s=s+b1[i][9];x1[i]=x1[i]*(1-w)+w*s;}for(i=0;i<9;i++){if (i==5)printf("\n");printf("x%d=%-10.6f",i,x1[i]);}printf("\n");u[0]=a[0][0]; /*以下是消去法*/y[0]=b[0];for(i=1;i<9;i++){q[i]=a[i][i-1]/u[i-1];u[i]=a[i][i]-q[i]*a[i-1][i];y[i]=b[i]-q[i]*y[i-1];}x[ge]=y[ge]/u[ge];for(i=ge-1;i>=0;i--)x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i];for(i=0;i<9;i++){if (i==5)printf("\n");printf(" x%d=%-10.6f",i,x[i]);}}int sign(double x){int z;z=(x>=(1e-6)?1:-1);return(z);}2.打印结果:Household:12.38412 -4.89308 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000-4.89308 25.39842 6.49410 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000.00000 6.49410 20.61150 8.24393 0.00000 0.00000 0.00000 0.00000 0.000000.00000 0.00000 8.24393 23.42284 -13.880070.00000 0.00000 0.00000 -0.000000.00000 0.00000 0.00000 -13.8800729.69828 4.53450 0.00000 0.00000 0.000000.00000 0.00000 0.00000 0.00000 4.53450 16.00612 4.88144 0.00000 0.000000.00000 0.00000 0.00000 0.00000 0.00000 4.88144 26.01332 -4.50363 -0.000000.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -4.50363 21.25406 4.504500.00000 0.00000 0.00000 -0.00000 0.00000 0.00000 -0.00000 4.50450 14.53412x0=1.073409 x1=2.272579 x2= -2.856601x3=2.292514 x4=2.112165 x5= -6.422586x6=1.357802 x7=0.634259 x8= -0.587042x0=1.075799 x1=2.275744 x2= -2.855515x3=2.293099 x4=2.112634 x5= -6.423838x6=1.357923 x7=0.634244 x8= -0.587266四、问题讨论:此程序具有很好的通用性。