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

实验名称实验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。

实验五矩阵的lu分解法，雅可比迭代法班级：学号：：实验五 矩阵的LU 分解法，雅可比迭代一、目的与要求：➢ 熟悉求解线性方程组的有关理论和方法；➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序；➢ 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。

二、实验内容：➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例➢ 列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示1) 消元过程对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。

③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i k j k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算kk ik ik a a l /=对j=l+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。

②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算ii ni j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下：#include <iostream>#include <cmath>using namespace std;void disp(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void disp(double* q,int n){cout<<"====================================="<<endl; for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMax(double** p,int start,int end){int max=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[max][start]))max=i;}return max;}void swapRow(double** p,int one,int other,int col){double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMax(p,i,row); //找列主元行号if(p[flag][i]==0) return false;swapRow(p,i,flag,col); //交换行for(int j=i+1;j<row;j++){double elem=p[j][i]/p[i][i]; //消元因子 p[j][i]=0;for(int k=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}return true;}double sumRow(double** p,double* q,int row,int col){ double sum=0;for(int i=0;i<col-1;i++){if(i==row)continue;sum+=(q[i]*p[row][i]);}return sum;}void back(double** p,int row,int col,double* q){for(int i=row-1;i>=0;i--){q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; }}int main(){cout<<"Input n:";int n; //方阵的大小cin>>n;double **p=new double* [n];for(int i=0;i<n;i++){p[i]=new double [n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){cout<<"奇异"<<endl;return 0;}double* q=new double[n];for(int i=0;i<n;i++)q[i]=0;back(p,n,n+1,q);disp(q,n);delete[] q;for(int i=0;i<n;i++)delete[] p[i];delete[] p;}1. 用列主元消去法解方程⎪⎩⎪⎨⎧=++-=++-=++035.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.1833.555x 2.304x 0.101x 321321321例2 解方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧=++++=++++=++++=++++=++++-12.041.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E2.92D 1.46C3.53B -20.071.0F 1.10E4.48D 1.21C 4.93B -32.041.0F 1.55E5.66D 2.40C 8.77B 计算结果如下B=-1.461954C= 1.458125D=-6.004824E=-2.209018F= 14.719421➢矩阵直接三角分解法算法：将方程组A x=b 中的A分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组A x=b化为解2个方程组Ly=b，Ux=y，具体算法如下：①对j=1,2,3,…,n计算jjau11=对i=2,3,…,n计算1111/aalii=②对k=1,2,3,…,n:a. 对j=k,k+1,…,n计算∑-=-=11kqqjkqkjkjulaub. 对i=k+1,k+2,…,n计算kkkqqkiqikikuulal/)(11∑-=-=③11by=，对k=2,3,…,n计算∑-=-=11kqqkqkkylby④nnnnuyx/=,对k=n-1,n-2,…,2,1计算kknkqqkqkkuxuyx/)(1∑+=-=注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n n n n a a a a a a a a a a a a 施行算法②，③，此时U 的第n+1列元素即为y 。

1.//////////////////////////////////////////////////////////////////////2.// 求实对称矩阵特征值与特征向量的雅可比法3.//4.// 参数：5.// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值6.// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与7.// 数组dblEigenValue中第j个特征值对应的特征向量8.// 3. int nMaxIt - 迭代次数，默认值为609.// 4. double eps - 计算精度，默认值为0.00000110.//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);//cos71.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.}。

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%。

1、用jacobi方法计算对称矩阵A的特征值和对应的特征向量。

