雅克比法求矩阵特征值特征向量.docx
《雅克比法求矩阵特征值特征向量.docx》由会员分享，可在线阅读，更多相关《雅克比法求矩阵特征值特征向量.docx（17页珍藏版）》请在冰豆网上搜索。
雅克比法求矩阵特征值特征向量
C语言课程设计报告
课程名称：
计算机综合课程设计
学院：
土木工程学院
设计题目：
矩阵特征值分解
级别：
B
学生姓名：
学号：
同组学生：
无
学号：
无
指导教师：
2012年9月5日
C语言课程设计任务书
（以下要求需写入设计报告书）
学生选题说明:
Ø以所发课程设计要求为准，请同学们仔细阅读；
Ø本任务书提供的设计案例仅供选题参考；也可自选，但难易程度需难度相当；
Ø鼓励结合本专业（土木工程、力学）知识进行选题，编制程序解决专业实际问题。
Ø限2人选的题目可由1-2人完成（A级）；限1人选的题目只能由1人单独完成（B级）；
设计总体要求:
Ø采用模块化程序设计；
Ø鼓励可视化编程；
Ø源程序中应有足够的注释；
Ø学生可自行增加新功能模块（视情况可另外加分）；
Ø必须上机调试通过；
Ø注重算法运用，优化存储效率与运算效率；
Ø需提交源程序（含有注释）及相关文件（数据或数据库文件）；
（cpp文件、txt或dat文件等）
Ø提交设计报告书，具体要求见以下说明。
设计报告格式：
1．课程设计任务书（功能简介、课程设计要求）；
2．系统设计（包括总体结构、模块、功能等，辅以程序设计组成框图、流程图解释）；
3．模块设计（主要模块功能、源代码、注释（如函数功能、入口及出口参数说明，函数调用关系描述等）；
4．调试及测试：
（调试方法，测试结果的分析与讨论，截屏、正确性分析）；
5．设计总结：
（编程中遇到的问题及解决方法）；
6．心得体会及致谢；
参考文献
1.课程设计任务书
功能简介：
a）输入一个对称正方矩阵A，从文本文件读入；
b）对矩阵A进行特征值分解，将分解结果：
即U矩阵、S矩阵输出至文本文件；
c）将最小特征值及对应的特征向量输出至文本文件；
d）验证其分解结果是否正确。
提示：
A=USUT,具体算法可参考相关文献。
功能说明：
矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。
注：
以三阶对称矩阵为例
2.系统设计
3.模块设计
#include
#include
#include
intmain（）
{
FILE*fp;
inttezheng（double*a,intn,double*s,double*u,doubleeps,intitmax）;//函数调用声明
inti,j,p,itmax=1000;//itmax为最大循环次数
doubleeps=1e-7,s[3][3],u[3][3];//eps为元素精度，s为对角矩阵S，u为矩阵U
doublea[9];//a为待分解矩阵A
i=tezheng（a,3,s,u,eps,1000）;
if（i>0）//i对应函数中的返回值it
{
if（（fp=fopen（"juzhen.txt","w"））==NULL）//打开待输入txt文件
{
printf（"无法打开文件.\n"）;
return;
}
printf（"U矩阵为:
\n"）;//下几句分别向屏幕和txt文件输入矩阵U
fprintf（fp,"U矩阵为:
\n"）;
for（i=0;i<3;i++）
{
for（j=0;j<3;j++）
{
printf（"%10.6f",u[i][j]）;
fprintf（fp,"%10.6f",u[i][j]）;
}
printf（"\n"）;
fprintf（fp,"\n"）;
}
printf（"S对角矩阵为:
\n"）;//下几句分别向屏幕和txt文件输入矩阵S
fprintf（fp,"S对角矩阵为:
\n"）;
for（i=0;i<3;i++）
{
for（j=0;j<3;j++）
{
printf（"%10.6f",s[i][j]）;
fprintf（fp,"%10.6f",s[i][j]）;
}
printf（"\n"）;
fprintf（fp,"\n"）;
}
p=0;
for（i=0;i<3;i++）//下面几句为求最小特征值及其对应特征向量，并输出到屏幕和txt文件中
if（s[i][i]p=i;
printf（"最小特征值为:
%10.6f\n",s[p][p]）;
fprintf（fp,"最小特征值为:
%10.6f\n",s[p][p]）;
printf（"对应特征向量为:
\n"）;
fprintf（fp,"对应特征向量为:
\n"）;
for（i=0;i<3;i++）
{
printf（"%10.6f\n",u[i][p]）;
fprintf（fp,"%10.6f\n",u[i][p]）;
}
}
}
inttezheng（double*a,intn,doubles[3][3],doubleu[3][3],doubleeps,intitmax）
{
FILE*A;
charline[1000];//存放从文件juzhenA.txt读入的数据
char*k="";
inti,j,p,q,it;
doublesint,cost,sin2t,cos2t,d,tmp,t1,t2,t3,fm;
it=0;
if（（A=fopen（"juzhenA.txt","r"））==NULL）
{
printf（"无法打开文件\n"）;
return;
}
fgets（line,1000,A）;//从文件指针A指向的文件，即juzhenA.txt中读入字符数据到数组line中
strtok（line,k）;// 原型char*strtok（char*s,constchar*delim）;说明：
strtok（）用来将字符串分割成一个个片段。
参数s指向欲分割的字符串，参数delim则为分割字符串，当strtok（）在参数s的字符串中发现到参数delim的分割字符时则会将该字符改为\0字符。
在第一次调用时，strtok（）必需给予参数s字符串，往后的调用则将参数s设置成NULL。
每次调用成功则返回被分割出片段的指针。
for（i=0;ia[i]=atof（strtok（NULL,k））;//atof函数功能:
把字符串转换成浮点数
fclose（A）;//释放A指针
for（i=0;i{
u[i][i]=1.0;
for（j=0;jif（i!
=j）
u[i][j]=0.0;
}
while（it{
it++;
d=0.0;
for（i=1;ifor（j=0;j{
tmp=fabs（a[i*n+j]）;
if（tmp>d）
{
d=tmp;
p=i;q=j;
}
}
if（d{
for（i=0;ifor（j=0;js[i][j]=a[i*n+j];
return（it）;
}
sint=2*a[p*n+q];//以下几句为递推公式,求S矩阵主对角元素
sin2t=2*a[p*n+q]/（sqrt（4*a[p*n+q]*a[p*n+q]+（a[q*n+q]-a[p*n+p]）*（a[q*n+q]-a[p*n+p]）））;
if（a[q*n+q]-a[p*n+p]<0.0）
sin2t=-sin2t;
cos2t=sqrt（1.0-sin2t*sin2t）;
sint=sin2t/（sqrt（2.0*（1.0+cos2t）））;
cost=sqrt（1.0-sint*sint）;
tmp=a[p*n+p];
t1=tmp*cost*cost;
t2=a[q*n+q]*cost*cost;
t3=a[p*n+q]*sin2t;
a[p*n+p]=t1+a[q*n+q]-t2-t3;
a[q*n+q]=tmp-t1+t2+t3;
a[p*n+q]=0.0;//置该非对角元素为0.0，下次循环最大便不是它了
a[q*n+p]=0.0;//同上
for（j=0;jif（（j!
=p）&&（j!
=q））
{
tmp=a[p*n+j];
a[p*n+j]=tmp*cost-a[q*n+j]*sint;
a[q*n+j]=tmp*sint+a[q*n+j]*cost;
}
for（i=0;iif（（i!
=p）&&（i!
=q））
{
a[i*n+p]=a[p*n+i];
a[i*n+q]=a[q*n+i];
}
for（i=0;i{
fm=u[i][p];
u[i][p]=fm*cost-u[i][q]*sint;
u[i][q]=fm*sint+u[i][q]*cost;
}
}
return（0）;
}
4.调试及测试
1.02.03.0
矩阵A=2.02.05.0
3.05.01.0
1.屏幕输出如下：
2.文本文件输出见文件“juzhen.txt”。
3.结果正确性分析
利用数学软件Mathematica8.0计算特征值（即矩阵S主对角线元素）及对应特征向量组U如下截图。
（Eigenvalues[]为计算矩阵特征值函数，Eigenvectors[]为计算矩阵特征向量函数）
5.设计总结
拿到这个题目，我首先想到用解方程组的方法来求解矩阵U和S,但是后来发现：
假如这是一个三阶矩阵，那么通过A=USUT可以得到含9个方程的方程组，而矩阵U和S中共12个未知数，故通过这种方法建模被否定了。
后查阅资料，理解了这是一个矩阵特征分解的问题，故转化为两部分求解：
特征值和特征向量。
可以通过雅克比变换来求解。
从txt文本文件中读入数据又是一大障碍，本以为存在文本文件中的浮点数据只需通过浮点格式字符串就能将其读出，结果运行后屏幕显示一团乱。
后查阅资料得知，文本文件中的数据都是字符数据。
故先是将其中的字符数据读出，然后用字符与浮点之间的转换函数还原浮点数据。
最后碰到的一个问题是主函数与被调用函数之间的衔接问题，就是一个指针所指的实参和形参，要么都是一维数组，要么都是二维数组，否则就会出错。
6.心得体会及致谢
起先听说我们这是编程课，感觉很激动，心想上学期所学内容终于有用武之地。
但拿到题目后，我像立刻被泼了一盆冷水似的，感觉这个题目无从下手。
咬着牙挣扎了两天仍然是毫无头绪，于是乎我想到了放弃，找到老师表达想换题目的想法。
起初对老师的拒绝感到很不爽，可老师建议我去查阅资料，先把矩阵论中没学过的只是给搞懂然后再去编程。
我先查阅了一些关于算法的书籍，结果其中的程序看不懂，故又借了关于矩阵的书一起研究。
经过一个多星期的努力，终于懂了，编出来了。
在此，谢谢老师的建议。
我想这对我以后自学任何东西都是有帮助的。
参考文献：
1、《基于MATLAB的实用数值计算》石辛民郝整清编著清华大学出版社·北京交通大学出版社
2006年2月第1版
2、《C语言算法速查手册》程晓旭耿鲁静张海王勇编著人民邮电出版社2009年10月第1版
3、《Mathematica4教程》丁大正编著电子工业出版社2002年3月第1版
4、XX百科
附：
所用原理
所用函数：
欢迎您的下载，
资料仅供参考！
致力为企业和个人提供合同协议，策划案计划书，学习资料等等
打造全网一站式需求
升级会员 