雅可比迭代法与矩阵的特征值

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.}。

准备工作➢算法设计矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。

分析一：由题目中所给条件λ1≤λ2≤…≤λn，可得出λ1、λn按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现|λs|<λ1|<|λn |或|λs|<λn|<|λ1 |的情况，导致按幂法和反幂法无法求解λ1或λn二者中的一者；分析二：题目要求求解与数μk =λ1+k(λn-λ1)/40最接近的特征值λik(k=1,2,3…39)，这个问题其实可以转换为求A-μk 按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得λ1和λn，所以第二个问题暂先搁浅；分析三：cond(A) 2 = ||A|| * ||A-1|| =|λ|max * |λ|min,这可以用幂法和反幂法求得，det(A) =λ1 *λ2 * … *λn，这需要求得矩阵A的所有特征值。

由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。

所以该题可以用Jacobi法求解。

➢模块设计由➢数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。

但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。

基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。

完整代码如下（编译环境windows10 + visual studio2010）：完整代码// math.cpp : 定义控制台应用程序的入口点。

//#include "stdafx.h"#include<stdio.h>#include<math.h>#include<time.h>#define N 501#define V (N+1)*N/2+1#define e 2.6630353#define a(i) (1.64 - 0.024 * (i)) * sin(0.2 * (i)) - 0.64 * pow(e , 0.1 / (i)) #define b 0.16#define c -0.064#define eps pow((double)10.0,-12)#define PFbits "%10.5f "#define PFrols 5#define PFe %.11e#define FK 39int p;int q;double cosz;double sinz;double MAX;int kk;//#define PTS pts#ifdef PTSvoid PTS(double *m){printf("-----------------------------------------------------------------------\n");printf(" 迭代第%d次\n",kk);for(int i = 1 ; i <= PFrols ; i++){for( int j = (i-1)*i/2+1 ; j <= (i+1)*i/2 ; j++){printf(PFbits,m[j]);}putchar(10);}for(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 1 ; i <= PFrols+2 ; i++){printf(" ... ");}putchar(10);}#elsevoid PTS(double *m){}#endifvoid recounti(int i , int *pp, int *qq){for(int j = 0 ; j <= N-1 ; j++){if( (i - (j+1)*j/2) <= j+1){*pp = j+1;*qq = i - (j+1)*j/2;break;}}}void refreshMetrix(double *m){int ipr,ipc,iqr,iqc;m[(p+1)*p/2] = m[(p+1)*p/2] * pow(cosz,2) + m[(q+1)*q/2] * pow(sinz,2) + 2 * m[(p-1)*p/2+q] * cosz * sinz;m[(q+1)*q/2] = m[(p+1)*p/2] * pow(sinz,2) + m[(q+1)*q/2] * pow(cosz,2) - 2 * m[(p-1)*p/2+q] * cosz * sinz;for(int i = 1; i <= N ;i++){if(i != p && i != q){if(i > p){ipr = i;ipc = p;}else{ipr = p;ipc = i;}if(i > q){iqr = i;iqc = q;}else{iqr = q;iqc = i;}m[(ipr-1)*ipr/2+ipc] = m[(ipr-1)*ipr/2+ipc] * cosz + m[(iqr-1)*iqr/2+iqc] * sinz;m[(iqr-1)*iqr/2+iqc] = -m[(ipr-1)*ipr/2+ipc] * sinz + m[(iqr-1)*iqr/2+iqc] * cosz;}}m[(p-1)*p/2+q] = 0;PTS(m);}//void calCosSin(double *m){double app = m[(p+1)*p/2];double aqq = m[(q+1)*q/2];double apq = m[(p-1)*p/2+q];cosz = cos(atan(2 * apq / (app - aqq))/2);sinz = sin(atan(2 * apq / (app - aqq))/2); }//void find_pq(double *m){double max = 0.0;int pp = 0;int qq = 0;for(int i = 1 ; i <= V ; i++){if(fabs(m[i]) > max){recounti(i,&pp,&qq);if(pp != qq){max = fabs(m[i]);p = pp;q = qq;}}}MAX = max;}void init(double *m){for(int i = 1 ; i <= N ;i++)m[(i+1)*i/2] = a(i);for(int i = 2 ; i <= N ; i++)m[(i-1)*i/2+i-1] = b;for(int i = 3 ; i <= N ; i++)m[(i-1)*i/2+i-2] = c;PTS(m);}void calFinal(double *m){printf("---------------------------------------------------------------------------------------------------\n");printf("结果输出:\n\n");double conda;double deta = 1.0;double minlumda = pow((double)10.0,12);double maxlumda = pow((double)10.0,-12);double absminlumda = pow((double)10.0,12);for(int i = 1 ; i <=N ;i++){if(m[(i+1)*i/2] > maxlumda)maxlumda = m[(i+1)*i/2];if(m[(i+1)*i/2] < minlumda)minlumda = m[(i+1)*i/2];if(fabs(m[(i+1)*i/2]) < absminlumda)absminlumda = fabs(m[(i+1)*i/2]);deta *= m[(i+1)*i/2];}if(fabs(minlumda) < fabs(maxlumda))conda = fabs(maxlumda) / absminlumda;elseconda = fabs(minlumda) / absminlumda;printf(" Lumda(1)=%.11e Lumda(%d)=%.11e Lumda(s)=%.11e\n",minlumda,N,maxlumda,absminlumda);printf(" Cond(A)=%.11e\n",conda);printf(" Det(A)=%.11e\n\n",deta);for(int i = 1 ; i <= FK ; i++){double muk = minlumda + i * (maxlumda - minlumda) / 40;double lumdak = 0.0;double tempabsmin = pow((double)10.0,12);for(int j = 1 ; j <= N ;j++){if(fabs(muk - m[(j+1)*j/2]) < tempabsmin){lumdak = m[(j+1)*j/2];tempabsmin = fabs(muk - m[(j+1)*j/2]);}}printf(" Lumda(i%d)=%.11e ",i,lumdak);if(i%3==0)putchar(10);}putchar(10);printf("------------------------------------------------------------------------------------------------------\n");putchar(10);putchar(10);}int _tmain(int argc, _TCHAR* argv[]){double m[(N+1)*N/2+1] = {0.0};kk=0;MAX=1.0;time_t t0,t1;t0 = time( &t0);init(m);#ifndef PTSprintf("正在计算...\n\n");#endifwhile(true){kk++;find_pq(m);if(MAX<eps)break;#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n"); #endifcalCosSin(m);refreshMetrix(m);}#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n");#endifprintf("矩阵最终形态...\n");for(int i = 1 ; i <= PFrols ; i++){for( int j = (i-1)*i/2+1 ; j <= (i+1)*i/2 ; j++){printf(PFbits,m[j]);}putchar(10);}for(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 1 ; i <= PFrols+2 ; i++){printf(" ... ");}putchar(10);t1 = time(&t1);#ifdef PTSprintf("计算并输出用时%.2f秒\n\n",difftime(t1,t0));#elseprintf("迭代次数%d,计算用时%.2f秒\n\n",kk,difftime(t1,t0)); #endifcalFinal(m);return 0;}运行结果如下：中间运行状态如下：结果分析➢有效性分析1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2.代码中有定义预编译宏//#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G存，64位windows10操作系统）。

