《数值分析》计算实习作业《一》
北航
第一题 设有501?501的矩阵
?a1b??ba2?cb?A????????其中
cba3ccbbcc????a499bba500cb0.1?????? c??b?a501??ai?(1.64?0.024i)sin(0.2i)?0.64ei(i?1,2,?501);b?0.16,c??0.064.矩阵的特征值?i(i?1,2,?,501)满足?1??2????501,|?s|?试求 1.
1?i?501min|?i|
?1，?501和?s的值
?501??140最接近的特征值?i?(??1,2,?,39)
2. 的与数?k??1??3. 的（谱范数）条件数cond(A)2和行列式detA 要求
1. 算法的设计方案（A的所有零元素都不能存储） 2. 全部源程序（详细注释）。变量为double，精度??10-12，输出为e型12位有效数字
3. 特征值?1，?501，?s和?i?(??1,2,?,39)以及cond(A)2,det4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因
A的值
解答：
1. 算法设计
对于?s满足|?s|?求得。
对于?1，?501，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设|1?i?501min|?i|，所以?s是按模最小的特征值，直接运用反幂法可
?1|?|?501|，然后运用幂法，看能否求得一个特征值，如果
可以求得一个，证明A是收敛的，求得的结果是正确的，然后对A进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是?501，较小的特征值就是?1。如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A不收敛，则需要将A进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。 2. 源程序（见附录A）
由于A一律采用带状存储，因此关于A的函数，都是针对带状存储的矩阵的运算。 3. 输出结果：
?1?-1.07001136150e+001
?501?9.72463409878e+000 ?s?-5.55791079423e-003
cond(A)2?1.92520427390e+003
detA?2.77278614175e+118
的与数?k??1???501??140最接近的特征值?i?(??1,2,?,39)
-1.01829340331e+001 -9.58570742507e+000 -9.17267242393e+000
-8.65228400790e+000 -8.09348380868e+000 -7.65940540769e+000 -7.11968464869e+000 -6.61176433940e+000 -6.06610322660e+000 -5.58510105263e+000 -5.11408352981e+000 -4.57887217687e+000 -4.09647092626e+000 -3.55421121575e+000 -3.04109001813e+000 -2.53397031113e+000 -2.00323076956e+000 -1.50355761123e+000 -9.93558606008e-001 -4.87042673885e-001 2.23173624957e-002
5.32417474207e-001 1.05289896269e+000 1.58944588188e+000 2.06033046027e+000 2.55807559707e+000 3.08024050931e+000 3.61362086769e+000 4.09137851045e+000 4.60303537828e+000 5.13292428390e+000 5.59490634808e+000 6.08093385703e+000 6.68035409211e+000 7.29387744813e+000 7.71711171424e+000 8.22522001405e+000 8.64866606519e+000 9.25420034458e+000
4. 讨论初始向量对计算结果的影响
如果初始迭代向量?0选取使得其在特征值?1上的分量为零，但是由于计算中存在误差，迭代次数充分大后理论上使得?k在特征值?1上的分量不为零于是以?k为新的初始迭代向量仍可以得到结果。
当增加迭代初始向量中非零元素的个数时，计算出的结果会越接近准确值。实际上，当初始向量?0中的0元素较多时，可能ai=0的情况较为普遍，许多a都有可能等于0，此时计算出的结果便与最大特征值差距较大。这是因为实际中由于精度的要求和算法允许的计算量无法做到迭代次数充分大，于是在实际选取中初始迭代向量最好使得其在按模最大特征值的分量上不为零。
附录A
#include
double mifa(double **a,int r,int s,int n,double *y); /*幂法函数*/ double fanmifa(double **a,int r,int s,int n,double *y); /*反幂法函数*/ int LU(double **A,double *B,int r,int s,int n,double *x); /*LU分解法*/ int det(double **A,int r,int s,int n,double *pdet); /*求矩阵行列式的值函数*/ double **dmatrix(int m,int n); /*生成坐标从1开始的矩阵*/ void free_dmatrix(double **a); /*释放由dmatrix申请的空间*/
int main() { int r=2,s=2,n=501,m=r+s+1; int i,j,k; double **A,*B; double lambda_1,lambda_s,lambda_501,temp; double uk[40],lambda_uk[40],y_uk[40][502]; double condA2,detA; double *y_1,*y_s,*y_501,*p_temp;/*特征值对应的特征向量*/ y_1=(double *)malloc((n+1)*sizeof(double)); y_s=(double *)malloc((n+1)*sizeof(double)); y_501=(double *)malloc((n+1)*sizeof(double)); B=(double *)malloc((n+1)*sizeof(double)); /*****初始化矩阵A*****/ A=dmatrix(m,n); for(i=1;i<=r+s+1;i++) { if(i == 1 || i == 5) for(j=1;j<=n;j++)
A[i][j]=-0.064; else if(i == 2 || i == 4) for(j=1;j<=n;j++) A[i][j]=0.16; else for(j=1;j<=n;j++) A[i][j]=B[j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j); }
A[1][1]=A[1][2]=A[2][1]=A[4][501]=A[5][500]=A[5][501]=0; /*****求lambda_s,用反幂法求按模最小的特征值*****/ lambda_s=fanmifa(A,r,s,n,y_s);
printf(\,lambda_s);
/*****用幂法求按模最大的特征值，所求为lambda_1或lambda_501*****/ lambda_1=mifa(A,r,s,n,y_1); /*带原点平移的~幂法~*/ for(j=1;j<=n;j++) A[s+1][j]=B[j]-lambda_1; lambda_501=mifa(A,r,s,n,y_501); lambda_501+=lambda_1;
if(lambda_501 - lambda_1 < 1e-12) /***注意浮点数比较大小的不可靠性***/ { temp=lambda_1;/*交换*/ lambda_1=lambda_501; lambda_501=temp; p_temp=y_1; y_1=y_501; y_501=p_temp; }
printf(\,lambda_1,lambda_501); /*****求和Uk最接近的特征值*****/ for(k=1;k<=39;k++) { uk[k]=lambda_1+k*(lambda_501-lambda_1)/40; for(j=1;j<=n;j++) A[s+1][j]=B[j]-uk[k]; lambda_uk[k]=fanmifa(A,r,s,n,y_uk[k]); lambda_uk[k]+=uk[k]; /**加上原点的偏移量**/ }
for(i=1;i<=39;i++) { printf(\,i,lambda_uk[i]); }
printf(\);
/*****求A的(谱范数)条件数cond(A)2*****/ condA2=fabs(lambda_1/lambda_s); if(condA2 < 1) condA2=1/condA2;
printf(\,condA2); /*****求A的行列式的值*****/ for(i=1;i<=n;i++)/*恢复矩阵A*/ A[s+1][i]=B[i];
if(det(A,r,s,n,&detA)) printf(\,detA); free_dmatrix(A); /*扫尾工作*/ free(B); free(y_1); free(y_501); free(y_s); return 0; }
/*------------------------------------------------------ ***函数名：mifa(double **a,int r,int s,int n,double *y) ***功 能：计算矩阵A按模最大的特征值和对应的特征向量 ***参 数：矩阵A,y为所求特征值对应的特征向量
***返回值：特征值(double) ***说 明：~~使用的是2范数~~ by lidezheng 日期：2013-10-19
------------------------------------------------------*/ double mifa(double **a,int r,int s,int n,double *y) { double beta1,beta2,max; int k,i,j,up,down; double *u; u=(double *)malloc((n+1)*sizeof(double)); /*选取非零向量u*/ for(i=1;i<=n;i++) { /***--->###这里要讨论初始向量对结果的影响###<---***/ *(u+i)=1; } /*幂法迭代*/ k=0; beta2=0; do{ k++; beta1=beta2; max=0; for(i=1;i<=n;i++)/*求出向量的2范数*/ { max+=pow(u[i],2); } max=sqrt(max); /*求y(k-1)*/ for(i=1;i<=n;i++) { y[i]=u[i]/max; } /*求Uk*/ for(i=1;i<=n;i++) { u[i]=0; down=((i-r))>1?(i-r):1; up=((i+s)>n)?n:(i+s); for(j=down;j<=up;j++) u[i]+=a[i-j+s+1][j]*y[j];/*****此处下标最重要******/ }