《北航数值分析-实习作业1(C语言详细注释)》由会员分享，可在线阅读，更多相关《北航数值分析-实习作业1(C语言详细注释)（10页珍藏版）》请在金锄头文库上搜索。
1、数值分析计算实习作业一北航第一题 设有的矩阵其中矩阵的特征值满足试求 1. 和的值2. 的与数最接近的特征值3. 的（谱范数）条件数和行列式要求1. 算法的设计方案（A的所有零元素都不能存储）2. 全部源程序（详细注释）。变量为double，精度，输出为e型12位有效数字3. 特征值和以及的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于满足，所以是按模最小的特征值，直接运用反幂法可求得。对于，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A是收敛的，求得的结果
2、是正确的，然后对A进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是，较小的特征值就是。如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A不收敛，则需要将A进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。2. 源程序（见附录A）由于A一律采用带状存储，因此关于A的函数，都是针对带状存储的矩阵的运算。3. 输出结果：的与数最接近的特征值-1.01829340331e+001 -9.58570742507e+000 -9.17267242393e+000-8.65228400790e+000 -8.093
3、48380868e+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.87042673
4、885e-001 2.23173624957e-0025.32417474207e-001 1.05289896269e+000 1.58944588188e+0002.06033046027e+000 2.55807559707e+000 3.08024050931e+0003.61362086769e+000 4.09137851045e+000 4.60303537828e+0005.13292428390e+000 5.59490634808e+000 6.08093385703e+0006.68035409211e+000 7.29387744813e+000 7.717111714
5、24e+0008.22522001405e+000 8.64866606519e+000 9.25420034458e+0004. 讨论初始向量对计算结果的影响如果初始迭代向量选取使得其在特征值上的分量为零，但是由于计算中存在误差，迭代次数充分大后理论上使得在特征值上的分量不为零于是以为新的初始迭代向量仍可以得到结果。当增加迭代初始向量中非零元素的个数时，计算出的结果会越接近准确值。实际上，当初始向量中的0元素较多时，可能=0的情况较为普遍，许多a都有可能等于0，此时计算出的结果便与最大特征值差距较大。这是因为实际中由于精度的要求和算法允许的计算量无法做到迭代次数充分大，于是在实际选取中初始迭
6、代向量最好使得其在按模最大特征值的分量上不为零。附录A#include /*编译环境gcc或VS2012均可通过*/#include #include #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
7、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 uk40,lambda_uk40,y_uk40502;double condA2,detA;double *y_1,*y_s,*y
8、_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+)Aij=-0.064;else if(i = 2 | i = 4)
9、for(j=1;j=n;j+)Aij=0.16;elsefor(j=1;j=n;j+)Aij=Bj=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);A11=A12=A21=A4501=A5500=A5501=0;/*求lambda_s,用反幂法求按模最小的特征值*/lambda_s=fanmifa(A,r,s,n,y_s);printf(lambda_s=%.11en,lambda_s);/*用幂法求按模最大的特征值，所求为lambda_1或lambda_501*/lambda_1=mifa(A,r,s,n,y_1);/*带原点平移的幂法*/for(j=1;
10、j=n;j+)As+1j=Bj-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=%.11enlambda_501=%.11enn,lambda_1,lambda_501);/*求和Uk最接近的特征值*/for(k=
11、1;k=39;k+)ukk=lambda_1+k*(lambda_501-lambda_1)/40;for(j=1;j=n;j+)As+1j=Bj-ukk;lambda_ukk=fanmifa(A,r,s,n,y_ukk);lambda_ukk+=ukk;/*加上原点的偏移量*/for(i=1;i=39;i+)printf(lambda_u%d=%.11en,i,lambda_uki);printf(n);/*求A的(谱范数)条件数cond(A)2*/condA2=fabs(lambda_1/lambda_s);if(condA2 1)condA2=1/condA2;printf(cond(A
12、)2=%.11en,condA2);/*求A的行列式的值*/for(i=1;i=n;i+)/*恢复矩阵A*/As+1i=Bi;if(det(A,r,s,n,&detA)printf(det(A)=%.11en,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)
13、*说 明：使用的是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#这里要讨论初始向量对结果的影响#-*/*(u+i)=1;/*幂法迭代*/k=0;beta2=0;dok+;beta1=beta2;max=0;for(i=1;i=n;i+)/*求出向量的2范数*/max+=pow(ui,2);max=sqrt
