北航数值分析大作业一

1数值分析大作业一1、算法设计方案由于算法要求所有零元素都不存储，因此采用带状矩阵来存储矩阵A,主要子程序为LU分解子程序、求解方程组子程序、带偏移量幂法求特征值子程序、带偏移量反幂法求特征值子程序具体算法设计思路如下：①利用幂法求按模最大的特征值max②使用带原点平移的幂法求maxAI的按模最大的特征值max，则maxmax为另外一特征值③比较max与maxmax的大小，较大值即为501，较小者即为1④使用反幂法求A的按模最小的特征值s，利用LU分解求取每一次迭代后的ku⑤使用带偏移量的反幂法，偏移量分别为5011140kk(1,2,,39)k，即可求得分别与之相距最近的特征值1239,,,iii⑥A为非奇异的是对称矩阵，因此12()ncondA，其中1和n分别为按模最大的特征值和按模最小的特征值⑦本题中LU分解采用doolittle方法，分解后所得U矩阵即为对A进行初等行列式变换后所得矩阵，切变换过程中行列式值不变，因此5011detdet()(1501)niiiAUui2、源程序#includemath.h#includestdio.h#includestdlib.h#defineN501//矩阵为501*501的矩阵#defines2//上半带宽#definer2//下半带宽3个宏定义为方便LU分解及求解方程组过程/**********全局变量定义**********************/doubleA[5][501];2doubleu[501],y[501];doublelambda[41];//lambda[0]为λ1及最小特征值，lambda[40]为λ501及最大特征值doublelambda_s,lambda_m;//按模最小(大)的特征值；doubleLU[5][501];/********子函数声明**************************/voidinit_A(doubleA[][501]);//初始化A矩阵doublemodule_value_u(doublett[]);//求u[501]的模值voidinit_u(doublett[]);//初始化迭代初始向量u0doublepower_method(doubleoffset);//带原点偏移的幂法，返回值为特征值doubleinverse_power_method(doubleoffset);//带原点偏移的反幂法，返回值为特征值，子函数中打印出偏移量，求得的特征值，误差，迭代次数voidLU_decomposition(doublec[][501]);//参数c为矩阵LU[5][501]首地址，程序进行完后，保存分解后的L和U缩减后的矩阵voidsolve(doublec[][501],doubleb[],doublex[]);////第1个参数为LU分解完的LU矩阵，第2个参数b为已知的右端值，第3个参数x为求得的解向量存储位置intmax_2(inta,intb);intmax_3(inta,intb,intc);intmin_2(inta,intb);voidmain()//主程序{inti;doubleuk;//偏移量doubledet;//A的行列式的值init_A(A);//初始化矩阵Alambda_m=power_method(0);//求按模最大的特征值lambda[40]=power_method(lambda_m);//求相反方向另一个端点的特征值if(lambda_mlambda[40])//若大小反向，则交换两个元素中的值，得到λ1和λ501{lambda[0]=lambda[40];lambda[40]=lambda_m;}else{3lambda[0]=lambda_m;}lambda_s=inverse_power_method(0);det=1;for(i=0;iN;i++)det=det*LU[s][i];for(i=1;i40;i++){uk=lambda[0]+(lambda[40]-lambda[0])*i/40;lambda[i]=inverse_power_method(uk);}printf(-----------Theresultsareasfollows-------------\n);printf(λ[1]=%1.11e\nλ[501]=%1.11e\n,lambda[0],lambda[40]);printf(λs=%1.11e\n,lambda_s);printf(cond(A)=%1.11e\n,fabs(lambda_m/lambda_s));printf(detA=%1.11e\n,det);for(i=1;i40;i++){printf(λ[i%d]=%1.11e,i,lambda[i]);if(i%3==0)printf(\n);}printf(\n);}voidinit_A(doubleA[][501])//初始化A矩阵{inti,j;for(i=0;i5;i++){if(abs(i-2)==0){for(j=0;j501;j++)A[i][j]=(1.64-0.024*(j+1))*sin(0.2*j+0.2)-0.64*exp(0.1/(j+1));}elseif(abs(i-2)==1){4for(j=0;j501;j++)A[i][j]=0.16;}elsefor(j=0;j501;j++)A[i][j]=-0.064;}}voidinit_u(doublett[])//初始化迭代初始向量,自变量为初始向量的数组{inti;for(i=0;i501;i++){tt[i]=1;}}doublemodule_value_u(doublett[])//求u[501]的模值{inti=501;doublet=0;for(i=0;i501;i++){t=t+tt[i]*tt[i];}returnsqrt(t);}doublepower_method(doubleoffset)//带原点偏移的幂法，返回值为特征值{doubleb=0,b0=0,e;doublem;//求得的向量模值inti,j,k=0;//k表示迭代次数init_u(u);do{k++;5b0=b;m=module_value_u(u);for(i=0;i501;i++){y[i]=u[i]/m;}for(i=0;i501;i++){u[i]=0;for(j=max_2(0,i-s);jmin_2(i+s+1,N);j++){u[i]=A[i-j+2][j]*y[j]+u[i];if(i==j)u[i]=u[i]-offset*y[j];}}b=0;for(i=0;i501;i++){b=b+y[i]*u[i];}}while(fabs(e=b-b0)1e-12);b=b+offset;printf(λ=%fe=%ek=%d\n,b,e,k);//分别为特征值，迭代误差，迭代次数return(b);}voidLU_decomposition(doublec[][501])//参数c为矩阵LU[5][501]首地址，程序进行完后，保存分解后的L和U缩减后的矩阵{intj,k,t;for(k=0;kN;k++)//求U矩阵的第k行以及L矩阵的第k列{for(j=k;jmin_2(k+s,N-1)+1;j++)//{for(t=max_3(0,k-r,j-s);tk;t++)//求U矩阵的行中各个元素的循环{c[k-j+s][j]=c[k-j+s][j]-c[k-t+s][t]*c[t-j+s][j];6}//每次计算完U的元素后才能计算L的元素，所以上下两个循环不能合并在一起if(c[s][k]==0){printf(a[%d][%d]=0,Unabletoshowthesolutionofequations\n,k,k);exit(1);}if(jmin_2(k+s,N-1)+1)//因为矩阵L没有N行，因此这里要加一个判断，符合条件才能进行下面的循环语句{for(t=max_3(0,j-r+1,k-s);tk;t++)//求L矩阵的行中各个元素的循环{c[j-k+s+1][k]=c[j-k+s+1][k]-c[j-t+s+1][t]*c[t-k+s][k];}c[j-k+s+1][k]=c[j-k+s+1][k]/c[s][k];}}}}voidsolve(doublec[][501],doubleb[],doublex[])//第1个参数为分解完的LU矩阵，第2个参数b为已知的右端值，第3个参数x为求得的解向量存储位置{inti,t;for(i=0;iN;i++)//此循环求出向量y{x[i]=b[i];for(t=max_2(0,i-r);ti;t++){x[i]=x[i]-c[i-t+s][t]*x[t];}}for(i=N-1;i=0;i--){for(t=i+1;tmin_2(i+s,N-1)+1;t++)7{x[i]=x[i]-c[i-t+s][t]*x[t];}x[i]=x[i]/c[s][i];}}doubleinverse_power_method(doubleoffset)//带原点偏移的反幂法，返回值为特征值{inti,j,k=0;doubleb=0,b0=0,e;for(i=0;i5;i++){if(i!=2)for(j=0;j501;j++)LU[i][j]=A[i][j];elsefor(j=0;j501;j++)LU[i][j]=A[i][j]-offset;}LU_decomposition(LU);//LU分解，只分解一次即可init_u(u);do{k++;b0=b;for(i=0;i501;i++)y[i]=u[i]/module_value_u(u);solve(LU,y,u);b=0;for(i=0;i501;i++){b=b+y[i]*u[i];}b=1/b;b=b+offset;}while(fabs(e=b-b0)1e-12);printf(offset=%1.11eλ=%1.11ee=%1.11ek=%d\n,offset,b,e,k);//分别为特征值，迭代误差，迭代次数returnb;}intmax_2(inta,intb)8{if(ab)a=b;returna;}intmax_3(inta,intb,intc){if(ab)a=b;if(ac)a=c;returna;}intmin_2(inta,intb){if(ab)a=b;returna;}3、计算结果如下所示9在幂法和反幂法中偏移量、特征值、误差及计算迭代次数也在计算过程中打印出来4、迭代过程中的现象、原因及解决方法迭代的初始向量对计算结果可能产生极为重要的影响，如果选取初始向量不恰当，则可能得到错误的结果。以最大（或最小）特征值——也即按模最大特征值为例：分别取一下两组初始向量：1.）（501,2,1;1][iiu2.501)2,3,(i0,1)(i,1][iu10当初始向量为）（501,2,1;1][iiu时可以得到按模最大的特征向量为：当初始向量为501)2,3,(i0,1)(i,1][iu时可以得到按模最大的特征向量为：由以上可以看出，初始向量对于计算结果有着十分正要的影响。原因分析：初始向量可以表示成：nnxxxu22110；假设||||||n21，则采用幂法迭代k次后有：))()((12122111nknnkkkxxxu，只有当01，时，才有111xukk；即幂法能够正常使用的前提是01，当初始向量中有若干0时，极有可能出现01的情况，此时有两种可能：一种可能是由于舍入误差影响，迭代若干次后x1