1、1 数值分析大作业 1、算法设计方案 1、矩阵初始化 矩阵 的下半带宽 r=2,上半带宽 s=2，设置矩阵 ，501ijaA 501srC 在矩阵 C 中检索矩阵 A 中的带内元素 的方法是： 。这样所需要ijajsiijca,1 的存储单元数大大减少，从而极大提高了运算效率。 2、利用幂法求出 501与 幂法迭代格式： 0111nTkkkkTuRyAu非 零 向 量 当 时，迭代终止。120/kk 首先对于矩阵 A 利用幂法迭代求出一个 ，然后求出矩阵 B，其中 ( 为单位矩阵） ，对矩阵 B 进行幂法迭代，求出 ，之后令IB ，比较 ，大者为 ，小者为 。与5011 3、利用反幂法求出 i
2、ks, 反幂法迭代格式： 0111nTkkkTkuRyAu非 零 向 量 当 时，迭代终止， 。120/kk 1sk 每迭代一次都要求解一次线性方程组 ，求解过程为：1yAu （1）作分解 LUA 对于 执行nk,.2 2 sknrki cccsj ksktkittksisi jstksjrttjskjsk );,m(,.21/)(:),.: ,1,1ax,1, ,11),max(,1,1 （2）求解 （数组 b 先是存放原方程组右端向量，后来存放中间向yUbL 量 y) )1,.2(/)(:/ ),.32(: ,1),min(11),ax(, nicbcistsittintiritstii
3、使用反幂法，直接可以求得矩阵按模最小的特征值 。s 求与数 最接近的特征值 ，对矩阵)39,.2(40151kk ik 实行反幂法，即可求出对应的 。IAk kikk,/1 4、求出 A 的条件数和行列式 根据 ，其中分子分母分别对应按模最大和最小的特征值。max2()scond 的计算：由于 ,其中 为下三角矩阵，且对角线元素为 1，故etALU ，所以有 ，又 为上三角矩阵，故 为对其对角线d()1L det()U 上各元素的乘积，最后可得 。det()t() 3 2、程序源代码 （1）定义所需要的函数： #include #include #include #define N 501 #
4、define R 2 #define S 2 int min(int a,int b); / 求最小值 int max(int a,int b,int c); / 求最大值 double Fan_two(double xN);/计算二范数 void FenjieLU(double (*C)N);/解线性方程组的LU分解过程 void Solve(double (*C)N, double *b,double *x);/解线性方程组的求解过程 double PowerMethod(double CN,double uN,double yN,double bta,double D);/ 幂法 dou
5、ble InversePowerMethod(double CN,double uN,double yN,double bta,double D);/反幂法 ; （2）程序的主函数，Main.cpp 代码如下： void main() double CR+S+1N; double uN; double yN; double miu39; double C1R+S+1N; double bta = 1.0; double Namda1,Namda501,NamdaS; double Namda39; double CondA2; double detA = 1.0; double D = 1.0
6、e-12; int i, j, k; FILE * fp; fp = fopen(“Namda.txt“,“w“); /对数组进行初始化/ int i, j; for (i = 0; i c ? temp : c; double Fan_two(double xN) double sum = 0.0; int i; for (i = 0; i N; i+) sum += pow(xi,2); return sqrt(sum); void FenjieLU(double (*C)N) double sum = 0; int i, j, k,t; for (k = 0; k N; k+) j =
7、k; i = k + 1; while (1) if (j = min(k + S + 1, N) break; for (t = max(0, k - R, j - S); t = k - 1; t+) sum += Ck-t+St * Ct-j+Sj; Ck-j+Sj = Ck-j+Sj - sum; sum = 0.0; j+; if (k = N-1) 7 break; if (i = min(k + R + 1, N) break; for (t = max(0, i - R,k - S); t = k - 1; t+) sum += Ci-t+St * Ct-k+Sk; Ci-k+
8、Sk = (Ci-k+Sk - sum) / CSk; sum = 0; i+; void Solve(double (*C)N, double *b,double *x) double sum = 0; int i, t; sum = 0; for (i = 1; i N; i+) for (t = max(0, i - R); t = 0; i-) for (t = i+1; t D) temp = bta; ita = Fan_two(u); for (i = 0; i N; i+) yi = ui / ita; for (i = 0; i N; i+) for (j = max(0,i
9、 - R); j min(i + S + 1,N); j+) sum += Ci - j + Sj * yj; ui = sum; sum = 0; for (i = 0; i D) temp = bta; ita = Fan_two(u); for (i = 0; i N; i+) yi = ui / ita; 9 /用到临时存储数组TC和ty 是因为函数Solve执行过程中会改变 A和 y for (i = 0; i R + S + 1; i+) for (j = 0; j N; j+) TCij = Cij; for (i = 0; i N; i+) tyi = yi; Solve(C,
10、 y, u); for (i = 0; i R+S+1; i+) for (j = 0; j N; j+) Cij = TCij; for (i = 0; i N; i+) yi = tyi; for (i = 0; i N; i+) sum += yi * ui; bta = sum; sum = 0; k+; bta = 1.0 / bta; return bta; 10 3、程序运行结果 下图为主程序运行结果 其中 的结果输出在 Namda.txt 文件中，结果如下：ik 11 四、分析迭代初始向量对计算结果的影响 选择不同的初始向量 可能会得到不同的特征值。uN 选取 时，运行结果如下
11、：1,0,u 12 选取 时，运行结果如下：1,uN 选取 时（i=int(N/2)时为 0） ，运行结果1,0,uN 如下： 选取 时（i=int(N/2)时为 1） ，运行结果0,1uN 如下： 通过以上类似的实验可以大致看出这样的规律： 的值趋近于 有两种情况：11.0736501e （1）当 的元素中， 1 的个数较多时;uN （2）在 1 的个数相同的条件下，1 的分布越靠中后段， 13 观察 对应的特征向量可以发现： （1）随着 i 的增加，特征向量元素的绝对值不断增大，即绝对值较大的数 集中于中后位置。因此，如果初始向量的非零元素集中在中后段，该初始向量 会更容易逼近对应的特征向
12、量，得到的结果也越准确。 对于，初始向量的非零元素集中在前半段的情况进行实验，会发现当算法 中不考虑给定的精度水平，强制性执行足够高次数（大约在 300 多次以上）的 迭代，运算结果也会趋近于 。这就说明，程序之前没有1.0736501e 得到准确结果的原因，是因为迭代次数不够。当迭代次数在 100 到 200 次左右 时，每一次迭代所造成的相对误差小于给定的精度水平，因此，如果由精度水 平来控制循环迭代的次数，程序将错误地判断已经收敛，但实际上，当继续迭 代到 300 次以上时，运算结果会突然变化，直至最终稳定在 。1.0736501e 由此，可以得出结论，当迭代次数足够高（300 次以上）时，得到的结果 会趋于稳定，不同的初始向量和选定的精度水平，决定着程序是否出现以及何 时出现假收敛。当所选取初始向量的非零元素越多，以及非零元素的位置越靠 后时，收敛会更加迅速、准确。
