北航数值分析B第一次上机作业算法作业.docVIP

一、算法方案设计 （1）矩阵A: 矩阵A501×501为大型带状矩阵，上半带宽S=2，下半带宽R=2，可将其用循环语句存储在一个5行501列的二维数组A[][501]中，使得矩阵A的第j列存放数组A 的第j列带状元素，并使得矩阵的主对角元素存放在数组的第三行。检索时只要按照aij=数组a[i-j+3][j]即可。 （2）求λ1，λ501，λs： 用幂法求出模最大特征，再对其判断正负，如果是正数，则该特征值为λ501，如果是负数，则该特征值为λ1。 第二步，用反幂法求λS。 由于λ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、求行列式： 在用反幂法时需要对矩阵A进行三角分解，A=LU，根据det(A)=det（LU）=det（L）*det（U），det（）=1，det(A)即为矩阵的对角线元素乘积。 b、求A的谱范数条件数cond（A）2: 由于A是非奇异实对称阵，从而cond（A）2 =∣λ1∣/∣λs∣。取初始向量为 u_0[501]={1,1…1}，计算结果截图如下： （2）讨论初始迭代向量取值对计算结果的影响 在编程实现算法的过程中遇到了很多问题，多次尝试才得以解决。例如，对各个函数的定义后需要对矩阵A进行初始化；为简化程序，方便改变变量值，应该尽可能地将某些多级变量写成函数的形式，只需要对初级变量赋值即可。 猜想：幂法迭代的精度跟初始向量 u_0[501]的取值有关。 ①取初始向量为[1,1,…,0]T（共50个1），计算结果如下： 对比初始向量为 [1,1,…,1]T 计算结果可发现，最小特征值和按模最小特征值均发生较大变化。 ②取初始向量为[1,1,…,0]T（共100个1），计算结果如下： 计算结果与初始向量为 [1,1,…,1]T 时的计算结果基本一样。 ③取初始向量为[1,1,…,0]T（共300个1），计算结果如下： 对比初始向量为 [1,1,…,1]T 计算结果可发现，最大、最小特征值以及按模最小特征值都发生变化，其中最小特征值和按模最小特征值相较初始向量中含0较多的①②情况更接近初始向量为 [1,1,…,1]T 的计算结果，但最大特征值变化比较明显。 ④取初始向量为[10,10,…,10]T（元素全为10），计算结果如下： 计算结果跟初始向量为 [1,1,…,1]T 时的计算结果一样，没有变化。 ⑤取初始向量为[0,0,…,0]T（元素全为0），计算结果如下： 计算出错。 总结： 通过以上五组计算结果对比，可以初步发现： 采用幂法迭代时，选取的初始向量中0元素的多少对于计算结果和运算次数有很大的影响0元素越少，结果越准确。由于对初始向量进行了单位化，所以初始向量的模的大小对于运算结果没有太大影响。α1=0，由于计算过程中的舍入误差影响，必然会使得迭代过程中的某一步出现x1方向分量不为零的情况，相当于重新进行迭代，增大计算量和误差。（参考课本49页） 附：C程序代码 #includestdio.h #includemath.h #define N 502 /*数组记法：从a1,1]到a501,501]*/ #define accuracy 1.0e-12 #define r 2 #define s 2 double c[6][N]; /*定义c矩阵存储压缩后的带状矩阵*/ double fuzhi(); /*赋值函数*/ void (); /*Doolittle三角分解程序*/ int max(int a,int b); /*求两个数较大值函数*/ int min(int a,int b); /*求两个数较小值函数*/ double mifa(); /*幂法计算*/ double fanmifa(); /*反幂法计算*/ double fuzhi() /*定义赋值程序，将带状函数逆时针旋转°后，按行赋值，行从到，列从到，未赋值区域值为*/ { int i,j; i=1; for(j=3;jN;j++) {c[i][j]=-0.064;} i=2; for(j=2;jN;j++) {c[i][j]=0