雅可⽐算法求矩阵的特征值和特征向量⽬的求⼀个实对称矩阵的所有特征值和特征向量。

前置知识对于⼀个实对称矩阵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%。

类矩阵两种迭代法的收敛性比较引言：在科学计算中，线性方程组的求解是很普遍的问题。

尤其是在大型科学计算中，线性方程组的求解是最重要的任务之一。

线性方程组的求解有很多种方法，例如高斯消元法、LU分解法、迭代法等等，其中迭代法是一种高效的方法。

迭代法的思想是从一个初值解开始，逐步改进解的准确度，直到满足误差要求。

在本文中，我们将讨论两种类矩阵迭代法的收敛性比较，即雅可比迭代法和高斯-赛德尔迭代法。

1.雅可比迭代法（Jacobi Iterative Method）：雅可比迭代法是最简单的迭代法之一。

它是基于线性方程组的矩阵形式 Ax=b，将 A 分解成 A=D-L-U（D为A的对角线元素，L为A的下三角矩阵，U为A的上三角矩阵），其中 D 为对角线元素，L为严格下三角矩阵，U 为严格上三角矩阵。

则有如下迭代关系式： x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b (1)其中，x^{(k)} 为 k 次迭代后的解，x^{(0)} 为初始解。

雅可比迭代法的迭代矩阵为M = D^{-1}(L+U)。

以下是雅可比迭代法的收敛性分析：定理1：若矩阵 A 为对称正定矩阵，则雅可比迭代法收敛。

证明：由于 A 为对称正定矩阵，所以存在唯一的解。

假设迭代后得到的解为 x^{(k)}，则我们可以用误差向量 e^{(k)} = x-x^{(k)} 表示剩余项，则有 Ax^{(k)}-b = e^{(k)}。

对 (1) 式两边同时乘以 A^-1，得：x^{(k+1)}=x^{(k)}-A^{-1}e^{(k)}。

(2)将 (2) 式代入 Ax^{(k)}-b = e^{(k)} 中，得：Ax^{(k+1)}-b = Ae^{(k)}.(3)由于 A 为对称正定矩阵，则存在 A=Q\\Lambda Q^{-1}，其中Q 为正交矩阵，\\Lambda 为对角矩阵。

因此，我们可以将 (3) 式转化为：\\| x^{(k+1)}-x \\|_{A} =\\| Q^{-1}A^{-1}Qe^{(k)}\\|_{\\Lambda} \\leq \\rho (Q^{-1}A^{-1}Q)\\|e^{(k)}\\|_{A}。

