搜索
《数值分析B》大作业一一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501].由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素aij的方法是：A的带内元素aij=C中的元素ci-j+2,j2.求解λ1，λ501，λs①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs；如果λmax0,则λ501=λmax；如果λmax0,则λ1=λmax。②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max，如果λmax0,则λ1=λ，max+p；如果λmax0,则λ501=λ，max+p。3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik（k=1,2，…，39）。使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。4.求解A的（谱范数）条件数cond（A）2和行列式detA。①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。二．源程序#includestdio.h#includeiostream.h#includestdlib.h#includemath.h#includefloat.h#includeiomanip.h#includetime.h#defineE1.0e-12/*定义全局变量相对误差限*/intmax2(inta,intb)/*求两个整型数最大值的子程序*/{if(ab)returna;elsereturnb;}intmin2(inta,intb)/*求两个整型数最小值的子程序*/{if(ab)returnb;elsereturna;}intmax3(inta,intb,intc)/*求三整型数最大值的子程序*/{intt;if(ab)t=a;elset=b;if(tc)t=c;return(t);}voidassignment(doublearray[5][501])/*将矩阵A转存为数组C[5][501]*/{inti,j,k;//所有元素归零for(i=0;i=4;){for(j=0;j=500;){array[i][j]=0;j++;}i++;}//第0,4行赋值for(j=2;j=500;){k=500-j;array[0][j]=-0.064;array[4][k]=-0.064;j++;}//第1,3行赋值for(j=1;j=500;){k=500-j;array[1][j]=0.16;array[3][k]=0.16;j++;}//第2行赋值for(j=0;j=500;){k=j;j++;array[2][k]=(1.64-0.024*j)*sin((double)(0.2*j))-0.64*exp((double)(0.1/j));}}doublemifa(doubleu[501],doublearray[5][501],doublep)/*带原点平移的幂法*/{inti,j;/*u[501]为初始迭代向量*/doublea,b,c=0;/*array[5][501]为矩阵A的转存矩阵*/doubley[501];/*p为平移量*/for(;;){a=0;b=0;/*选用第一种迭代格式*///求ηk-1for(i=0;i=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求yk-1for(i=0;i=500;i++){y[i]=u[i]/a;}//求ukfor(i=0;i=500;i++){u[i]=0;for(j=max2(i-2,0);j=min2(i+2,500);j++){u[i]+=array[i-j+2][j]*y[j];}u[i]=u[i]-p*y[i];/*引入平移量*/}//求βkfor(i=0;i=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)=E)/*达到精度水平，迭代终止*/break;c=b;}return(b+p);/*直接返回A的特征值*/}voidchuzhi(doublea[])/*用随机数为初始迭代向量赋值*/{inti;srand((int)time(0));for(i=0;i=500;i++){a[i]=(10.0*rand()/RAND_MAX);/*生成0~10的随机数*/}}voidchuzhi2(doublea[],intj)/*令初始迭代向量为ei*/{inti;for(i=0;i=500;i++){a[i]=0;}a[j]=1;}voidLU(doublearray[5][501])/*对矩阵A进行Doolittle分解*/{/*矩阵A转存在C[5][501]中*/intj,k,t;/*分解结果L，U分别存在C[5][501]的上半部与下半部*/for(k=0;k=500;k++){for(j=k;j=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t=(k-1);t++){array[k-j+2][j]-=array[k-t+2][t]*array[t-j+2][j];}}if(k500)for(j=k+1;j=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t=(k-1);t++){array[j-k+2][k]-=array[j-t+2][t]*array[t-k+2][k];}array[j-k+2][k]=array[j-k+2][k]/array[2][k];}}}doublefmifa(doubleu[501],doublearray[5][501],doublep){/*带原点平移的反幂法*/inti,j;doublea,b,c=0;doubley[501];//引入平移量for(i=0;i=500;i++){array[2][i]-=p;}//先将矩阵Doolittle分解LU(array);for(;;){a=0;b=0;//求ηk-1for(i=0;i=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求yk-1for(i=0;i=500;i++){y[i]=u[i]/a;}//回带过程，求解ukfor(i=0;i=500;i++){u[i]=y[i];}for(i=1;i=500;i++){for(j=max2(0,(i-2));j=(i-1);j++){u[i]-=array[i-j+2][j]*u[j];}}u[500]=u[500]/array[2][500];for(i=499;i=0;i--){for(j=i+1;j=min2((i+2),500);j++){u[i]-=array[i-j+2][j]*u[j];}u[i]=u[i]/array[2][i];}//求βkfor(i=0;i=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)=E)/*达到精度要求，迭代终止*/break;c=b;}return(p+(1/b));/*直接返回距离原点P最接近的A的特征值*/}//主函数main(){inti;doubled1,d501,ds,d,a;doubleu[501];doubleMatrixC[5][501];printf(《数值分析》计算实习题目第一题\n);printf(SY1103120朱舜杰\n);//将矩阵A转存为MatrixCassignment(MatrixC);//用带原点平移的幂法求解λ1，λ501chuzhi(u);d=mifa(u,MatrixC,0);chuzhi(u);a=mifa(u,MatrixC,d);if(d0){d1=d;d501=a;}else{d501=d;d1=a;}printf(λ1=%.12e\n,d1);printf(λ501=%.12e\n,d501);//用反幂法求λschuzhi(u);ds=fmifa(u,MatrixC,0);printf(λs=%.12e\n,ds);//用带原点平移的反幂法求λikfor(i=1;i=39;i++){a=d1+(i*(d501-d1))/40;assignment(MatrixC);chuzhi(u);d=fmifa(u,MatrixC,a);printf(与μ%02d=%+.12e最接近的特征值λi%02d=%+.12e\n,i,a,i,d);}//求A的条件数d=fabs((d1/ds));printf(A的（谱范数）条件数condA2=%.12e\n,d);//求detAassignment(MatrixC);LU(MatrixC);a=1;for(i=0;i=500;i++){a*=MatrixC[2][i];}printf(行列式detA=%.12e\n,a);//测试不同迭代初始向量对λ1计算结果的影响。printf(改变迭代初始向量对λmax计算结果的测试如下：\n);assignment(MatrixC);for(i=0;i=500;i++){chuzhi2(u,i);d1=mifa(u,MatrixC,0);printf(u%03d,λmax=%+e,i,d1);if(((i+1)%3)==0)printf(\n);}printf(Pressanykeytocontinue\n);getchar();}三.程序结果：四．分析初始向量选择对计算结果的影响矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到的特征值。以幂法为例（反幂法原理相同），选取初始迭代向量ui=ei(i=0,1，…500)；即u[j]=0,j≠i；u[j]=1,j=i。测试结果如下：试验结果发现只有当i取特定的一些值时才能得到正确的结果，即得到按摸最大的特征值。i取不同值时，即取不同的初始向量时，可能得到不同的特征值。这是因为以A的n个线性无关的特征向量为一组基，将初始向量线性表出时，按摸最大的特征值λ1对应的特征向量x1的系数α1如果为0，就无法求出特征值λ1。如果按摸第二大的特征值λ2对应的特征向量x2的系数α2不为0，则求出该特征值。若为0，则以此类推。