雅克比法求矩阵特征值特征向量

1、C语言课程设计报告课程名称： 计算机综合课程设计 学 院： 土木工程学院 设计题目： 矩阵特征值分解 级 别： B 学生姓名： 学 号： 同组学生： 无 学 号： 无 指导教师： 2012年 9 月 5 日C语言课程设计任务书（以下要求需写入设计报告书）学生选题说明: 以所发课程设计要求为准，请同学们仔细阅读； 本任务书提供的设计案例仅供选题参考；也可自选，但难易程度需难度相当； 鼓励结合本专业（土木工程、力学）知识进行选题，编制程序解决专业实际问题。 限2人选的题目可由1-2人完成（A级）；限1人选的题目只能由1人单独完成（B级）； 设计总体要求: 采用模块化程序设计； 鼓励可视化编程； 源

2、程序中应有足够的注释； 学生可自行增加新功能模块（视情况可另外加分）； 必须上机调试通过； 注重算法运用，优化存储效率与运算效率； 需提交源程序(含有注释)及相关文件(数据或数据库文件)； (cpp文件、txt或dat文件等) 提交设计报告书，具体要求见以下说明。设计报告格式：目录1 课程设计任务书（功能简介、课程设计要求）；2 系统设计（包括总体结构、模块、功能等，辅以程序设计组成框图、流程图解释）；3 模块设计（主要模块功能、源代码、注释（如函数功能、入口及出口参数说明，函数调用关系描述等）；4 调试及测试：（调试方法，测试结果的分析与讨论，截屏、正确性分析）；5 设计总结：（编程中遇到的

3、问题及解决方法）；6 心得体会及致谢；参考文献 1.课程设计任务书功能简介：a） 输入一个对称正方矩阵A，从文本文件读入；b） 对矩阵A进行特征值分解，将分解结果：即U矩阵、S矩阵输出至文本文件；c） 将最小特征值及对应的特征向量输出至文本文件；d） 验证其分解结果是否正确。提示：A=USUT,具体算法可参考相关文献。功能说明：矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。注：以三阶对称矩阵为例总体结构 2.系统设计递推求矩阵U递推求矩阵S找矩阵A中非对角元素中的最大值，并记下其位置从文本文件中读入数组A向屏幕和txt文件输入矩阵S

4、求最小特征值及其对应特征向量，并输出到屏幕和txt文件中向屏幕和txt文件输入矩阵Utezheng函数main函数3. 模块设计#include#include#includeint main()FILE *fp;int tezheng(double *a,int n,double *s,double *u,double eps,int itmax); /函数调用声明int i,j,p,itmax=1000; /itmax为最大循环次数double eps=1e-7,s33,u33; /eps为元素精度，s为对角矩阵S，u为矩阵U double a9;/a为待分解矩阵Ai=tezheng(a,

5、3,s,u,eps,1000);if(i0) /i对应函数中的返回值itif(fp=fopen(juzhen.txt,w)=NULL) /打开待输入txt文件printf(无法打开文件.n);return;printf(U矩阵为:n); /下几句分别向屏幕和txt文件输入矩阵Ufprintf(fp,U矩阵为:n);for(i=0;i3;i+)for(j=0;j3;j+)printf(%10.6f,uij);fprintf(fp,%10.6f,uij);printf(n);fprintf(fp,n);printf(S对角矩阵为:n); /下几句分别向屏幕和txt文件输入矩阵Sfprintf(fp

6、,S对角矩阵为:n);for(i=0;i3;i+)for(j=0;j3;j+)printf(%10.6f,sij);fprintf(fp,%10.6f,sij); printf(n);fprintf(fp,n); p=0;for(i=0;i3;i+)/下面几句为求最小特征值及其对应特征向量，并输出到屏幕和txt文件中if(siis00)p=i;printf(最小特征值为:%10.6fn,spp);fprintf(fp,最小特征值为:%10.6fn,spp);printf(对应特征向量为:n);fprintf(fp,对应特征向量为:n);for(i=0;i3;i+)printf(%10.6fn,

7、uip);fprintf(fp,%10.6fn,uip);int tezheng(double *a,int n,double s33,double u33,double eps,int itmax)FILE *A;char line1000; /存放从文件juzhen A.txt读入的数据char *k= ;int i,j,p,q,it;double sint,cost,sin2t,cos2t,d,tmp,t1,t2,t3,fm;it=0;if(A=fopen(juzhen A.txt,r)=NULL)printf(无法打开文件n);return; fgets(line,1000,A); /

