雅克比过关法求特征值和特征向量
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.}。
《实验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;}。