function [k,Bk,V,D,Wc]=jacobite(A,jd,max1[n,n]=size(A;P0=eye(n;Vk=eye(n;Bk=A;k=1;state=1;while ((k<=max1&(state==1aij=abs(Bk-diag(diag(Bk;[m1 i]=max(abs(aij;[m2 j]=max(m1;i=i(j;Aij=(Bk-diag(diag(Bk;mk=m2*sign(Aij(i,j;Wc=m2;Dk=diag(diag(Bk;Pk=P0;c=(Bk(j,j-Bk(i,i/(2*Bk(i,j;t=sign(c/(abs(c+sqrt(1+c^2;pii=1/(sqrt(1+t^2;pij=t/(sqrt(1+t^2;Pk(i,i=pii;Pk(i,j=pij;Pk(j,j=pii;Pk(j,i=-pij;B1=Pk'*Bk;B2=B1*Pk;Vk=Vk*Pk;Bk=B2;state=0;if (Wc>jdstate=1;endk,k=k+1;Pk,Vk,Bk=B2,Wc,endif (k>max1disp('迭代次数k已经达到最大迭代次数ma1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下:'elsedisp('迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下:'endk=k-1;Bk=B2;V=Vk;D=diag(diag(Bk;Wc;[V1,D1]=eig(A,'nobalance'>> A=[12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12];>> [k,Bk,V,D,Wc]=jacobite(A,0.0001,100k =1Pk =0.7227 0.6912 0 0-0.6912 0.7227 0 00 0 1.0000 00 0 0 1.0000 Vk =0.7227 0.6912 0 0-0.6912 0.7227 0 00 0 1.0000 00 0 0 1.0000 Bk =65.5558 0 0.7858 -0.7227-0.0000 -46.5558 3.5189 -0.69120.7858 3.5189 5.0000 1.0000-0.7227 -0.6912 1.0000 12.0000 Wc = 56k =2Pk =1.0000 0 0 00 0.9977 0.0678 00 -0.0678 0.9977 00 0 0 1.0000 Vk =0.7227 0.6896 0.0468 0-0.6912 0.7210 0.0490 00 -0.0678 0.9977 00 0 0 1.0000 Bk =65.5558 -0.0533 0.7840 -0.7227-0.0533 -46.7948 0 -0.75740.7840 0.0000 5.2391 0.9509-0.7227 -0.7574 0.9509 12.0000 Wc = 3.5189k =3Pk =1.0000 0 0 00 1.0000 0 00 0 0.9906 0.13670 0 -0.1367 0.9906 Vk =0.7227 0.6896 0.0464 0.0064-0.6912 0.7210 0.0485 0.00670 -0.0678 0.9883 0.13640 0 -0.1367 0.9906Bk =65.5558 -0.0533 0.8754 -0.6088-0.0533 -46.7948 0.1035 -0.75020.8754 0.1035 5.1079 -0.0000-0.6088 -0.7502 -0.0000 12.1312 Wc = 0.9509k =4Pk =0.9999 0 -0.0145 00 1.0000 0 00.0145 0 0.9999 00 0 0 1.0000 Vk =0.7233 0.6896 0.0359 0.0064-0.6904 0.7210 0.0585 0.00670.0143 -0.0678 0.9882 0.1364-0.0020 0 -0.1367 0.9906 Bk =65.5685 -0.0518 -0.0000 -0.6087-0.0518 -46.7948 0.1043 -0.7502-0.0000 0.1043 5.0952 0.0088-0.6087 -0.7502 0.0088 12.1312 Wc = 0.8754k =5Pk =1.0000 0 0 00 0.9999 0 -0.01270 0 1.0000 00 0.0127 0 0.9999 Vk =0.7233 0.6896 0.0359 -0.0024-0.6904 0.7211 0.0585 -0.00250.0143 -0.0660 0.9882 0.1372-0.0020 0.0126 -0.1367 0.9905 Bk = 65.5685 -0.0595 -0.0000 -0.6080-0.0595 -46.8044 0.1044 0.0000-0.0000 0.1044 5.0952 0.0075-0.6080 0.0000 0.0075 12.1407 Wc = 0.7502k =6Pk =0.9999 0 0 0.01140 1.0000 0 00 0 1.0000 0-0.0114 0 0 0.9999 Vk =0.7233 0.6896 0.0359 0.0059-0.6903 0.7211 0.0585 -0.01030.0127 -0.0660 0.9882 0.1374-0.0132 0.0126 -0.1367 0.9905 Bk = 65.5754 -0.0595 -0.0001 -0.0000-0.0595 -46.8044 0.1044 -0.0007-0.0001 0.1044 5.0952 0.0075-0.0000 -0.0007 0.0075 12.1338 Wc =0.6080k =7Pk =1.0000 0 0 00 1.0000 0.0020 00 -0.0020 1.0000 00 0 0 1.0000 Vk =0.7233 0.6895 0.0373 0.0059-0.6903 0.7209 0.0600 -0.01030.0127 -0.0680 0.9881 0.1374-0.0132 0.0129 -0.1366 0.9905 Bk = 65.5754 -0.0595 -0.0002 -0.0000-0.0595 -46.8046 -0.0000 -0.0007-0.0002 -0.0000 5.0954 0.0075-0.0000 -0.0007 0.0075 12.1338 Wc = 0.1044k =8Pk =1.0000 0.0005 0 0-0.0005 1.0000 0 00 0 1.0000 00 0 0 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9881 0.1374-0.0133 0.0129 -0.1366 0.9905 Bk = 65.5754 0.0000 -0.0002 0.0000-0.0000 -46.8046 -0.0000 -0.0007-0.0002 -0.0000 5.0954 0.00750.0000 -0.0007 0.0075 12.1338 Wc = 0.0595k =9Pk =1.0000 0 0 00 1.0000 0 00 0 1.0000 0.00110 0 -0.0011 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9880 0.1384-0.0133 0.0129 -0.1377 0.9903 Bk = 65.5754 0.0000 -0.0002 0.0000-0.0000 -46.8046 0.0000 -0.0007-0.0002 0.0000 5.0954 0.00000.0000 -0.0007 -0.0000 12.1338 Wc = 0.0075k =10Pk =1.0000 0 0 00 1.0000 0 -0.00000 0 1.0000 00 0.0000 0 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9880 0.1384-0.0133 0.0129 Bk = 65.5754 0.0000 0.0000 -46.8046 -0.0002 0.0000 0.0000 -0.0000 Wc = 6.9206e-004 k= 11 Pk = 1.0000 0 0 1.0000 -0.0000 0 0 0 Vk = 0.72290.6899 -0.6907 0.7206 0.0128 -0.0680 -0.0133 0.0129 Bk = 65.5754 -0.0000 -0.0000 -46.8046 -0.0000 0.0000 0.0000 -0.0000 Wc = 2.0482e-004 k= 12 Pk = 1.0000 0 0 1.0000 0 -0.0000 0 0 Vk = 0.7229 0.6899 -0.6907 0.7206 0.0128 -0.0680 -0.0133 0.0129 Bk = 65.5754 -0.0000 -0.0000 -46.8046 -0.0000 0.0000 0.0000 -0.0000 -0.1377 -0.00020.0000 5.0954 -0.0000 0.9903 0.0000 0.0000 -0.0000 12.1338 0.0000 0 1.0000 0 0.0373 0.0600 0.9880 -0.1377 0.0000 0.0000 5.0954 -0.0000 0 0 0 1.0000 0.0059 -0.0103 0.1384 0.9903 0.0000 0.0000 -0.0000 12.1338 0 0.0000 1.0000 0 0.0373 0.0600 0.9880 -0.1377 0.0000 -0.0000 5.0954 -0.0000 0 0 0 1.0000 0.0059 -0.0103 0.1384 0.9903 0.0000 0.0000 -0.0000 12.1338Wc = 6.2740e-007 迭代次数 k，对称矩阵 Bk，以特征向量为列向量的矩阵 V，特征值为对角元的对角矩阵 D 如下： V1 = 0.6899 -0.0373 0.0059 -0.7229 0.7206 -0.0600 -0.0103 0.6907 -0.0680 -0.9880 0.1384 -0.0128 0.0129 0.1377 0.9903 0.0133 D1 = -46.8046 0 0 0 0 5.0954 0 0 0 0 12.1338 0 0 0 0 65.5754 k= 12 Bk = 65.5754 -0.0000 0.0000 0.0000 -0.0000 -46.8046 -0.0000 0.0000 -0.0000 0.0000 5.0954 -0.0000 0.0000 -0.0000 -0.0000 12.1338 V= 0.7229 0.6899 0.0373 0.0059 -0.6907 0.7206 0.0600 -0.0103 0.0128 -0.0680 0.9880 0.1384 -0.0133 0.0129 -0.1377 0.9903 D= 65.5754 0 0 0 0 -46.8046 0 0 0 0 5.0954 0 0 0 0 12.1338 Wc = 6.2740e-007。

雅克⽐（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方法是一种通过迭代对矩阵进行相似变换，使得原矩阵逐步转化为对角矩阵的方法。

通过数值迭代，可以逐步逼近矩阵的特征值和对应的特征向量。

其基本原理如下：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。

#ifnd‎e f ks‎f loat‎‎#defi‎n e ks‎f loat‎ dou‎b le#‎e ndif‎#if‎n def ‎i nt16‎‎#defi‎n e in‎t16 ‎i nt#‎e ndif‎#if‎n def ‎u int1‎6‎#def‎i ne u‎i nt16‎unsi‎g ned ‎i nt#‎e ndif‎#if‎n def ‎k snew‎‎#def‎i ne k‎s new(‎s,t) ‎((t*)‎m allo‎c((s)‎*size‎o f(t)‎))#e‎n dif‎#ifn‎d ef k‎s free‎‎#defi‎n e ks‎f ree(‎p) if‎(p!=0‎){fre‎e(p);‎p=0;}‎#end‎i f/‎*‎ksJa‎c obiG‎e tRea‎l SymM‎a trix‎F eatu‎r eVal‎u e‎用ja‎c obi方‎法求取*实‎对称*矩阵‎的特征值和‎特征向量‎‎ aM‎a trix‎[in]:‎ n 阶‎实对称阵‎‎ si‎d e[in‎]: 矩阵‎阶大小‎‎ eps‎[in]:‎控制精度‎‎ m‎a xTim‎e s[in‎]: 最大‎迭代次数‎返‎回值：‎‎前si‎d e*si‎d e个单元‎存储特征向‎量，一列为‎一特征向量‎；后sid‎e个存储单‎元为特征值‎，‎‎也就是总大‎小为 si‎z e = ‎s ide*‎(side‎+1) 存‎储单元*‎/ksf‎l oat ‎*ksJa‎c obiG‎e tRea‎l SymM‎a trix‎F eatu‎r eVal‎u e( k‎s floa‎t aMa‎t rix[‎] , u‎i nt16‎side‎,ks‎f loat‎eps ‎, uin‎t16 m‎a xTim‎e s )‎{‎ ksf‎l oat ‎*t, *‎s, *s‎1, *v‎e ctor‎, *q1‎;‎ flo‎a t a1‎, c ‎, s2 ‎, t1 ‎, max‎E ps =‎0 , ‎d1 , ‎d2, r‎= 1;‎‎u int1‎6 i ,‎j , ‎p , q‎, m ‎, tim‎e s = ‎0 , l‎e n = ‎0;‎ ui‎n t16 ‎s tart‎I ndex‎= 0 ‎, end‎I ndex‎;‎ end‎I ndex‎= si‎d e + ‎s tart‎I ndex‎;‎side‎+= s‎t artI‎n dex;‎‎len ‎= sid‎e*sid‎e;‎ t ‎= ksn‎e w(le‎n , k‎s floa‎t);‎ s ‎= ks‎n ew(l‎e n , ‎k sflo‎a t);‎ s‎1 = k‎s new(‎l en ,‎ksfl‎o at);‎‎q1 = ‎k snew‎(len ‎, ksf‎l oat)‎;‎vect‎o r = ‎k snew‎(len ‎+ sid‎e , k‎s floa‎t);‎ me‎m set ‎( vec‎t or ,‎0 , ‎(len+‎s ide)‎*size‎o f(ks‎f loat‎));‎ fo‎r( i ‎= sta‎r tInd‎e x ; ‎i < e‎n dInd‎e x ; ‎i++)‎ {‎‎ v‎e ctor‎[i*si‎d e+i]‎= 1;‎//将初始‎Q[i][‎j]置为单‎位矩阵‎ }‎‎w hile‎(r >=‎eps ‎&& ti‎m es <‎maxT‎i mes)‎‎{‎‎p = q‎= st‎a rtIn‎d ex; ‎m axEp‎s = 0‎;‎‎f or(i‎= st‎a rtIn‎d ex; ‎i < e‎n dInd‎e x; i‎++)‎‎ {‎‎‎for(‎j = s‎t artI‎n dex;‎j < ‎e ndIn‎d ex; ‎j++)‎‎‎ {‎‎‎‎i f( (‎j != ‎i) &&‎fabs‎(aMat‎r ix[i‎*side‎+j]) ‎> max‎E ps) ‎‎‎‎ {‎‎‎‎ m‎a xEps‎= (f‎l oat)‎f abs(‎a Matr‎i x[i*‎s ide+‎j]);/‎/找非对角‎元素绝对值‎最大的元‎‎‎‎‎p = i‎; q =‎j;‎‎‎‎}‎‎ }‎‎ }‎‎‎r = m‎a xEps‎;//重置‎当前误差‎‎ a1‎= (f‎l oat)‎((aMa‎t rix[‎q*sid‎e+q]-‎a Matr‎i x[p*‎s ide+‎p])/(‎2*aMa‎t rix[‎p*sid‎e+q])‎);‎‎if(a‎1 >= ‎0)‎‎{‎‎‎t1 = ‎(floa‎t)(1/‎(fabs‎(a1)+‎s qrt(‎1+a1*‎a1)))‎;‎‎}‎‎e lse‎‎ {‎‎‎ t1 ‎= (fl‎o at)(‎-1/(f‎a bs(a‎1)+sq‎r t(1+‎a1*a1‎)));‎‎ }‎‎ c ‎= (fl‎o at)(‎1/sqr‎t(1+t‎1*t1)‎);‎‎s2 =‎t1*c‎;‎‎mems‎e t ( ‎s , 0‎, le‎n*siz‎e of(k‎s floa‎t));‎‎ fo‎r( i ‎= sta‎r tInd‎e x ; ‎i < e‎n dInd‎e x ; ‎i++)‎‎ {‎‎‎ s[i‎*side‎+i] =‎1;//‎将初始Q[‎i][j]‎置为单位矩‎阵‎‎}‎‎s[p*‎s ide+‎p] = ‎c ;s‎[p*si‎d e+q]‎=s2;‎‎ s[‎q*sid‎e+p] ‎= -s2‎;s[q*‎s ide+‎q]=c;‎‎ f‎o r(i ‎= sta‎r tInd‎e x; i‎< en‎d Inde‎x; i+‎+)‎‎{‎‎‎f or(j‎= st‎a rtIn‎d ex; ‎j < e‎n dInd‎e x; j‎++)‎‎‎{‎‎‎ s‎1[i*s‎i de+j‎] = s‎[j*si‎d e+i]‎;//将矩‎阵s1[i‎][j]化‎为s[i]‎[j]的转‎置‎‎ }‎‎ }‎‎‎f or(i‎= st‎a rtIn‎d ex; ‎i < e‎n dInd‎e x; i‎++)‎‎ {‎‎‎for(‎j = s‎t artI‎n dex;‎j < ‎e ndIn‎d ex; ‎j++)‎‎‎ {‎‎‎‎d1 = ‎0;‎‎‎ f‎o r(m ‎= sta‎r tInd‎e x; m‎< en‎d Inde‎x; m+‎+)‎‎‎ {‎‎‎‎‎d1 +‎= (fl‎o at)(‎s1[i*‎s ide+‎m]*aM‎a trix‎[m*si‎d e+j]‎);//计‎算s1[i‎][j]*‎a Matr‎i x[i]‎[j] ‎‎‎‎}‎‎‎ t[‎i*sid‎e+j] ‎= d1;‎‎‎ }‎‎ }‎‎ f‎o r(i ‎= sta‎r tInd‎e x; i‎< en‎d Inde‎x; i+‎+)‎‎{‎‎‎f or(j‎= st‎a rtIn‎d ex; ‎j < e‎n dInd‎e x; j‎++)‎‎‎{‎‎‎ d‎1 = d‎2 = 0‎;‎‎‎ fo‎r(m =‎star‎t Inde‎x; m ‎< end‎I ndex‎; m++‎)‎‎‎ {‎‎‎‎‎d1 +=‎(flo‎a t)(t‎[i*si‎d e+m]‎*s[m*‎s ide+‎j]);/‎/计算t[‎i][j]‎*s[i]‎[j]‎‎‎‎ d‎2 += ‎(floa‎t)(ve‎c tor[‎i*sid‎e+m]*‎s[m*s‎i de+j‎]);//‎计算Q[i‎][j]*‎s[i][‎j]‎‎‎ }‎‎‎‎ aMa‎t rix[‎i*sid‎e+j] ‎= d1;‎‎‎‎ q1[‎i*sid‎e+j] ‎= d2;‎‎‎ }‎‎ }‎‎ m‎e mcpy‎( ve‎c tor ‎, q1 ‎, len‎*size‎o f(ks‎f loat‎));‎‎ tim‎e s +=‎1;‎ }‎‎// 对特‎征值与特征‎向量进行整‎合‎for(‎i = s‎t artI‎n dex;‎i < ‎e ndIn‎d ex; ‎i++)‎ {‎‎ v‎e ctor‎[len+‎i] = ‎a Matr‎i x[i*‎s ide+‎i];‎ }‎‎k sfre‎e ( t‎);‎ ks‎f ree ‎( s )‎;‎ksfr‎e e ( ‎s1 );‎‎k sfre‎e ( q‎1 );‎‎i f(ti‎m es >‎maxT‎i mes)‎‎{‎‎k sfre‎e ( v‎e ctor‎); ‎‎ ret‎u rn N‎U LL;‎ }‎‎retu‎r n ve‎c tor;‎}‎。

雅可比方法求特征值matlab雅可比方法是一种用于求解矩阵的特征值和特征向量的有效算法。

在这个方法中，我们首先将矩阵表示为一个对称矩阵，并通过迭代的方式逐步逼近特征值和特征向量的精确解。

我们需要将矩阵进行对称化处理。

这可以通过将矩阵与其转置相乘来实现。

得到的对称矩阵将具有相同的特征值和特征向量。

接下来，我们需要选择一个初始的特征向量估计值。

这个初始估计可以是任意的非零向量。

然后，我们通过将矩阵与初始估计向量相乘，得到一个新的向量。

这个新向量将与初始估计向量正交，并且具有相同的特征值。

然后，我们需要对新向量进行归一化处理，以确保它具有单位长度。

这可以通过将新向量除以其范数来实现。

接下来，我们需要计算新向量与矩阵的乘积。

这将得到一个新的特征向量估计值。

然后，我们需要对新的特征向量进行归一化处理，并计算其与前一次迭代的特征向量之间的夹角。

如果夹角的差异小于某个预定的阈值，我们就可以认为特征向量已经收敛，并得到了矩阵的特征值和特征向量。

我们可以通过对特征值进行排序，以及相应地对特征向量进行排列，来得到特征值的精确解。

雅可比方法是一种非常强大的求解特征值和特征向量的方法。

它在各种领域中都有广泛的应用，包括物理学、工程学和计算机科学等。

在MATLAB中，我们可以使用eig函数来实现雅可比方法，并得到矩阵的特征值和特征向量。

雅可比方法是一种求解特征值和特征向量的有效算法。

它通过迭代的方式逐步逼近精确解，并在MATLAB中有着广泛的应用。

无论是在学术研究还是实际应用中，雅可比方法都是一个非常有用的工具。

雅可比迭代公式矩阵求法
在图形图像中很多地方用到求矩阵的特征值和特征向量，比如主成分分析、OBB包围盒等。

编程时一般都是用数值分析的方法来计算，这里介绍一下雅可比迭代法求解特征值和特征向量。

雅可比迭代法的原理，网上资料很多，详细可见参考资料1.这里我们简单介绍求解矩阵S特征值和特征向量的步骤：
1、初始化特征向量为对角阵V，即主对角线的元素都是1.其他元素为0.
2、在S的非主对角线元素中，找到绝对值最大元素Sij。

3、用下式计算tan2θ，求cosθ、sinθ及旋转矩阵Gij。

实验五矩阵的lu 分解法，雅可比迭代法班级：**：：实验五 矩阵的LU 分解法，雅可比迭代一、目的与要求：➢ 熟悉求解线性方程组的有关理论和方法；➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序；➢ 通过实际计算，进一步了解各种方法的优缺点，选择适宜的数值方法。

二、实验内容：➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例➢ 列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示1) 消元过程对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序完毕；否则执行③。

③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i k j k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算对j=l+1, ┅,n+1计算2) 回代过程①假设0=nn a ，则矩阵A 奇异，程序完毕；否则执行②。

②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算程序与实例程序设计如下：#include <iostream>#include <cmath>using namespace std;void disp(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void disp(double* q,int n){cout<<"====================================="<<endl; for(int i=0;i<n;i++)cout<<"*["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMa*(double** p,int start,int end){int ma*=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[ma*][start]))ma*=i;}return ma*;}void swapRow(double** p,int one,int other,int col){double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMa*(p,i,row); //找列主元行号if(p[flag][i]==0) return false;swapRow(p,i,flag,col); //交换行for(int j=i+1;j<row;j++){double elem=p[j][i]/p[i][i]; //消元因子 p[j][i]=0;for(int k=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}return true;}double sumRow(double** p,double* q,int row,int col){ double sum=0;for(int i=0;i<col-1;i++){if(i==row)continue;sum+=(q[i]*p[row][i]);}return sum;}void back(double** p,int row,int col,double* q){for(int i=row-1;i>=0;i--){q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; }}int main(){cout<<"Input n:";int n; //方阵的大小cin>>n;double **p=new double* [n];for(int i=0;i<n;i++){p[i]=new double [n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){cout<<"奇异"<<endl;return 0;}double* q=new double[n];for(int i=0;i<n;i++)q[i]=0;back(p,n,n+1,q);disp(q,n);delete[] q;for(int i=0;i<n;i++)delete[] p[i];delete[] p;}1. 用列主元消去法解方程例2 解方程组计算结果如下➢ 矩阵直接三角分解法算法：将方程组A *=b 中的A 分解为A =LU ，其中L 为单位下三角矩阵，U 为上三角矩阵，则方程组A *=b 化为解2个方程组Ly =b ，U*=y ，具体算法如下：①对j=1,2,3,…,n 计算对i=2,3,…,n 计算②对k=1,2,3,…,n:a. 对j=k,k+1,…,n 计算b. 对i=k+1,k+2,…,n 计算③11b y =，对k=2,3,…,n 计算④nn n n u y x /=,对k=n-1,n-2,…,2,1计算注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n n n n a a a a a a a a a a a a 施行算法②，③，此时U 的第n+1列元素即为y 。

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方法是求对称矩阵的全部特征值以及相应的特征向量的一种方法，它是基于以下两个结论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,重复上述过程。

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迭代法的原理、应用以及个人观点和理解。

### 1. jacobi迭代法的原理和概念jacobi迭代法是一种通过不断相似变换将矩阵对角化的方法，它可以被用于求解实对称矩阵的特征值和特征向量，而在这篇文章中，我们将着重讨论其在求解复矩阵时的应用。

### 2. jacobi迭代法的算法步骤在使用jacobi迭代法求解复矩阵特征值和特征向量时，我们需要经历一系列的算法步骤。

我们可以通过对角线元素的绝对值大小来判断矩阵是否已经对角化，然后进行迭代，直到满足精度要求为止。

### 3. jacobi迭代法的实际应用在实际应用中，jacobi迭代法除了可以求解复矩阵的特征值和特征向量外，还可以在解决其他涉及特征值和特征向量的问题时发挥重要作用。

通过简单的算法步骤和迭代过程，我们可以有效地得到复矩阵的特征值和特征向量，为进一步的分析和计算提供便利。

### 4. 个人观点和理解从个人的角度来看，jacobi迭代法在求解复矩阵特征值和特征向量时具有一定的优势，尤其在算法实现的过程中，我们可以通过简单的迭代步骤快速得到结果。

然而，对于大规模复矩阵的计算，可能还需要考虑其他更高效的方法或并行计算的应用。

### 结论通过本文的深入探讨，我们对jacobi迭代法求解复矩阵特征值和特征向量有了更深入的了解。

在实际应用中，我们需要灵活运用不同的方法和算法，以便更好地解决实际问题。

总结来说，jacobi迭代法是一种常用的求解复矩阵特征值和特征向量的方法，它通过简单的算法步骤和迭代过程，能够快速有效地得到结果。

然而，在实际应用中，我们还需要综合考虑不同的因素，以便获得更好的计算效果。

通过本文的阐述，希望读者能够更加深入地理解jacobi迭代法以及其在求解复矩阵特征值和特征向量中的应用，为进一步的学习和研究打下良好的基础。