雅克比过关法求特征值和特征向量
《实验8：Jacobi法求实对称矩阵的特征值及特征向量》
实验名称实验8实验地点6A-XXX 实验类型设计实验学时 2 实验日期20 /X/X ★撰写注意：版面格式已设置好（不得更改），填入内容即可。
一、实验目的
1.Jacobi法求实对称矩阵的特征值及特征向量
二、实验内容
1.实验任务
1.Jacobi法求实对称矩阵的特征值及特征向量
2.程序设计
1）数据输入（输入哪些数据、个数、类型、来源、输入方式）
double a[N][N], int n
2）数据存储（输入数据在内存中的存储）
函数
void Jacobi(double a[N][N], int n)
3）数据处理（说明处理步骤。
若不是非常简单，需要绘制流程图）
1.输入要处理的数据进入变量中
2.进行函数处理
3.输出函数处理结果
4）数据输出（贴图：程序运行结果截图。
图幅大小适当，不能太大）
三、实验环境
1.操作系统：WINDOWS 7及以上
2.开发工具：VS 2015
3.实验设备：PC。
雅克比法求矩阵特征值特征向量
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.系统设计总体结构main函数tezheng函数从文本文找矩阵A递推求矩递推求矩向屏幕和向屏幕和求最小特阵S件中读入阵U中非对角txt文件输txt文件输征值及其数组A入矩阵U入矩阵S元素中的对应特征最大值，向量，并并记下其输出到屏位置幕和txt文件中3.模块设计#include<stdio.h>#include<stdlib.h>#include<math.h>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为待分解矩阵Ai=tezheng(a,3,s,u,eps,1000);WORD格式if(i>0)//i对应函数中的返回值it{if((fp=fopen("juzhen.txt","w"))==NULL)//打开待输入txt文件{printf("无法打开文件.\n");return;}printf("U矩阵为:\n");//下几句分别向屏幕和txt文件输入矩阵Ufprintf(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文件输入矩阵Sfprintf(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]<s[0][0])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()用来将字符串分割成一个个片段。
雅可比矩阵求特征方程
雅可比矩阵求特征方程雅可比矩阵是一种将多元函数的偏导数矩阵以及变量向量形式组合起来的矩阵。
在数学中，雅可比矩阵的特征方程是对于该矩阵进行特征值分解之后得到的特征向量满足的特殊方程。
本文将详细介绍雅可比矩阵的概念、特征值与特征向量的计算方法，以及特征方程的推导过程。
为了更好地理解雅可比矩阵，我们首先给出它的定义。
设有多元函数f(x1, x2, ..., xn)，其中x1, x2, ..., xn为n个变量。
那么该函数关于这n个变量的雅可比矩阵J可以表示为：J = ( ∂f1/∂x1, ∂f1/∂x2, ..., ∂f1/∂xn;∂f2/∂x1, ∂f2/∂x2, ..., ∂f2/∂xn;...∂fn/∂x1, ∂fn/∂x2, ..., ∂fn/∂xn )其中，∂fj/∂xi表示f对变量xi的偏导数。
这个矩阵的维度为n×n，每个元素都是一个偏导数。
现在我们考虑如何求解雅可比矩阵的特征方程。
设J是一个n阶方阵，特征值为λ，特征向量为v，那么有以下的特征方程：Jv=λv对于特征向量v的每一个分量vi，我们可以写作：Jv = (j1, j2, ..., jn)=λv= (λv1, λv2, ..., λvn)根据矩阵与向量的乘法规则，有：ðf1/∂x1 * v1 + ðf1/∂x2 * v2 + ... + ðf1/∂xn * vn = λv1ðf2/∂x1 * v1 + ðf2/∂x2 * v2 + ... + ðf2/∂xn * vn = λv2 ...ðfn/∂x1 * v1 + ðfn/∂x2 * v2 + ... + ðfn/∂xn * vn = λvn 将上述方程用矩阵的形式表示，可以写为：Jv=λv⇒ ( ∂f1/∂x1, ∂f1/∂x2, ..., ∂f1/∂xn;∂f2/∂x1, ∂f2/∂x2, ..., ∂f2/∂xn;...∂fn/∂x1, ∂fn/∂x2, ..., ∂fn/∂xn )* (v1, v2, ..., vn)= (λv1, λv2, ..., λvn)那么可以得到以下的方程组：∂f1/∂x1 * v1 + ∂f1/∂x2 * v2 + ... + ∂f1/∂xn * vn = λv1∂f2/∂x1 * v1 + ∂f2/∂x2 * v2 + ... + ∂f2/∂xn * vn = λv2 ...∂fn/∂x1 * v1 + ∂fn/∂x2 * v2 + ... + ∂fn/∂xn * vn = λvn以上方程可化简为：(∂f1/∂x1 - λ)v1 + ∂f1/∂x2 * v2 + ... + ∂f1/∂xn * vn = 0∂f2/∂x1 * v1 + (∂f2/∂x2 - λ)v2 + ... + ∂f2/∂xn * vn = 0...∂fn/∂x1 * v1 + ∂fn/∂x2 * v2 + ... + (∂fn/∂xn - λ)vn = 0注意到方程左侧的形式与行列式的形式相似，我们可以进一步将方程化简为：(∂f1/∂x1 - λ)v1 + ∂f1/∂x2 * v2 + ... + ∂f1/∂xn * vn = 0∂f2/∂x1 * v1 + (∂f2/∂x2 - λ)v2 + ... + ∂f2/∂xn * vn = 0...∂fn/∂x1 * v1 + ∂fn/∂x2 * v2 + ... + (∂fn/∂xn - λ)vn = 0写成矩阵的形式：(∂f1/∂x1 - λ, ∂f1/∂x2, ..., ∂f1/∂xn;∂f2/∂x1, (∂f2/∂x2 - λ), ..., ∂f2/∂xn;...∂fn/∂x1, ∂fn/∂x2, ..., (∂fn/∂xn - λ) )* (v1, v2, ..., vn)=(0,0, 0这是一个关于λ和v的齐次线性方程组，若存在非零解v，则其中必然存在一个非零特征值λ。
雅可比算法求矩阵的特征值和特征向量
雅可⽐算法求矩阵的特征值和特征向量⽬的求⼀个实对称矩阵的所有特征值和特征向量。
前置知识对于⼀个实对称矩阵A ，必存在对⾓阵D 和正交阵U 满⾜D =U T AUD 的对⾓线元素为A 的特征值，U 的列向量为A 的特征向量。
定义n 阶旋转矩阵G (p ,q ,θ)=1⋯0 ⋱ 1 cos θ−sin θ 10 ⋱ 0即在单位矩阵的基础上，修改a pp =a qq =cos θ,a qp =−a pq =sin θ对于n 阶向量α，α⋅G (p ,q ,θ)的⼏何意义是把α在与第p 维坐标轴和第q 维坐标轴平⾏的平⾯内旋转⾓度θ，并且旋转后的模长保持不变。
算法原理⼤概思路使通过旋转变换使⾮对⾓线上的元素不断变⼩，最后得到与原矩阵相似的对⾓矩阵。
每次找到矩阵A 绝对值最⼤的的⾮对⾓线元素，设为a pq ，令U =G (p ,q ,θ)，将A 变换为U T AU变换后的值为通过令b p ,q =0解得θ=12arctan 2a pq a qq−a pp 特别地当a qq =a pp 时θ=π4注意到旋转操作并不会改变每个⾏向量或列向量的模长，即矩阵A 的F-范数||A ||F =∑i ∑j a 2ij 是不变的，并且通过计算可以得出$$b_{ip}2+b_{iq}2=a_{ip}2+a_{iq}2$$从⽽可以得知⾮对⾓线元素的平⽅和变⼩，对⾓线上元素的平⽅和增⼤，故⾮主对⾓线上元素的平⽅和收敛。
算法流程（1）令矩阵T =E ，即初始化单位矩阵（2）找到A 中绝对值最⼤的⾮对⾓选元素a pq（3）找到对应的⾓度θ，构造矩阵U =G (p ,q ,θ)（4）令A =U T AU ,T =TU（5）不停地重复（2）到（4），直到a pq <ϵ或迭代次数超过某个限定值，则A 的对⾓线元素近似等于A 的特征值，T 的列向量为A 的特征向量代码#include<bits/stdc++.h>using namespace std;const int N=1005;const double eps=1e-5;const int lim=100;int n,id[N];[√double key[N],mat[N][N],EigVal[N],EigVec[N][N],tmpEigVec[N][N];bool cmpEigVal(int x,int y){return key[x]>key[y];}void Find_Eigen(int n,double (*a)[N],double *EigVal,double (*EigVec)[N]){for (int i=1;i<=n;i++)for (int j=1;j<=n;j++)EigVec[i][j]=0;for (int i=1;i<=n;i++) EigVec[i][i]=1.0;int count=0;while (1){//统计迭代次数count++;//找绝对值最⼤的元素double mx_val=0;int row_id,col_id;for (int i=1;i<n;i++)for (int j=i+1;j<=n;j++)if (fabs(a[i][j])>mx_val) mx_val=fabs(a[i][j]),row_id=i,col_id=j;if (mx_val<eps||count>lim) break;//进⾏旋转变换int p=row_id,q=col_id;double Apq=a[p][q],App=a[p][p],Aqq=a[q][q];double theta=0.5*atan2(-2.0*Apq,Aqq-App);double sint=sin(theta),cost=cos(theta);double sin2t=sin(2.0*theta),cos2t=cos(2.0*theta);a[p][p]=App*cost*cost+Aqq*sint*sint+2.0*Apq*cost*sint;a[q][q]=App*sint*sint+Aqq*cost*cost-2.0*Apq*cost*sint;a[p][q]=a[q][p]=0.5*(Aqq-App)*sin2t+Apq*cos2t;for (int i=1;i<=n;i++)if (i!=p&&i!=q){double u=a[p][i],v=a[q][i];a[p][i]=u*cost+v*sint;a[q][i]=v*cost-u*sint;u=a[i][p],v=a[i][q];a[i][p]=u*cost+v*sint;a[i][q]=v*cost-u*sint;}//计算特征向量for (int i=1;i<=n;i++){double u=EigVec[i][p],v=EigVec[i][q];EigVec[i][p]=u*cost+v*sint;EigVec[i][q]=v*cost-u*sint;}}//对特征值排序for (int i=1;i<=n;i++) id[i]=i,key[i]=a[i][i];std::sort(id+1,id+n+1,cmpEigVal);for (int i=1;i<=n;i++){EigVal[i]=a[id[i]][id[i]];for (int j=1;j<=n;j++)tmpEigVec[j][i]=EigVec[j][id[i]];}for (int i=1;i<=n;i++)for (int j=1;j<=n;j++)EigVec[i][j]=tmpEigVec[i][j];//特征向量为列向量}int main(){scanf("%d",&n);for (int i=1;i<=n;i++)for (int j=1;j<=n;j++)scanf("%lf",&mat[i][j]);Find_Eigen(n,mat,EigVal,EigVec);printf("EigenValues = ");for (int i=1;i<=n;i++) printf("%lf ",EigVal[i]);printf("\nEigenVector =\n");for (int i=1;i<=n;i++)for (int j=1;j<=n;j++)printf("%lf%c",EigVec[i][j],j==n?'\n':' ');return 0;}Processing math: 100%。
用共轭梯度法解方程,用Jacobi方法求矩阵的全部特征值和特征向量
end
A(1,1)=4;
A(n,n)=4;
A(n,n-1)=1;
X=zeros(n,1);
fori=1:n
X(i,1)=1;
end
%用共轭梯度法求解方程
fprintf('方程的精确解\n');
X
fprintf('用共轭法求解方程\n');
x=cg(A,b)
%用方法求解方程的特征值和特征向量
-0.1859 -0.0000 0.0000 0.0000 5.6180 -0.0000 -0.0000
-0.2236 0.0000 -0.0000 0.0000 0.0000 5.4142 0.0000
-0.2558 0.0000 -0.0000 0.0000 0.0000 0.0000 5.1756
0.1859 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000
0.1436 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000
-0.0977 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
0.6634
1.1794
0.6188
1.3453
用方法Jacobi求解矩阵的全部特征值及特征向量
D =
Columns 1 through 7
4.0000 0 0 0 0 0 0
0.1382 5.9021 0.0000 0.0000 -0.0000 0.0000 0.0000
-0.2629 0.0000 5.6180 0.0000 0.0000 -0.0000 -0.0000
雅克比法解特征值
A1 仍然是实对称阵，且 A1与 A 的特征值相同。
把式(3.12)的右端乘开，并与左端比较，得到A1 的元素计算公式：
A
(1)
a a (p11) a (1) q1 a (1) n1
(1) 11
a a a a
(1) 1p
4
sgn(a pq ).
2 tan 1 2 1 tan c
4
， 故 tan 1,
故可取 tan
sgn(c ) c c2 1
t
t sin t cos 1 t 2 则 1 1 cos 1 tan 2 1 t 2
i2 ( B) B .
2 i 1
F
n
对于 n 维列向量 x， Upq x 相当于把坐标轴Oxp和 Oxq于所在平面内旋转角度 。 记矩阵 A = [aij]n×n ，对A作一次正交相似变换， 得到矩阵A1。
记
A1= UpqT A Upq = [aij(1)]n×n
(3.12)
故可选择角θ，使
c12 c21 (a22 a11 ) sin cos a12 (cos sin ) 0
2 2
c12 c21 (a22 a11 ) sin cos a12 (cos2 sin 2 ) 0
即 tan 2 2a12 ， a11 a22
1 2 a aij n(n 1) i j
2 pq
从而
a
i j
(1) 2 ij
2 2 1 aij n(n 1) i j
jacobi旋转法,Jacobi过关法
实验名称：项目一 Jacobi 旋转法，Jacobi 过关法实验题目：用Jacobi 旋转法，Jacobi 过关法计算下面矩阵A 的全部特征值以及特征向量2100121001210012A -⎡⎤⎢⎥--⎢⎥=⎢⎥--⎢⎥-⎣⎦ 实验过程：1、 Jacobi 旋转法程序如下：function [EigV alMat,EigV ecMat]=JocobiRot(A,eps)% 本程序是根据jacobi 旋转法求实对称矩阵的全部特征值和特征向量% 输入变量:A 为对称实矩阵,eps 为允许误差.% 输出变量:EigValMat 为A 的特征值,EigVecMat 为特征向量n=size(A);n=n(1); % n 为矩阵A 的阶数P=eye(n); % P 为旋转矩阵,赋初值为单位阵trc=1; % trc 为矩阵A 的非对角元素的平方和,赋初值为1; num=0; % num 设置为累加器,记录迭代的次数while abs(trc)>=eps % 进行正交变换A=PAP'将A 的两个绝对值最大的非对角元素化零,直到所有非对角元素的平方和小于给定的eps,则结束循环.MaxMes=FinMax(A); % 寻找绝对值最大的非对角元素的位置及所有非对角元素的平方和l=MaxMes(1); % l 为绝对值最大的非对角元素的行号s=MaxMes(2); % s 为绝对值最大的非对角元素的列号trc=MaxMes(3); % trc 为矩阵A 的非对角元素的平方和RotMat=ComRotMat(A,l,s); % 计算此次旋转的平面旋转矩阵RotMatA=RotMat*A*RotMat'; % 对当前A 进行一次旋转变换将A 的两个绝对值最大的非对角元素化零,并仍记为AP=RotMat*P; % 记录到目前为止所有旋转矩阵的乘积num=num+1; % 记录已经进行旋转的次数endEigValMat=A;EigVecMat=P;numfunction MaxMes=FinMax(A)% 对上三角进行扫描，寻找矩阵A的绝对值最大的非对角元素的位置及所有非对角元素的平方和% 输入量:实对称矩阵A% 输出量:MaxMes记录矩阵A的绝对值最大的非对角元素的行号列号及所有非对角元素的平方和n=size(A);n=n(1);trc=0;MaxVal=0; % MaxVal记录矩阵A的绝对值最大的非对角元素赋初值为0for i=1:n-1for j=i+1:ntrc=trc+2*A(i,j)^2;if abs(A(i,j))>MaxValMaxVal=abs(A(i,j));l=i; % l为绝对值最大的非对角元素的行号s=j; % s为绝对值最大的非对角元素的列号endendendMaxMes=[l,s,trc];function RotMat=ComRotMat(A,l,s)% 计算平面旋转矩阵% 输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s% 输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);if A(l,l)==A(s,s)c1=1/sqrt(2); s1=1/sqrt(2);elsetan2=2*A(l,s)/(A(l,l)-A(s,s));c2=1/sqrt(1+(tan2)^2);s2=tan2*c2;c1=sqrt((1+c2)/2);s1=s2/(2*c1);endRotMat=eye(n);RotMat(l,l)=c1;RotMat(s,s)=c1;RotMat(l,s)=s1;RotMat(s,l)=-s1;2、Jacobi过关法程序如下：function [EigV alMat,EigV ecMat]=JocobiGG(A,eps)% 本程序是根据jacobi过关法求实对称矩阵的全部特征值和特征向量% 输入变量:A为对称实矩阵,eps为允许误差.% 输出变量:EigValMat为A的特征值,EigVecMat为特征向量n=size(A);n=n(1); % n为矩阵A的阶数P=eye(n); % P为旋转矩阵,赋初值为单位阵trc=0; % trc为矩阵A的非对角元素的平方和,赋初值为0;num=0; % num设置为累加器,记录迭代的次数for i=1:n-1for j=i+1:ntrc=trc+2*A(i,j)^2; % 计算矩阵A的非对角元素平方和endendr0=trc^(1/2);r1=r0/n;while r1>=epsfor i=1:n-1for j=i+1:nif abs(A(i,j))>=r1 % 对abs(A(i,j))>=r1的元素进行旋转变换，将A(i,j)化为0，其余元素视为“过关”，不作相应的变换l=i;s=j;RotMat=ComRotMat(A,l,s); % 计算此次旋转的平面旋转矩阵RotMatA=RotMat*A*RotMat';P=RotMat*P; % 记录到目前为止所有旋转矩阵的乘积num=num+1; % 记录已经进行旋转的次数endendendr1=r1/n; % 当所有非对角元素绝对值都小于r1后，缩小阀值endEigValMat=A;EigVecMat=P;numfunction RotMat=ComRotMat(A,l,s)% 计算平面旋转矩阵% 输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s% 输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);if A(l,l)==A(s,s)c1=1/sqrt(2); s1=1/sqrt(2);elsetan2=2*A(l,s)/(A(l,l)-A(s,s));c2=1/sqrt(1+(tan2)^2);s2=tan2*c2;c1=sqrt((1+c2)/2);s1=s2/(2*c1);endRotMat=eye(n);RotMat(l,l)=c1;RotMat(s,s)=c1;RotMat(l,s)=s1;RotMat(s,l)=-s1;实验结果：1、Jacobi旋转法运行结果如下：>>A=[2 -1 0 0;-1 2 -1 0;0 -1 2 -1;0 0 -1 2];eps=0.0001;[EigValMat,EigV ecMat]=JocobiRot(A,eps) num =7EigValMat =0.3820 0 0 -0.00000 3.6180 0 00 -0.0000 1.3820 00 0 0 2.6180EigVecMat =0.3717 0.6015 0.6015 0.3717-0.3717 0.6015 -0.6015 0.3717-0.6015 -0.3717 0.3717 0.60150.6015 -0.3717 -0.3717 0.60152、Jacobi过关法运行结果如下：>> [EigValMat,EigVecMat]=JocobiGG(A,eps)num =16EigValMat =0.3820 0.0000 -0.0000 -0.00000.0000 3.6180 -0.0000 0.0000-0.0000 -0.0000 1.3820 0.00000.0000 0.0000 0.0000 2.6180EigVecMat =0.3717 0.6015 0.6015 0.3717-0.3718 0.6015 -0.6015 0.3717-0.6015 -0.3717 0.3717 0.60150.6015 -0.3717 -0.3718 0.6015实验名称：项目二 Romberg 方法求数值积分 实验题目：用Romberg 方法计算下列积分：（1）30sin x e xdx ⎰，要求准确到610-（2）210x e dx -⎰，要求准确到610-实验过程：function I=romberg(fun,a,b,eps)m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while 1N=2^(m-1);h=h/2;I=I/2;for i=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;while M>1T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;endif abs(T(k,k)-T(k,k-1))<epsbreak;endm=m+1;endI=T(k,k);实验结果：>> fun=inline('exp(x)*sin(x)');a=0;b=3;eps=0.000001;I=romberg(fun,a,b,eps)I =11.8595>> fun=inline('exp(-x^2)');a=0;b=1;eps=0.000001;I=romberg(fun,a,b,eps)I =0.7468。
雅克比（Jacobi）方法
雅克⽐（Jacobi）⽅法可以⽤来求解协⽅差矩阵的特征值和特征向量。
雅可⽐⽅法(Jacobian method)求全积分的⼀种⽅法，把拉格朗阶查⽪特⽅法推⼴到求n个⾃变量⼀阶⾮线性⽅程的全积分的⽅法称为雅可⽐⽅法。
雅克⽐迭代法的计算公式简单，每迭代⼀次只需计算⼀次矩阵和向量的乘法，且计算过程中原始矩阵A始终不变，⽐较容易并⾏计算。
考虑线性⽅程组Ax=b时，⼀般当A为低阶稠密矩阵时，⽤主元消去法解此⽅程组是有效⽅法。
但是，对于由⼯程技术中产⽣的⼤型稀疏矩阵⽅程组（A的阶数很⾼，但零元素较多，例如求某些偏微分⽅程数值解所产⽣的线性⽅程组），利⽤迭代法求解此⽅程组就是合适的，在计算机内存和运算两⽅⾯，迭代法通常都可利⽤A中有⼤量零元素的特点。
雅克⽐迭代法就是众多迭代法中⽐较早且较简单的⼀种，其命名也是为纪念普鲁⼠著名数学家雅可⽐。
原理【收敛性】设Ax=b，其中A=D+L+U为⾮奇异矩阵，且对⾓阵D也⾮奇异，则当迭代矩阵J的谱半径ρ（J）<1时，雅克⽐迭代法收敛。
⾸先将⽅程组中的系数矩阵A分解成三部分，即：A = L+D+U，其中D为对⾓阵，L为下三⾓矩阵，U为上三⾓矩阵。
之后确定迭代格式，X^(k+1) = B*X^(k) +f ，（这⾥^表⽰的是上标，括号内数字即迭代次数），其中B称为迭代矩阵，雅克⽐迭代法中⼀般记为J。
（k = 0,1，......）再选取初始迭代向量X^(0)，开始逐次迭代。
【优缺点】雅克⽐迭代法的优点明显，计算公式简单，每迭代⼀次只需计算⼀次矩阵和向量的乘法，且计算过程中原始矩阵A始终不变，⽐较容易并⾏计算。
然⽽这种迭代⽅式收敛速度较慢，⽽且占据的存储空间较⼤，所以⼯程中⼀般不直接⽤雅克⽐迭代法，⽽⽤其改进⽅法。
实现通过雅克⽐(Jacobi)⽅法求实对称矩阵的特征值和特征向量操作步骤：S′=G T SG,其中G是旋转矩阵，S′和S均为实对称矩阵，S′和S有相同的Frobenius norm，可以⽤⼀个最简单的3维实对称矩阵为例，根据公式进⾏详细推导(参考 )：通过旋转矩阵将对称矩阵转换为近似对⾓矩阵，进⽽求出特征值和特征向量，对⾓矩阵中主对⾓元素即为S近似的实特征值。
jacobi方法求特征值和特征向量 例题
一、引言Jacobi方法是一种用于计算矩阵特征值和特征向量的迭代数值方法。
它是数值线性代数中的重要算法之一，广泛应用于科学计算、工程技术和金融领域。
本文将通过一个例题来介绍Jacobi方法的原理和求解过程，并分析其在实际问题中的应用。
二、Jacobi方法的原理Jacobi方法是一种通过迭代对矩阵进行相似变换，使得原矩阵逐步转化为对角矩阵的方法。
通过数值迭代，可以逐步逼近矩阵的特征值和对应的特征向量。
其基本原理如下：1. 对称矩阵特征值问题：对于对称矩阵A，存在一个正交矩阵P，使得P^T * A * P = D，其中D为对角矩阵，其对角线上的元素为A的特征值。
所以我们可以通过迭代找到P，使得P逼近正交矩阵，从而逼近A的特征值和特征向量。
2. Jacobi迭代：Jacobi方法的基本思想是通过正交相似变换，逐步将矩阵对角化。
具体来说，对于矩阵A，找到一个旋转矩阵G，使得A' = G^T * A * G为对角矩阵，然后递归地对A'进行相似变换，直到达到精度要求。
三、Jacobi方法求解特征值和特征向量的例题考虑以下矩阵A：A = [[4, -2, 2],[-2, 5, -1],[2, -1, 3]]我们将通过Jacobi方法来计算矩阵A的特征值和特征向量。
1. 对称化矩阵我们需要对矩阵A进行对称化处理。
对称化的思路是找到正交矩阵P，使得P^T * A * P = D，其中D为对角矩阵。
我们可以通过迭代找到逼近P的矩阵序列，直到达到一定的精度。
2. Jacobi迭代在Jacobi迭代的过程中，我们需要找到一个旋转矩阵G，使得A' =G^T * A * G为对角矩阵。
具体的迭代过程是：找到矩阵A中绝对值最大的非对角元素a[i][j]，然后构造一个旋转矩阵G，将a[i][j]置零。
通过迭代地对A'进行相似变换，最终使得A'的非对角元素逼近零，即达到对角化的目的。
3. 计算特征值和特征向量经过一定次数的Jacobi迭代后，得到了对称矩阵A的对角化矩阵D和正交矩阵P。
jacobi方法求取实对称矩阵的特征值和特征向量
#ifnde f ksf loat#defin e ksf loat doub le#e ndif#ifn def i nt16#defin e int16 i nt#e ndif#ifn def u int16#defi ne ui nt16unsig ned i nt#e ndif#ifn def k snew#defi ne ks new(s,t) ((t*)m alloc((s)*sizeo f(t)))#en dif#ifnd ef ks free#defin e ksf ree(p) if(p!=0){free(p);p=0;}#endi f/*ksJac obiGe tReal SymMa trixF eatur eValu e用jac obi方法求取*实对称*矩阵的特征值和特征向量 aMa trix[in]: n 阶实对称阵 sid e[in]: 矩阵阶大小 eps[in]:控制精度 ma xTime s[in]: 最大迭代次数返回值：前sid e*sid e个单元存储特征向量，一列为一特征向量；后side个存储单元为特征值，也就是总大小为 siz e = s ide*(side+1) 存储单元*/ksfl oat *ksJac obiGe tReal SymMa trixF eatur eValu e( ks float aMat rix[] , ui nt16side,ksf loateps , uint16 ma xTime s ){ ksfl oat *t, *s, *s1, *ve ctor, *q1; floa t a1, c , s2 , t1 , maxE ps =0 , d1 , d2, r= 1;u int16 i ,j , p , q, m , time s = 0 , le n = 0; uin t16 s tartI ndex= 0 , endI ndex; endI ndex= sid e + s tartI ndex;side+= st artIn dex;len = side*side; t = ksne w(len , ks float); s = ksn ew(le n , k sfloa t); s1 = ks new(l en ,ksflo at);q1 = k snew(len , ksfl oat);vecto r = k snew(len + side , ks float); mem set ( vect or ,0 , (len+s ide)*sizeo f(ksf loat)); for( i = star tInde x ; i < en dInde x ; i++) { ve ctor[i*sid e+i]= 1;//将初始Q[i][j]置为单位矩阵 }w hile(r >=eps && tim es <maxTi mes){p = q= sta rtInd ex; m axEps = 0;f or(i= sta rtInd ex; i < en dInde x; i++) {for(j = st artIn dex;j < e ndInd ex; j++) {i f( (j != i) &&fabs(aMatr ix[i*side+j]) > maxE ps) { ma xEps= (fl oat)f abs(a Matri x[i*s ide+j]);//找非对角元素绝对值最大的元p = i; q =j;} } }r = ma xEps;//重置当前误差 a1= (fl oat)((aMat rix[q*side+q]-a Matri x[p*s ide+p])/(2*aMat rix[p*side+q]));if(a1 >= 0){t1 = (float)(1/(fabs(a1)+s qrt(1+a1*a1)));}e lse { t1 = (flo at)(-1/(fa bs(a1)+sqr t(1+a1*a1))); } c = (flo at)(1/sqrt(1+t1*t1));s2 =t1*c;memse t ( s , 0, len*size of(ks float)); for( i = star tInde x ; i < en dInde x ; i++) { s[i*side+i] =1;//将初始Q[i][j]置为单位矩阵}s[p*s ide+p] = c ;s[p*sid e+q]=s2; s[q*side+p] = -s2;s[q*s ide+q]=c; fo r(i = star tInde x; i< end Index; i++){f or(j= sta rtInd ex; j < en dInde x; j++){ s1[i*si de+j] = s[j*sid e+i];//将矩阵s1[i][j]化为s[i][j]的转置 } }f or(i= sta rtInd ex; i < en dInde x; i++) {for(j = st artIn dex;j < e ndInd ex; j++) {d1 = 0; fo r(m = star tInde x; m< end Index; m++) {d1 += (flo at)(s1[i*s ide+m]*aMa trix[m*sid e+j]);//计算s1[i][j]*a Matri x[i][j] } t[i*side+j] = d1; } } fo r(i = star tInde x; i< end Index; i++){f or(j= sta rtInd ex; j < en dInde x; j++){ d1 = d2 = 0; for(m =start Index; m < endI ndex; m++) {d1 +=(floa t)(t[i*sid e+m]*s[m*s ide+j]);//计算t[i][j]*s[i][j] d2 += (float)(vec tor[i*side+m]*s[m*si de+j]);//计算Q[i][j]*s[i][j] } aMat rix[i*side+j] = d1; q1[i*side+j] = d2; } } me mcpy( vec tor , q1 , len*sizeo f(ksf loat)); time s +=1; }// 对特征值与特征向量进行整合for(i = st artIn dex;i < e ndInd ex; i++) { ve ctor[len+i] = a Matri x[i*s ide+i]; }k sfree ( t); ksf ree ( s );ksfre e ( s1 );k sfree ( q1 );i f(tim es >maxTi mes){k sfree ( ve ctor); retu rn NU LL; }retur n vec tor;}。
4-2Jocabi方法
a 1 a
2
, d sin bc
(4)计算 A ( k 1) 2 (k ) 2 (k ) (k ) aii c aii d a jj 2cdaij
( k 1)
a
( k 1) jj
d a
2
(k ) ii
c a
2
(k ) jj
2cda
(k ) ij
( k 1) ( k 1) (k ) (k ) ail ali cail da jl a ( k 1) a ( k 1) da ( k ) ca ( k ) lj il jl jl
( k 1) 2 lm
(6)若E ( A
( k 1)
) , 则a
(0) (1)
( k 1) 11
,a
(k )
( k 1) 22
, , a
( k 1) nn
为
特征值，Q V V V 的各列为相应 的特征向量。否则，k+1 k,返回(2)。
一般情况下, Jacobi法不能在有限步内将A 化为对角阵,但有以下收敛性结果:
a
k l ( k ,l ) ( i , j )
2 2 2 akl 2 a a ik jk k i , j
2 E ( A) 2aij E ( A)
4.2.2 古典Jacobi方法
(1) 令k =0,A A. (k ) (k ) (2) 求整数i, j，使得 aij max alm 。
k 1
E A
2 因为1 1, 所以, n(n 1) lim E A
n
雅克比法求矩阵特征值特征向量
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=USU T,具体算法可参考相关文献。
功能说明：矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。
注：以三阶对称矩阵为例2.系统设计#include<>#include<>int main(){FILE *fp;int tezheng(double *a,int n,double *s,double *u,double eps,int itmax); n");return;}printf("U矩阵为:\n"); 幕输出如下：2.文本文件输出见文件“”。
jacobi迭代法求复矩阵特征值和特征向量
jacobi迭代法求复矩阵特征值和特征向量Jacobi迭代法是一种经典的求解复矩阵特征值和特征向量的方法。
在数值分析领域，特征值和特征向量的求解是一个十分重要且常见的问题。
它不仅在理论上有重要意义，还在实际应用中有着广泛的应用，比如在物理、工程、金融等领域。
Jacobi迭代法的提出，极大地简化了这个复杂问题的求解过程，为研究人员和工程师提供了一个高效、可靠的数值计算工具。
我们需要了解什么是特征值和特征向量。
对于一个n阶方阵A，如果存在数λ和一个非零向量x，使得Ax=λx成立，则称λ是A的特征值，x是对应于特征值λ的特征向量。
特征值和特征向量的求解十分重要，因为它们包含了矩阵A的重要特性和信息，对于矩阵的对角化、矩阵的稳定性、矩阵的特征分解等问题有着重要的作用。
接下来，让我们来介绍Jacobi迭代法的基本思想和步骤。
Jacobi迭代法的核心思想是通过一系列相似变换，将原始矩阵对角化，从而得到其特征值和特征向量。
具体步骤如下：1. 我们选择一个n阶方阵A，将其初始化为对角矩阵D，将初始的特征向量矩阵初始化为单位矩阵I。
2. 我们选择两个不同的下标i和j(1≤i，j≤n，i≠j)，使得矩阵A的元素aij为非零元素，即aij≠0。
这两个下标表示我们要进行的相似变换的维度。
3. 我们构造一个旋转矩阵P，使得通过P的相似变换，可以将aij对应的元素变为0。
这一步是Jacobi迭代法的核心步骤，旋转矩阵P的构造涉及到对称双射矩阵的变换和特征值的迭代计算。
4. 我们通过P的相似变换，更新矩阵A和特征向量矩阵I，得到新的对角矩阵D和新的特征向量矩阵。
5. 我们检查新得到的对角矩阵D的非对角线元素是否足够小，如果满足要求，则停止迭代，否则继续进行第2步的操作。
通过这样一系列的迭代操作，我们可以逐步地将矩阵A对角化，并得到其特征值和特征向量。
Jacobi迭代法以其简洁、直观的特点，在复矩阵特征值和特征向量的求解中得到了广泛的应用。
雅克比法求矩阵特征值特征向量
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=USU T,具体算法可参考相关文献。
功能说明：矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。
注：以三阶对称矩阵为例2.系统设计3.模块设计#include<stdio.h>#include<stdlib.h>#include<math.h>int 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,s[3][3],u[3][3]; //eps为元素精度，s为对角矩阵S，u为矩阵Udouble a[9];//a为待分解矩阵Ai=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文件输入矩阵Ufprintf(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文件输入矩阵Sfprintf(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]<s[0][0])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]);}}}int tezheng(double *a,int n,double s[3][3],double u[3][3],double eps,int itmax){FILE *A;char line[1000]; //存放从文件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); //从文件指针A指向的文件，即juzhen A.txt中读入字符数据到数组line 中strtok(line,k); //原型char *strtok(char *s, const char *delim);说明：strtok()用来将字符串分割成一个个片段。
jacobi迭代法求复矩阵特征值和特征向量
题目：深入探究jacobi迭代法求复矩阵特征值和特征向量上线性代数的学习过程中，我们经常会遇到求解复矩阵的特征值和特征向量的问题。
而jacobi迭代法则是一种被广泛应用的方法之一。
本文将深入探讨jacobi迭代法的原理、应用以及个人观点和理解。
### 1. jacobi迭代法的原理和概念jacobi迭代法是一种通过不断相似变换将矩阵对角化的方法，它可以被用于求解实对称矩阵的特征值和特征向量，而在这篇文章中，我们将着重讨论其在求解复矩阵时的应用。
### 2. jacobi迭代法的算法步骤在使用jacobi迭代法求解复矩阵特征值和特征向量时，我们需要经历一系列的算法步骤。
我们可以通过对角线元素的绝对值大小来判断矩阵是否已经对角化，然后进行迭代，直到满足精度要求为止。
### 3. jacobi迭代法的实际应用在实际应用中，jacobi迭代法除了可以求解复矩阵的特征值和特征向量外，还可以在解决其他涉及特征值和特征向量的问题时发挥重要作用。
通过简单的算法步骤和迭代过程，我们可以有效地得到复矩阵的特征值和特征向量，为进一步的分析和计算提供便利。
### 4. 个人观点和理解从个人的角度来看，jacobi迭代法在求解复矩阵特征值和特征向量时具有一定的优势，尤其在算法实现的过程中，我们可以通过简单的迭代步骤快速得到结果。
然而，对于大规模复矩阵的计算，可能还需要考虑其他更高效的方法或并行计算的应用。
### 结论通过本文的深入探讨，我们对jacobi迭代法求解复矩阵特征值和特征向量有了更深入的了解。
在实际应用中，我们需要灵活运用不同的方法和算法，以便更好地解决实际问题。
总结来说，jacobi迭代法是一种常用的求解复矩阵特征值和特征向量的方法，它通过简单的算法步骤和迭代过程，能够快速有效地得到结果。
然而，在实际应用中，我们还需要综合考虑不同的因素，以便获得更好的计算效果。
通过本文的阐述，希望读者能够更加深入地理解jacobi迭代法以及其在求解复矩阵特征值和特征向量中的应用，为进一步的学习和研究打下良好的基础。
雅可比方法求特征值
雅可比方法求特征值
雅可比方法是求解特征值和特征向量的一种常用方法。
其基本思想是通过不断地施加正交变换，将原矩阵转化为对角矩阵。
具体步骤如下：
1. 确定初始矩阵A和精度要求ε。
2. 构造正交矩阵P1，使得P1^T * A * P1的a_ij (i ≠ j)的绝对值小于精度要求ε。
3. 计算A1 = P1^T * A * P1，此时A1的对角线元素就是A的特征值。
4. 如果A1的对角线元素已经满足精度要求，则停止计算，否则重复第2步和第3步。
5. 将每个特征值代入A-lambda*I，其中I是单位矩阵，得到A-lambda*I的解集，解集中非零解的个数就是该特征值的重数。
对每个重数，求解A-lambda*I的零空间，就可以得到该特征值对应的特征向量。
需要注意的是，雅可比方法只适用于对称矩阵，对于非对称矩阵可以通过将其转化为对称矩阵来处理。
- 1 -。
雅可比迭代公式矩阵求法
雅可比迭代公式矩阵求法
在图形图像中很多地方用到求矩阵的特征值和特征向量，比如主成分分析、OBB包围盒等。
编程时一般都是用数值分析的方法来计算，这里介绍一下雅可比迭代法求解特征值和特征向量。
雅可比迭代法的原理，网上资料很多，详细可见参考资料1.这里我们简单介绍求解矩阵S特征值和特征向量的步骤：
1、初始化特征向量为对角阵V，即主对角线的元素都是1.其他元素为0.
2、在S的非主对角线元素中，找到绝对值最大元素Sij。
3、用下式计算tan2θ，求cosθ、sinθ及旋转矩阵Gij。
计算方法5矩阵特征值和特征向量的计算
A1
RART
a（1） 11 0
a（0212），
RT
c
s
s
c
RT的两个列向量是相应的特征向量。
1 0 2
例
考虑三阶矩阵
A 0
0
2
1
2 1 1
将A 中（3,1）和（1,3）位置上的元素变成0， 0
0.707 0 0.707
取R 0
0
1
0
0.707 0 0.707
做正交相似变换后得到
A-1u
k
，
1
因为A1 的计算
比较麻烦，而且往往不能保持矩阵A 的一些好性质
（如稀疏性），因此，反幂法在实际计算时以求解
方程组
Auk
u
k
，代替迭代
1
uk
A-1uk1求得uk，每
迭代一次要解一线性方程组。 由于矩阵在迭代过
程中不变，故可对A 先进行三角分解，每次迭代只 要解两个三角形方程组。
反幂法规范后的计算格式
满足 1 2 n 0 ，所对应的 n 个特征向量
x1 , x2 ,L , xn线性无关。 任取非零的初始向量u0，构造向量序列uk Auk1
向量uk逼近A的主特征值（按模最大的）对应的
特征向量，1
uk j uk-1 j
存在不全为零的常数（i i 1, 2,L , n）,（这里假设1 0），
Jacobi算法的基本思想:
d1
Ak +1 = Rk L
R2 R1 AR1T R2T L
RkT
=
RART
d2 O
dn
记R Rk L R2 R1，RT的每个列向量是对应的特征向量。
二、 雅可比过关法
雅克比矩阵知识介绍
雅可比矩阵(Jacobi方法)Jacobi 方法Jacobi方法是求对称矩阵的全部特征值以及相应的特征向量的一种方法，它是基于以下两个结论1) 任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q,使得Q T AQ = diag(λ1 ,λ2,…,λn) (3.1)其中λi(i=1,2,…,n)是A的特征值，Q中各列为相应的特征向量。
2) 在正交相似变换下,矩阵元素的平方和不变。
即设A=(aij )n×n,Q交矩阵,记B=Q T AQ=(bij )n×n, 则Jacobi方法的基本思想是通过一次正交变换,将A中的一对非零的非对角化成零并且使得非对角元素的平方和减小。
反复进行上述过程,使变换后的矩阵的非对角元素的平方和趋于零，从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。
1 矩阵的旋转变换设A为n阶实对称矩阵，考虑矩阵易见 Vij(φ)是正交矩阵, 记注意到B=VijA的第i,j行元素以及的第i,j列元素为可得≠0，取φ使得则有如果aij对A(1)重复上述的过程，可得A(2) ,这样继续下去, 得到一个矩阵序列{A(k) }。
可以证明，虽然这种变换不一定能使矩阵中非对角元素零元素的个数单调增加，但可以保证非对角元素的平方和递减，我们以A与A(1)为例进行讨论。
设由式(3.4)可得这表明,在上述旋转变换下,非对角元素的平方和严格单调递减,因而由(3.2)可知，对角元素的平方和单调增加。
2. Jacobi方法通过一系列旋转变换将A变成A(k+1) ,求得A的全部特征值与特征向量的方法称为Jacobi方法。
计算过程如下1)令k=0, A(k) =A2) 求整数i,j, 使得3) 计算旋转矩阵4) 计算A(k+1)5) 计算6) 若E(A(k+1))<ε, 则为特征值，Q T = (V(0) V(1)…V(k+1))T的各列为相应的特征向量；否则，k+1=>k返回2,重复上述过程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1.//////////////////////////////////////////////////////////////////////
2.// 求实对称矩阵特征值与特征向量的雅可比法
3.//
4.// 参数：
5.// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
6.// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
7.// 数组dblEigenValue中第j个特征值对应的特征向量
8.// 3. int nMaxIt - 迭代次数，默认值为60
9.// 4. double eps - 计算精度，默认值为0.000001
10.//
11.// 返回值：BOOL型，求解是否成功
12.//////////////////////////////////////////////////////////////////////
13.BOOL CMatrix::JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
14.{
15.int i,j,p,q,u,w,t,s,l;
16.double fm,cn,sn,omega,x,y,d;
17.
18.if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
19.return FALSE;
20.
21.l=1;
22.for (i=0; i<=m_nNumColumns-1; i++)
23.{
24.mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
25.for (j=0; j<=m_nNumColumns-1; j++)
26.if (i!=j)
27.mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;//单位矩阵
28.}
29.
30.while (TRUE)
31.{
32.fm=0.0;
33.for (i=1; i<=m_nNumColumns-1; i++)
34.{
35.for (j=0; j<=i-1; j++)
36.{
37.d=fabs(m_pData[i*m_nNumColumns+j]);
38.if ((i!=j)&&(d>fm))
39.{
40.fm=d;
41.p=i;
42.q=j;
}//取绝对值最大的非对角线元素，并记住位置
43.}
44.}
45.
46.if (fm<eps)
47.{
48.for (i=0; i<m_nNumColumns; ++i)
49.dblEigenValue = GetElement(i,i);
50.return TRUE;
51.}
52.
53.if (l>nMaxIt)
54.return FALSE;
55.
56.l=l+1;
57.u=p*m_nNumColumns+q;
58.w=p*m_nNumColumns+p;
59.t=q*m_nNumColumns+p;
60.s=q*m_nNumColumns+q;
61.x=-m_pData;
62.y=(m_pData[s]-m_pData[w])/2.0;
63.omega=x/sqrt(x*x+y*y);
64.
65.if (y<0.0)
66.omega=-omega;
67.
68.sn=1.0+sqrt(1.0-omega*omega);
69.sn=omega/sqrt(2.0*sn);//sin
=sqrt(1.0-sn*sn);//cos
71.fm=m_pData[w];
72.m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData*omega;
73.m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData*omega;
74.m_pData=0.0;
75.m_pData[t]=0.0;
76.for (j=0; j<=m_nNumColumns-1; j++)
77.{
78.if ((j!=p)&&(j!=q))
79.{
80.u=p*m_nNumColumns+j; w=q*m_nNumColumns+j;
81.fm=m_pData;
82.m_pData=fm*cn+m_pData[w]*sn;
83.m_pData[w]=-fm*sn+m_pData[w]*cn;
84.}
85.}
86.for (i=0; i<=m_nNumColumns-1; i++)
87.{
88.if ((i!=p)&&(i!=q))
89.{
90.u=i*m_nNumColumns+p;
91.w=i*m_nNumColumns+q;
92.fm=m_pData;
93.m_pData=fm*cn+m_pData[w]*sn;
94.m_pData[w]=-fm*sn+m_pData[w]*cn;
95.}
96.}
97.
98.for (i=0; i<=m_nNumColumns-1; i++)
99.{
100.u=i*m_nNumColumns+p;
101.w=i*m_nNumColumns+q;
102.fm=mtxEigenVector.m_pData;
103.mtxEigenVector.m_pData=fm*cn+mtxEigenVector.m_pData[w]*sn; 104.mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn; 105.}
106.
107.}
108.
109.for (i=0; i<m_nNumColumns; ++i)
110.dblEigenValue = GetElement(i,i);
111.
112.return TRUE;
113.}。