8、从文件指针A指向的文件，即juzhen A.txt中读入字符数据到数组line中 strtok(line,k); /原型char *strtok(char *s, const char *delim);说明：strtok()用来将字符串分割成一个个片段。参数s指向欲分割的字符串，参数delim则为分割字符串，当strtok()在参数s的字符串中发现到参数delim的分割字符时则会将该字符改为0 字符。在第一次调用时，strtok()必需给予参数s字符串，往后的调用则将参数s设置成NULL。每次调用成功则返回被分割出片段的指针。 for(i=0;in*n;i+) ai = atof(strtok

9、(NULL,k); /atof函数功能: 把字符串转换成浮点数 fclose(A); /释放A指针for(i=0;in;i+) /以下几句初始化矩阵Uuii=1.0;for(j=0;jn;j+)if(i!=j)uij=0.0;while(ititmax)it+;d=0.0;for(i=1;in;i+) /该for循环找矩阵A中非对角元素中的最大值，并记下其位置for(j=0;jd)d=tmp;p=i;q=j;if(deps) /该if条件句用数组元素精度控制循环结束与否for(i=0;in;i+)for(j=0;jn;j+)sij=ai*n+j;return (it);sint=2*ap*n+

10、q; /以下几句为递推公式,求S矩阵主对角元素 sin2t=2*ap*n+q/(sqrt(4*ap*n+q*ap*n+q+(aq*n+q-ap*n+p)*(aq*n+q-ap*n+p);if(aq*n+q-ap*n+p0.0)sin2t=-sin2t;cos2t=sqrt(1.0-sin2t*sin2t);sint=sin2t/(sqrt(2.0*(1.0+cos2t);cost=sqrt(1.0-sint*sint);tmp=ap*n+p;t1=tmp*cost*cost;t2=aq*n+q*cost*cost;t3=ap*n+q*sin2t;ap*n+p=t1+aq*n+q-t2-t3;a

11、q*n+q=tmp-t1+t2+t3;ap*n+q=0.0; /置该非对角元素为0.0，下次循环最大便不是它了aq*n+p=0.0; /同上for(j=0;jn;j+) /以下两for语句求S矩阵非主对角元素if(j!=p)&(j!=q)tmp=ap*n+j;ap*n+j=tmp*cost-aq*n+j*sint;aq*n+j=tmp*sint+aq*n+j*cost;for(i=0;in;i+)if(i!=p)&(i!=q)ai*n+p=ap*n+i;ai*n+q=aq*n+i;for(i=0;in;i+) /该for语句递推求矩阵U fm=uip;uip=fm*cost-uiq*sint;

12、uiq=fm*sint+uiq*cost;return (0);4. 调试及测试 1.0 2.0 3.0矩阵A= 2.0 2.0 5.0 3.0 5.0 1.01.屏幕输出如下：2. 文本文件输出见文件“juzhen.txt”。3. 结果正确性分析利用数学软件Mathematica 8.0计算特征值（即矩阵S主对角线元素）及对应特征向量组U如下截图。（Eigenvalues 为计算矩阵特征值函数，Eigenvectors 为计算矩阵特征向量函数）5. 设计总结拿到这个题目，我首先想到用解方程组的方法来求解矩阵U和S,但是后来发现：假如这是一个三阶矩阵，那么通过A=USUT可以得到含9个方程的方

13、程组，而矩阵U和S中共12个未知数，故通过这种方法建模被否定了。后查阅资料，理解了这是一个矩阵特征分解的问题，故转化为两部分求解：特征值和特征向量。可以通过雅克比变换来求解。从txt文本文件中读入数据又是一大障碍，本以为存在文本文件中的浮点数据只需通过浮点格式字符串就能将其读出，结果运行后屏幕显示一团乱。后查阅资料得知，文本文件中的数据都是字符数据。故先是将其中的字符数据读出，然后用字符与浮点之间的转换函数还原浮点数据。最后碰到的一个问题是主函数与被调用函数之间的衔接问题，就是一个指针所指的实参和形参，要么都是一维数组，要么都是二维数组，否则就会出错。6. 心得体会及致谢起先听说我们这是编程课，感觉很激动，心想上学期所学内容终于有用武之地。但拿到题目后，我像立刻被泼了一盆冷水似的，感觉这个题目无从下手。咬着牙挣扎了两天仍然是毫无头绪，于是乎我想到了放弃，找到老师表达想换题目的想法。起初对老师的拒绝感到很不爽，可老师建议我去查阅资料，先把矩阵论中没学过的只是给搞懂然后再去编程。我先查阅了一些关于算法的书籍，结果其中的程序看不懂，故又借了关于矩阵的书一起研究。经过一个多星期的努力，