北航数值分析第一次大作业.doc

一、算法的设计方案：（一）各所求值得计算方法1、最大特征值501，最小特征值1，按模最小特征值s的计算方法 首先使用一次幂法运算可以得到矩阵的按模最大的特征值，必为矩阵A的最大或最小特征值，先不做判断。对原矩阵A进行一次移项，即（A-I），在进行一次幂法运算，可以得到另一个按模最大特征值0。比较和0的大小，较大的即为501，较小的即为1。 对矩阵A进行一次反幂法运算，即可得到按模最小特征值s。2、A与k值最接近的特征值ik的计算方法首先计算出k所对应的k值，对原矩阵A进行一次移项，即（A-kI），得到一个新的矩阵，对新矩阵进行一次反幂法运算，即可得到一个按模最小特征值i。则原矩阵A与k值最接近的特征值 ik=i+k。3、A的（谱范数）条件数cond（A）2的计算方法condA2=1n 其中1和n分别为矩阵A的按模最大和按模最小特征值。（二）程序编写思路。 由于算法要求A的零元素不存储，矩阵A本身为带状矩阵，所以本题的赋值，LU分解，反幂法运算过程中，均应采用Doolittle分解法求解带状线性方程组的算法思路。 幂法、反幂法和LU分解均是多次使用，应编写子程序进行反复调用。二、源程序：#include#include#include #include #include#include /*头文件*/ /*定义全局变量*/ #define N 502 /*取N为502，可实现从1到501的存储，省去角标变换的麻烦*/ #define epsilon 1.0e-12 /*定义精度*/ #define r 2 /*r，s为带状矩阵的半带宽，本题所给矩阵二者都是2*/ #define s 2doublec6N; /*定义矩阵存储压缩后的带状矩阵*/ double fuzhi(); /*赋值函数*/void LUfenjie();/*LU分解程序*/int max(int a,int b); /*求两个数字中较大值*/int min(int a,int b); /*求两个数字中较小值*/double mifa();/*幂法计算程序*/double fanmifa(); /*反幂法计算程序*/double fuzhi() /*赋值程序，按行赋值，行从1到5，列从1到501，存储空间 的第一行第一列不使用，角标可以与矩阵一一对应，方便书写程序*/int i,j;i=1;for(j=3;jN;j+)cij=-0.064;i=2;for(j=2;jN;j+)cij=0.16;i=3;for(j=1;jN;j+)cij=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);i=4;for(j=1;jN-1;j+)cij=0.16;i=5;for(j=1;jb)?a:b);int min(int a,int b) return(ab)?a:b);void LUfenjie() /*LU分解程序，采用的是带状矩阵压缩存储后的LU分解法*/ double temp; int i,j,k,t; for(k=1;kN;k+) for(j=k;j=min(k+s,N-1);j+) temp=0; for(t=max(1,max(k-r,j-s);t=(k-1);t+) temp=temp+ck-t+s+1t*ct-j+s+1j; ck-j+s+1j=ck-j+s+1j-temp; for(i=k+1;i=min(k+r,N-1);i+) temp=0; for(t=max(1,max(i-r,k-s);t=(k-1);t+) temp=temp+ci-t+s+1t*ct-k+s+1k; ci-k+s+1k=(ci-k+s+1k-temp)/cs+1k; double mifa() /*幂法计算程序*/double u0N,u1N;double temp,Lu,beta=0,beta0; int i,j; for(i=1;iN;i+) /*选取迭代初始向量*/u0i=1; do beta0=beta; temp=0; for(i=1;iN;i+) temp=temp+u0i*u0i; Lu=sqrt(temp); for(i=1;iN;i+) u1i=u0i/Lu; for(i=1;iN;i+) temp=0; for(j=max(i-1,1);j=min(i+2,N-1);j+) temp=temp+ci-j+s+1j*u1j; u0i=temp; temp=0; for(i=1;i=epsilon); /*迭代运行条件判断*/return(beta);double fanmifa() /*反幂法计算程序*/double u0N,u1N,u2N;double temp,Lu,beta=0,beta0; int i,t; for(i=1;iN;i+) /*选取迭代初始向量*/u0i=1;do beta0=beta; temp=0; for(i=1;iN;i+) temp=temp+u0i*u0i; Lu=sqrt(temp); for(i=1;iN;i+)u1i=u0i/Lu;u2i=u1i; fuzhi(); LUfenjie(); /*带状矩阵压缩存储并进行LU分解后，求解线性方程组得到迭代向量uk，即程序中的u0*/ for(i=2;iN;i+) temp=0;for(t=max(1,i-r);t=1;i-)temp=0;for(t=i+1;t=min(i+s,N-1);t+)temp=temp+ci-t+s+1t*u0t;u0i=(u2i-temp)/cs+1i;temp=0; for(i=1;i=epsilon); /*迭代运行条件判断*/return(beta);void main() double u40; /*定义数组，存放k值运算得到的k值*/double lambda1,lambda501,lambdak,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/fuzhi();for(j=1;jb) /*比较两个按模最大值大小，并相应输出最大特征值501和最小特征值1*/ lambda1=b; lambda501=a; printf(矩阵A最小的特征值lambda1=%13.11en,lambda1); printf(矩阵A最大的特征值lambda501=%13.11en,lambda501); else lambda1=a; lambda501=b; printf(矩阵A最小的特征值lambda1=%13.11en,lambda1); printf(矩阵A最大的特征值lambda501=%13.11en,lambda501); printf(矩阵A按模最小特征值lambdas=%13.11en,d); /*输出按模最小特征值s*/ for(k=1;k40;k+) /*对每一个进行移项反幂法运算，求出最接近k的特征值并输出*/ uk=(lambda501-lambda1)*k/40+lambda1; fuzhi(); for(j=1;jN;j+) c3j=c3j-uk; lambdak=fanmifa()+uk; i=k; printf(矩阵A最接近uk特征值lambdak%d=%13.11en,i,lambdak); cond=fabs(a/d); printf(A的条件数=%13.11en,cond); /*计算A条件数并输出*/ fuzhi(); /*计算A的行列式值并输出*/ LUfenjie(); det=1; for(i=1;iN;i+) det=det*c3i; printf(行列式的值detA=%13.11en,det); 三、程序的运行结果：四、初始向量的选取对计算结果的影响：（一）选取形式不变，数值变换1、取u0为0.5,0.5.0.5，运行结果如下：2、取u0为50,50.50，运行结果如下： 从运行结果来看，此类初始向量的选取对结果不会产生影响，即使选成0，结果也不变化。（二）改变选取形式取u0为1,0.0，运行结果如下：可以看出，改变初始向量的类型之后，运行结果出现了明显的变化，说明迭代初始向量的选取对迭代运算的结果是有影响的。为了观察初始向量与运行结果之间的联系，下面以幂法为例进行试验。1、取u0为1,0.0，运行幂法程序，运行结果为：对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第22个位置，且第1位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：2、取u0为01.0，运行幂法程序，1位于第200个位置运行结果为：对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第198个位置，且第200位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：3、取u0为01.0，运行幂法程序，1位于第300个位置运行结果为：对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第288个位置，且第300位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：4、取u0为01.0，运行幂法程序，1位于第400个位置运行结果为：对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第413个位置，且第400位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：5、取u0为01.0，运行幂法程序，1位于第500个位置运行结果为：对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第477个位置，且第500位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：综上可以看出，迭代初始向量的选取对迭代结果有明显的影响。输出的最终迭代向量（即矩阵计算出特征值所对应的特征向量）与所选的迭代初始向量形式上很接近。因此可以总结如下：当所选取的迭代初始向量与矩阵的特征向量uk最为接近时，那么，幂法迭代计算结果输出的应该就是特征向量uk对应所对应的特征值。以上试验只是考虑了初始向量u0为01.0的情况，如果