北航数值分析计算实习报告一.doc

北京航空航天大学数值分析计算实习报告第一大题学 院：自动化科学与电气工程学院专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日 北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有1.求，和的值。2.求A的与数最接近的特征值。3.求A的(谱范数)条件数和行列式detA。说明：1. 在所用的算法中，凡是要给出精度水平的，都取。2. 选择算法时，应使矩阵A的所有零元素都不储存。3. 打印以下内容：（1） 全部源程序；（2） 特征值以及的值。4.采用e型输出实型数，并且至少显示12位有效数字。一、算法设计方案1、求，和的值。由于，可知绝对值最大特征值必为和其中之一，故可用幂法求出绝对值最大的特征值，如果=0，则=，否则=。将矩阵A进行一下平移： （1）对用幂法求出其绝对值最大的特征值，则A的另一端点特征值或为+。为按模最小特征值，可对A使用反幂法求得。2、求A的与数最接近的特征值。计算-，其模值最小的值对应的特征值与最接近。因此对A进行平移变换： （2）对用反幂法求得其模最小的特征值，则=+。3、求A的(谱范数)条件数和行列式detA。由矩阵A为非奇异对称矩阵可得： （3）其中为按模最大特征值，为按模最小特征值，通过第一问我们求得的和可以很容易求得A的条件数。在进行反幂法求解时，要对A进行LU分解得到。因L为单位下三角阵，行列式为1，U为上三角阵，行列式为主对角线乘积，所以A的行列式等于U的行列式，为U的主对角线的乘积。二、 算法实现1、 矩阵存储原矩阵A为一个上、下半带宽都为2的501501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量，A的带外元素不给存储，值存储带内元素，如下C矩阵所示： （4）C是一个5501的矩阵，相比A大大节省了存储空间，在数组C中检索矩阵A的带内元素的方法是： （5）2、 幂法 幂法迭代公式如下： （6）其中，不断迭代当时即可认为其满足精度要求，令。在程序中计算时，根据A矩阵的特点，简化如下： （7）3、反幂法反幂法迭代公式如下： （8）当k足够大时，。在求解时，可先对A进行Doolittle分解，由于A是带状结构，所以分解出的L、U也是带状结构，利用C矩阵进行Doolittle分解并求向量的算法如下：（1）作分解A=LU 对于执行： （9） 由于C语言中数组下标是从0开始的，所以在程序中矩阵元素c的下标都减1。（2） 求解（数组b先是存放原方程组右端向量，后来存放中间向量y,在程序中b和y都保存在数组y501中。） （10）求出后，其他部分与幂法求解相同。3、 结果分析实验表明，本程序中，初始向量对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。（实验结果截图见附录）迭代次数迭代次数159381160381535573350592674611311611304611343611343611343611表1 不同初始向量对应的和及其迭代次数由表1可以得到如下结论：1. 不同的初始向量对本程序的影响大，对没有影响，都能保证收敛到正确值。2. 初始向量中必须保证中至少有一个为1才能保证收敛到正确值。3. 初始向量非零值的多少和大小对迭代次数并没有明显影响。4. 为解决初始向量对程序的影响，可以先对A做平移变换再求。4、 实验程序#include stdio.h#include math.h#include stdlib.hstatic double b=0.16,c=-0.064;#define Precision 1e-12void copy(double b501,double y501);double dianji(double x,double y);/计算两个向量内积void InitMatrix(double *p);/Init Matrix Adouble NeiJi(double a,double b);/get 2范数void get_y(double *y,double *u);/get yvoid get_u(double *u,double *y,double *a);/get uvoid Initu(double *p);/初始化初始向量udouble Get_Fabs_Eigenvalue(double *a,double *u,int *iterations);/循环迭代得到绝对值最大特征值void A_sub_minI(double *a,double min);/A-min*Ivoid InitC(double C5501);/初始化数组Cvoid DoolittleC(double C5501,int n,int s,int r);/进行Doolittle分解int min(int a,int b);/取返回a,b最小值int max(int a,int b);/求最大值void Doolittle_getx(double C5501,double y501,double u501,int n,int s,int r);double Get_min_Eigenvalue(double C5501,double *u,int n,int s,int r,int *iterations);double detA(double C5501);/求A的行列式struct FinalValuedouble min;/特征值最小值double max;/特征值最大值double abs_min;/模最小特征值double abs_max;/模最大特征值double detA;/A的行列式double cond2;/A的条件数int min_iterations;/最小值迭代次数int max_iterations;/最大值迭代次数int abs_min_iterations;/求绝对值最小的迭代次数;int main()FinalValue main_num= 0,0,0,0;double temp;/两值交换中间变量int temp1;double u501=0;/,y501;/为u0赋初值double Norm_u=0;/范数double C5501 = 0;InitC(C);/初始化Cint i=0,*iterations;Initu(u);iterations = main_num.min_iterations;main_num.min = Get_Fabs_Eigenvalue(C2,u,iterations);/将求得的绝对值最大特征值放到min变量中main_num.abs_max = main_num.min;A_sub_minI(C2,main_num.min);for(i=0;i501;i+)ui=i+1;iterations = main_num.max_iterations;main_num.max = Get_Fabs_Eigenvalue(C2,u,iterations);/将求得的绝对值最大特征值放到min变量中main_num.max += main_num.min;if(main_num.minmain_num.max)temp = main_num.min;main_num.min = main_num.max;main_num.max = temp;temp1 = main_num.min_iterations;main_num.min_iterations = main_num.max_iterations;main_num.max_iterations = temp1;printf(最小特征值为：%.12e,迭代次数为：%dn,main_num.min,main_num.min_iterations);printf(最大特征值为：%.12e,迭代次数为：%dn,main_num.max,main_num.max_iterations);/*/*以下利用反幂法求解模最小的特征值*/*/InitC(C);/初始化CInitu(u);DoolittleC(C,501,2,2);main_num.detA = detA(C);iterations = main_num.abs_min_iterations;main_num.abs_min = Get_min_Eigenvalue(C,u,501,2,2,iterations);main_num.cond2 = fabs(main_num.abs_max/main_num.abs_min);printf(绝对值最小特征值为：%.12e,迭代次数为：%dn,main_num.abs_min,main_num.abs_min_iterations);printf(A的行列式为：%.12e；,main_num.detA);printf(A的条件数cond(A)2为：%.12en,main_num.cond2);/*/*以下利用反幂法求与数列Uk中元素最相近的特征值*/*/double Uk39;/保存Uk的值并保存与Uki最接近的double B39;int diedai;iterations = diedai;for( i=1;i=39;i+)Uki-1 = main_num.min + i*(main_num.max - main_num.min)/40;InitC(C);Initu(u);for(int j =0;j501;j+)C2j -= Uki-1;DoolittleC(C,501,2,2);Bi-1 = Get_min_Eigenvalue(C,u,501,2,2,iterations);Bi-1 +=Uki-1;printf(%-2d=%-3.12e,与%-2d最相近的特征值=%-3.12en,i,Uki-1,i,Bi-1);return 0;double detA(double C5501)int i =0;double e =1;for(i =0;i501;i+)e*=C2i;return e;void InitMatrix(double *p)for(int i=1;i=501;i+)*p=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);p+;void Initu(double *p) for(int i=1;i=501;i+) if(i=500)*p=0; else *p=1;p+; / *p=1;void A_sub_minI(double *a,double min)for(int i=1;i=501;i+)*a=*a-min;a+;double NeiJi(double *a,double *b)double e=0.0;for(int i=0;i501;i+)e=e+(*a)*(*b);a+;b+;return e;void get_y(double *y,double *u)double Norm_u=sqrt(NeiJi(u,u);for(int i=0;i501;i+)*y = *u/Norm_u;y+;u+;void get_u(double *u,double *y,double *a)u0=a0*y0+b*y1+c*y2;u1=b*y0+a1*y1+b*y2+c*y3;for(int i=2;i499;i+)ui=c*yi-2+b*yi-1+ai*yi+b*yi+1+c*yi+2;u499=c*y497+b*y498+a499*y499+b*y500;u500=c*y498+b*y499+a500*y500;void InitC(double C5501)int i;for(i = 2;i 501;i+)C0i = -0.064;C4i-2 = -0.064;for(i = 1;i 501;i+)C1i = 0.16;C3i-1 = 0.16;for(i = 1;i = 501;i+)C2i-1 = (1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);double Get_Fabs_Eigenvalue(double *a,double *u,int *iterations)double y501,B_k0=0,B_k1=0;double wucha;int i=0;while(1)+i;get_y(y,u);get_u(u,y,a);B_k1=NeiJi(y,u);/是否判断B_K1是否为0？wucha=(fabs(B_k1-B_k0)/(fabs(B_k1);/get wuchaif(wucha=Precision)break;else if(i10000)printf(迭代次数超长，请更改初始向量n);break;else B_k0 = B_k1;*iterations = i;return B_k1;double Get_min_Eigenvalue(double C5501,double *u,int n,int s,int r,int *iterations)double y501 = 0,B_k0=0,B_k1=0;double b501 = 0;/储存y值的中间向量double wucha;int i=0;while(1)+i;get_y(y,u);copy(b,y);/保护y向量Doolittle_getx(C,y,u,501,2,2);copy(y,b);B_k1=NeiJi(y,u);/是否判断B_K1是否为0？wucha=(fabs(1/B_k1-1/B_k0)/(1/fabs(B_k1);/get wucha/wucha=(fabs(B_k1-B_k0)/(fabs(B_k1);/get wuchaif(wucha=Precision)break;else if(i10000)printf(迭代次数超长，请更改初始向量n);break;else B_k0 = B_k1;*iterations = i;return (1/B_k1);void Doolittle_getx(double C5501,double y501,double u501,int n,int s,int r)int i,t;double e;for(i = 2;i = n;i+)e = 0;for(t = max(1,i-r);t = i-1;t+)e+=Ci-t+s+1-1t-1*yt-1;yi-1 -= e;un-1 = yn-1/Cs+1-1n-1;for(i = n-1;i=1;i-)e = 0;for(t = i+1;t