矩阵特征值问题的计算方法特征值问题：A V=λV¾直接计算：A的阶数较小,且特征值分离得较好 特征值：det(λI-A)=0,特征向量：(λI-A)V=0¾迭代法：幂法与反幂法¾变换法：雅可比方法与QR方法内容：一、 特征值的估计及其误差问题二、 幂法与反幂法三、 雅可比方法四、 QR方法一、 特征值的估计及其误差问题 （一）特征值的估计结论 1.1：n 阶矩阵()ij n n A a ×=的任何一个特征值必属于复平面上的n 个圆盘：1,||||,1,2,ni ii ij j j i D z z a a i n =≠⎧⎫⎪⎪=−≤=⎨⎬⎪⎪⎩⎭∑"（10.1） 的并集。

结论1.2：若（10.1）中的m个圆盘形成一个连通区域D，且D与其余的n-m个圆盘不相连，则D中恰有A的m个特征值。

（二）特征值的误差问题结论1.3：对于n 阶矩阵()ij n n A a ×=，若存在n 阶非奇异矩阵H ，使得11(,,)n H AH diag λλ−=Λ="， （10.2）则11min ||||||||||||||i p p p i nH H A λλ−≤≤−≤∆ （10.3）其中λ是A A +∆的一个特征值，而(1,,)i i n λ="是A 的特征值，1,2,p =∞。

结论1.4：若n 阶矩阵A 是实对称的，则1min ||||||i p i nA λλ≤≤−≤∆。

（10.4）注：（10.4）表明，当A 是实对称时，由矩阵的微小误差所引起的特征值摄动也是微小的。

但是对于非对称矩阵而言，特别是对条件数很大的矩阵，情况未必如此。

二、 幂法与反幂法（一） 幂法：求实矩阵按模最大的特征值与特征向量假设n 阶实矩阵A 具有n 个线性无关的特征向量,1,iV i n ="，则对于任意的0nX R ∈，有 01ni ii X a V ==∑，从而有01111112((/))n nk k k i i i i ii i nk k i i i i A X a A V a V a V a V λλλλ======+∑∑∑.若A 的特征值分布如下：123||||||||n λλλλ>≥≥≥"，则有01111()k kk A X a V λλ→∞⎯⎯⎯→为对应的特征向量须注意的是，若1||1λ<，则10kλ→，出现“下溢”，若1||1λ>，则1kλ→∞，出现“上溢”，为避免这些现象的发生，须对0kA X 进行规范化。

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。

第4章 线性方程组和矩阵特征值的迭代解法线性代数计算方法中的迭代解法（即迭代法）是一类重要方法。

其基本思想是构造适当的矩阵序列或向量序列，使其逐步逼近所求问题的精确解，故又称矩阵迭代方法。

在求解阶数较高且零系数较多的大型稀疏线性代数方程组时，迭代法是很有效的。

矩阵特征值问题的求解通常也要用迭代法。

本章着重介绍求解线性代数方程组常用的简单迭代法及其收敛条件，并对计算矩阵特征值问题的雅可比方法和QR 方法作一些介绍。

4.1 线性代数方程组的迭代解法线性方程组(3.1)的迭代解法其基本思想与一元非线性方程的迭代解法类似，即构造适当的迭代公式，任选一个初始向量)0(x 进行迭代计算，使生成的向量序列,,)1()0(x x …，)(k x ,…收敛于方程组的精确解。

4.1.1 简单迭代法的一般形式设方程组(3.1)的系数矩阵非奇异，把它化为等价的方程组g Mx x += (4.1)其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n n nn n n n n x x x g g g m m m m m m m m m M M ΛΛΛΛ2121212222111211,,x g M按(4.1)构造迭代公式Λ,2,1,0,)()1(=+=+k k k g x M x (4.2)其中),1,0(],,,[T)()(2)(1)(ΛΛ==k x x x k n k k k x 。

任取初始向量)0(x ，用(4.2)逐次计算近似解向量,,x ,,x ,x ΛΛ)()2()1(k 这种方法称为简单迭代法，称(4.2)为简单迭代公式，M 为迭代矩阵。

公式的分量形式是⎪⎪⎩⎪⎪⎨⎧++++=++++=++++=+++nk n nn k n k n k nk n n k k k k n n k k k g x m x m x m x g x m x m x m x g x m x m x m x )()(22)(11)1(2)(2)(222)(121)1(21)(1)(212)(111)1(1ΛΛΛΛΛΛΛ即 ),2,1,0(,,,2,1,1)()1(ΛΛ==+=∑=+k n i g x m x i nj k j ij k i (4.3)如果)(k x的各分量存在极限n i x x i k i k ,,2,1,lim )(Λ==*∞→ (4.4)则称向量序列}{)(k x 收敛于向量T 21][****=n x x x ,,,x Λ，并记为 *∞→=x x )(lim k k (4.5)这时，称简单迭代法(4.2)是收敛的，否则就是发散的。