1、 数值分析第一次上机作业 一、算法方案设计 （ 1） 存储 矩阵 A（参考课本 25 页 ） : 矩阵 A501 501为大型带状矩阵，上半带宽 S=2，下半带宽 R=2， 参照课本，可将其用循环语句存储在一个 5 行 501 列的二维数组 A5501中，使得矩阵 A 的第 j 列存放数组 A 的第 j 列带状元素，并使得矩阵的主对角元素存放在数组的第三行。检索 矩阵的元素 时只要按照 公式： aij=数组 ai-j+3j即可。 （ 2）求 1， 501， s： 第一步， 用幂法求出 矩阵 A按 模最大特征 值 ，再对其判断正负，如果是正数，则该特征值为 501 ， 如果是负数，则该特征值 为
2、 1。 幂法实现过程具体为： 第二步，用反幂法 求 矩阵 A按模最小特征向量 S。 班级： ZY16131 学号： ZY16131 姓名： 数值分析第一次上机作业 第三步， 由于 1 2 . 501， 可 采用原点平移法对矩阵 A 平移 1，得到矩阵 （ - 1I+A）， 记为矩阵 B， 再用 用幂法求出 B 的模最大特征值 B，则 501= B+ 1。 （ 3） 距离 k最近的特征值： 还是用原点平移法，将矩阵 A 平移 k个单位，再用反幂法求出平移后矩阵模最小特征值 k，矩阵 A 与 k 最接近的特征值 为： ik = k + k = k + （ 4） A 的谱范数条件数与 detA： a
3、、求 矩阵 A的 行列式 值 ： 在用反幂法时需要对矩阵 A进行 Doolittle三角分解， A=LU，根据 det(A)=det（ LU） =det（ L） *det（ U）， 其中 det（ L） =1， det(A)即为 U矩阵的对角线元素 的 乘积。 数值分析第一次上机作业 b、求 A 的 （ 谱范数 ） 条件数 cond（ A） 2: 由于 A 是非奇异实对称阵，从而 cond（ A） 2 = 1/ s 。 二、 运行结果分析 （ 1） 取初始向量为 u_0501=1,1 1， 计算结果截图如下： 数值分析第一次上机作业 （ 2）讨论初始迭代向量取值对计算结果的影响 在编程实现算法
4、的过程中遇到了很多问题，多次尝试才得以解决。例如，对各个函数的定义后需要对矩阵 A 进行初始化；为简化程序，方便改变变量值，应该尽可能地将某些多级变量写成函数的形式，只 需要对初级变量赋值即可。 猜想：幂法迭代的精度跟初始向量 u_0501的取值有关。 取初始向量为 1,1, ,0T（共 50 个 1），计算结果如下： 对比初始向量为 1,1, ,1T 计算 结果可发现，最小特征值和按模最小特征值均发生较大变化。 取初始向量为 1,1, ,0T（共 100个 1），计算结果如下： 计算结果与 初始向量为 1,1, ,1T 时的计算结果基本一样。 取 初始向量为 1,1, ,0T（共 300个
5、1），计算结果如下： 对比 初始向量为 1,1, ,1T 计算结果可发现，最大、最小特征值以及按模最小特征值都发生变化，其中最小特征值和按模最小特征值相 较 初始向量中含 0较多的情况更接近 初始向量为 1,1, ,1T 的计算结果，但最大特征值变化比较明显。 数值分析第一次上机作业 取初始向量为 10,10, ,10T（元素全为 10），计算结果如下： 计算结果跟 初 始向量为 1,1, ,1T 时的计算结果一样，没有变化。 取初始向量为 0,0, ,0T（元素全为 0），计算结果如下： 计算出错。 总结： 通过以上五组计算结果对比，可以初步发现： 采用幂 法迭代时，选取的 初始向量中 0
6、元素的多少对于计算结果和运算次数有很大的影响 。基本规律是： 0 元素越少，结果越准确。由于对初始向量进行了单位化，所以初始向量的模的大小对于运算结果没有太大影响。 1 1 2 2 3 30u x x x 1221221112211nknkknknkkkxxxxAxAxAu这是因为，如果所选的 u0使得 1=0，由于计算过程中的舍入误差影响，必然会使得迭代过程中的某一步出现 x1 方向分量不为零的情况，相当于重新进行迭代，增大计 算量和误差。（参考课本 49 页） 数值分析第一次上机作业 附： C 程序代码 #include #include #define N 502 /*数组记法：从 a1
7、,1到 a501,501*/ #define accuracy 1.0e-12 #define r 2 #define s 2 double c6N; /*定义 c矩阵存储压缩后的带状矩阵 */ double fuzhi(); /*赋值函数 */ void Doolittle(); /*Doolittle三角 分解程序 */ int max(int a,int b); /*求两个数较大值函数 */ int min(int a,int b); /*求两个数较小值函数 */ double mifa(); /*幂法计算 */ double fanmifa(); /*反幂法计算 */ double f
8、uzhi() /*定义 赋值程序，将带状函数逆时针旋转 45后，按行赋值，行从 1到 5，列从 1到 501，未赋值区域值为 0*/ int i,j; i=1; for(j=3;jb)?a:b); int min(int a,int b) /*定义求最小值函数 */ return(a=accuracy); /*判断精度水平是否满足，不满足则继续迭代 */ return(beta); /*满足精度要求之后返回 beta值 */ double fanmifa() /*反幂法计算程序 */ double u0N,u1N,u2N; double temp,yita,beta=0,beta0; int
9、i,t; for(i=1;i=1;i-) temp=0; for(t=i+1;t=accuracy); /*迭代运行条件判断 */ return(beta); main() /*主函数 */ double u40; /*定义数组，用于计算第 2问 */ double lambda_1,lambda_501,lambda_s,a,b,d,cond,det; int i,j,k; fuzhi(); a=mifa(); /*幂法计算按模最大值 */ fuzhi(); for(j=1;jb) /*比较两个按模最大值大小，并相应输出最大特征值 和最小特征值 */ lambda_1=b; lambda_5
10、01=a; printf(“Lmd_1=%13.11en“,lambda_1); printf(“Lmd_501=%13.11en“,lambda_501); else lambda_1=a; lambda_501=b; printf(“Lmd_1=%13.11en“,lambda_1); printf(“Lmd_501=%13.11en“,lambda_501); fuzhi(); d=fanmifa(); /*反幂法计算按模最小值 */ printf(“Lmd_s=%13.11en“,d); /*输出按模最小特征值 s*/ printf(“n“); for(k=1;k40;k+) /*对每
11、一个进行移项反幂法运算，求出最接近 k的特征值并输出 */ 数值分析第一次上机作业 uk=(lambda_501-lambda_1)*k/40+lambda_1; fuzhi(); for(j=1;jN;j+) c3j=c3j-uk; lambda_s=fanmifa()+uk; i=k; printf(“u%d=%13.11e Lmd_i%d=%13.11en“,i,uk,i,lambda_s); printf(“n“); cond=fabs(a/d); printf(“cond(A)2=%13.11en“,cond); /*计算 A条件数并输出 */ fuzhi(); /*计算 A的行列式值并输出 */ Doolittle(); det=1; for(i=1;iN;i+) det=det*c3i; printf(“detA=%13.11en“,det); printf(“n“);
