北航研究生数值分析编程大作业1.docx

1、北航研究生数值分析编程大作业1数值分析大作业1、算法设计方案1、矩阵初始化 矩阵的下半带宽r=2,上半带宽s=2，设置矩阵，在矩阵C中检索矩阵A中的带内元素的方法是：。这样所需要的存储单元数大大减少，从而极大提高了运算效率。2、利用幂法求出幂法迭代格式： 当时，迭代终止。 首先对于矩阵A利用幂法迭代求出一个，然后求出矩阵B，其中(为单位矩阵），对矩阵B进行幂法迭代，求出，之后令，比较，大者为，小者为。3、利用反幂法求出反幂法迭代格式： 当时，迭代终止，。 每迭代一次都要求解一次线性方程组，求解过程为：（1）作分解对于执行（2）求解（数组b先是存放原方程组右端向量，后来存放中间向量y) 使用反幂

2、法，直接可以求得矩阵按模最小的特征值。 求与数最接近的特征值，对矩阵实行反幂法，即可求出对应的。4、求出A的条件数和行列式 根据，其中分子分母分别对应按模最大和最小的特征值。的计算：由于,其中为下三角矩阵，且对角线元素为1，故，所以有，又为上三角矩阵，故为对其对角线上各元素的乘积，最后可得。2、程序源代码（1）定义所需要的函数：#include #include #include #define N 501#define R 2#define S 2int min(int a,int b); / 求最小值int max(int a,int b,int c); / 求最大值double Fan_

3、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);/幂法 double InversePowerMethod(double CN,double uN,double yN,double bta,double D);/反幂法;（2）程序的主函数，Main.cpp代码如下：

4、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.0e-12; int i, j, k; FILE * fp; fp = fopen(Namda.txt,w); /对数组进行初始化/ int i, j; for (i = 0; i N; i+) ui = 1; f

5、or (i = 0;i R + S + 1;i+) for (j = 0;j N;j+) if (i=0|i=4) Cij=-0.064; else if (i=1|i=3) Cij=0.16; else if (i=2) Cij=(1.64-0.024*(j+1)*sin(0.2*(j+1) -0.64*exp(0.1/(j+1); /幂法求Namda1/ Namda1 = PowerMethod(C, u, y, bta, D); printf(n=n); printf(Namda1 = %12.11e, Namda1); printf(n=n); /幂法求Namda501/ bta =

6、1.0; for (i = 0; i R + S + 1; i+) for (j = 0; j N; j+) if (i = 2) C1ij = Cij -Namda1; else C1ij = Cij; Namda501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1; printf(n=n); printf(Namda501 = %12.11e, Namda501); printf(n=n); /反幂法求NamdaS/ bta = 1.0; NamdaS = InversePowerMethod(C, u, y, bta, D); prin

7、tf(n=n); printf(NamdaS = %12.11e, NamdaS); printf(n=n); /反幂法求Namdak/ printf(n=n); for (k = 0; k 39; k+) miuk = Namda1 + (k + 1) * (Namda501 - Namda1) / 40.0; bta = 1.0; for (i = 0; i R + S + 1; i+) for (j = 0; j N; j+) if (i = 2) C1ij = Cij - miuk; else C1ij = Cij; Namdak = InversePowerMethod(C1, u,

8、 y, bta, D) + miuk; fprintf(fp,与%12.11e最接近的特征值为:%12.11en,miuk,Namdak); printf(求与miuk最接近的Namdak的计算结果已经输出到文件Namda.txt中); printf(n=n); /求A的谱范数/ printf(n=n); printf(A的谱范数为:%12.11e, sqrt(Namda501); printf(n=n); /求A的条件数/ CondA2 = fabs( Namda1 / NamdaS); printf(n=n); printf(A的谱范数的条件数Cond(A)2为:%12.11e,CondA

9、2); printf(n=n); /求det(A)2的值/ for (j = 0; j N; j+) detA *= C2j; printf(n=n); printf(行列式A的值为:%12.11e,detA); printf(n=n); fclose(fp); _getch(); return;（3）成员函数的实现int min(int a,int b) return a b ? a : b; return temp c ? temp : c;double Fan_two(double xN) double sum = 0.0; int i; for (i = 0; i N; i+) sum

10、 += 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 = 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) br

11、eak; 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+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 =

12、 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 - 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; /用到临时存储数组TC和ty是因为函数Solve执行过程

13、中会改变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, 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;3、程序运

14、行结果下图为主程序运行结果其中的结果输出在Namda.txt文件中，结果如下：四、分析迭代初始向量对计算结果的影响选择不同的初始向量可能会得到不同的特征值。选取时，运行结果如下：选取时，运行结果如下：选取时（i=int(N/2)时为0），运行结果如下：选取时（i=int(N/2)时为1），运行结果如下：通过以上类似的实验可以大致看出这样的规律：的值趋近于有两种情况：（1）当的元素中，1的个数较多时;（2）在1的个数相同的条件下，1的分布越靠中后段，观察对应的特征向量可以发现：（1）随着i的增加，特征向量元素的绝对值不断增大，即绝对值较大的数集中于中后位置。因此，如果初始向量的非零元素集中在中后

15、段，该初始向量会更容易逼近对应的特征向量，得到的结果也越准确。对于，初始向量的非零元素集中在前半段的情况进行实验，会发现当算法中不考虑给定的精度水平，强制性执行足够高次数（大约在300多次以上）的迭代，运算结果也会趋近于。这就说明，程序之前没有得到准确结果的原因，是因为迭代次数不够。当迭代次数在100到200次左右时，每一次迭代所造成的相对误差小于给定的精度水平，因此，如果由精度水平来控制循环迭代的次数，程序将错误地判断已经收敛，但实际上，当继续迭代到300次以上时，运算结果会突然变化，直至最终稳定在。由此，可以得出结论，当迭代次数足够高（300次以上）时，得到的结果会趋于稳定，不同的初始向量和选定的精度水平，决定着程序是否出现以及何时出现假收敛。当所选取初始向量的非零元素越多，以及非零元素的位置越靠后时，收敛会更加迅速、准确。