第6章矩阵特征值及特征向量计算

1、第6章矩阵特征值及特征向量计算 第第6章章 矩阵的特征值及特征向矩阵的特征值及特征向 量的计算量的计算 第6章矩阵特征值及特征向量计算 6.1 引言引言 在实际的工程计算中，经常会遇到计算矩阵的特征值及特征向量在实际的工程计算中，经常会遇到计算矩阵的特征值及特征向量 的问题。在线性代数中我们知道求矩阵的问题。在线性代数中我们知道求矩阵 (6.1) 的特征值就是求代数方程的特征值就是求代数方程的根的根。 (6.2) nn2n1n n22221 n11211 aaa aaa aaa A 0 aaa aa aaa IA)( nn2n1n n22221 n11211 a 第6章矩阵特征值及特征向量计算

2、 是关于的次多项式，可以写成是关于的次多项式，可以写成 (6.3) (6.3) 其中其中 取决于矩阵取决于矩阵 A A 的元素，称的元素，称 是矩阵是矩阵 A A 的特征多的特征多 项式，项式， 为特征方程，它有为特征方程，它有n n个根（包括重根），称为个根（包括重根），称为 A A 的特征根或的特征根或 特征值。特征值。 当当 是矩阵是矩阵 A A 的特征值时，相应的方程组的特征值时，相应的方程组 的非零解的非零解 ，称为矩阵，称为矩阵 A A 关于关于 的特征向量。的特征向量。 从从(6.1)(6.1)式及式及(6.2)(6.2)式看，它只是代数方程求根及线性方程组求解的问题。当式看，它

3、只是代数方程求根及线性方程组求解的问题。当 很小时（如很小时（如 ），这种方法是可行的。但当），这种方法是可行的。但当 n n 稍大时，多项式方程稍大时，多项式方程 是一个高次方程，求解它是一个很困难的问题。是一个高次方程，求解它是一个很困难的问题。 本章主要介绍四种目前在计算机上比较常用的计算矩阵的特征值和特征向本章主要介绍四种目前在计算机上比较常用的计算矩阵的特征值和特征向 量的幂法、反幂法、雅可比法及雅可比过关法。量的幂法、反幂法、雅可比法及雅可比过关法。 )( 0)( 1 1 1 nn nn aaa )1,2,(n ,iai 0)( 0 xAI x 2,3,4n )( 第6章矩阵特征值

4、及特征向量计算 6.2 6.2 幂法和反幂法幂法和反幂法 6.2.1 6.2.1 幂法和反幂法的基本思想幂法和反幂法的基本思想 幂法基本思想幂法基本思想 幂法是计算矩阵按模最大的特征值和相应的特征向量幂法是计算矩阵按模最大的特征值和相应的特征向量 的数值计算方法。简单地说，对于某的数值计算方法。简单地说，对于某 阶矩阵阶矩阵 ，任取，任取 初始向量初始向量 ，进行如下迭代计算，进行如下迭代计算 得到迭代序列得到迭代序列 ，再分析，再分析 与与 之间的关之间的关 系，就可得到系，就可得到 的按模最大的特征值和相应的特征向量的按模最大的特征值和相应的特征向量 的近似解的近似解。 n (0) X )

5、(1)(kk AXX )(k X 1)( k X A )(k X 1)( k X A 第6章矩阵特征值及特征向量计算 反幂法的基本思想反幂法的基本思想 反幂法是计算矩阵按模最小的特征值和相应的特征向反幂法是计算矩阵按模最小的特征值和相应的特征向 量的数值计算方法。设某量的数值计算方法。设某 阶矩阵阶矩阵 可逆，可逆， 和和 分别分别 为为 的特征值和相应的特征向量，并设的特征值和相应的特征向量，并设 ， 对对 两边同乘两边同乘 ，得，得 ，可见，可见 和和 的的 特征值互为倒数，而且特征值互为倒数，而且 也是也是 的特征值的特征值 的特征向的特征向 量。量。 的按模最大的特征值正是的按模最大的

6、特征值正是 的按模最小的特征值的按模最小的特征值 的倒数，用幂法计算的倒数，用幂法计算 的按模最大的特征值而得到的按模最大的特征值而得到 的的 按模最小的特征值的方法，称为反幂法。按模最小的特征值的方法，称为反幂法。 n n, 2 , 1i , 0 i A A 1 1 A 1 A 1 A A 1 A 1 A 1 AA 1 AA 第6章矩阵特征值及特征向量计算 6.2.2 6.2.2 实现幂法和反幂法的基本步骤实现幂法和反幂法的基本步骤 实现幂法和反幂法的基本步骤如下实现幂法和反幂法的基本步骤如下: : （1 1） 输入矩阵的阶数输入矩阵的阶数 n n、矩阵矩阵A A、最大迭代次数最大迭代次数M

7、AX_NMAX_N以及容许误差以及容许误差； （2 2） 对对 求矩阵求矩阵A A的按模最大的特征值的按模最大的特征值 和相应的特征向量和相应的特征向量 （ 足够大）；利用反幂法求得其按模最小的特征值和相应的特征向量足够大）；利用反幂法求得其按模最小的特征值和相应的特征向量 （即求矩阵（即求矩阵 的按模最大的特征值的按模最大的特征值) )。其中。其中 是用直是用直 接分解法：用公式接分解法：用公式 计算计算 ； （3 3） 输出求得的按模最大和按模最小的特征值以及相应的特征向量；输出求得的按模最大和按模最小的特征值以及相应的特征向量； （4 4） 结束。结束。 ,k210 )(1)(kk Ay

8、x )1()1(1)(kkk xxy |x| k i ni )( 1 max )(k y k 1 A )(11)(kk yAx )(1)(kk yUxL )(1)(kk yAx 第6章矩阵特征值及特征向量计算 例例6.16.1 使用幂法计算下列矩阵的按模最大的特征值 和相应的特征向量，精度为 #include #include #define N 3 /* N为矩阵的阶数 */ #define eps 1e-05 /* 容许误差 */ #define MAX_N 100 /* 最大迭代次数 */ 163 4310 232 A 5 10 第6章矩阵特征值及特征向量计算 int mifa(doub

9、le aNN)int mifa(double aNN) double xN,nxN; double xN,nxN; double xmax,oxmax; double xmax,oxmax; int i,j,k,flag; int i,j,k,flag; flag=0; oxmax=0; flag=0; oxmax=0; for(i=0;iN;i+) xi=1.0; for(i=0;iN;i+) xi=1.0; for(i=0;iMAX_N;i+) for(i=0;iMAX_N;i+) for(j=0;jN;j+) / for(j=0;jN;j+) /* * 幂乘幂乘 * */ / nxj=0

10、;nxj=0; for(k=0;kN;k+) nxj=nxj+ajk for(k=0;kN;k+) nxj=nxj+ajk* *xk;xk; xmax=0.0; xmax=0.0; for(j=0;jN;j+) / for(j=0;jxmax) xmax=fabs(nxj); if(fabs(nxj)xmax) xmax=fabs(nxj); for(j=0;jN;j+) nxj/=xmax; for(j=0;jN;j+) nxj/=xmax; for(j=0;jN;j+) xj=nxj; for(j=0;jN;j+) xj=nxj; if(fabs(xmax-oxmax)eps) / if(

11、fabs(xmax-oxmax)eps) /* * 输出按模最大的特征值和特征向量输出按模最大的特征值和特征向量 * */ / printf(n Max EigenValuen%11lfn,xmax);printf(n Max EigenValuen%11lfn,xmax); printf(n Max EigenVectorn); printf(n Max EigenVectorn); for(i=0;iN;i+) printf(%11lfn,nxi); for(i=0;iN;i+) printf(%11lfn,nxi); flag=1; flag=1; break; break; oxmax

12、=xmax; oxmax=xmax; return(flag); return(flag); 第6章矩阵特征值及特征向量计算 main()main() double aNN=2.0,3.0,2.0,10.0,3.0,4.0,3.0,6.0,1.0; double aNN=2.0,3.0,2.0,10.0,3.0,4.0,3.0,6.0,1.0; int i,j,flag; int i,j,flag; clrscr(); clrscr(); printf( Matrixn); printf( Matrixn); for(i=0;iN;i+) / for(i=0;iN;i+) /* * 输出矩阵输

13、出矩阵 * */ / printf( );printf( ); for(j=0;jN;j+) for(j=0;jN;j+) printf(%11lf,aij); printf(%11lf,aij); printf(n); printf(n); printf(n); printf(n); flag=mifa(a); flag=mifa(a); if(flag=0) / if(flag=0) /* * 输出无解信息输出无解信息 * */ / printf(Return no solovedn);printf(Return no solovedn); 第6章矩阵特征值及特征向量计算 程序运行结果：程

14、序运行结果： MatrixMatrix 2.000000 3.000000 2.000000 2.000000 3.000000 2.000000 10.000000 3.000000 4.000000 10.000000 3.000000 4.000000 3.000000 6.000000 1.000000 3.000000 6.000000 1.000000 Max EigenValueMax EigenValue 11.000002 11.000002 Max EigenVectorMax EigenVector 0.500000 0.500000 1.000000 1.000000

15、0.750000 0.750000 第6章矩阵特征值及特征向量计算 例例6.26.2 使用反幂法计算下列矩阵的按模最小的特征 值和相应的特征向量，精度为 。 #include #include #define N 3 /* N为矩阵的阶数 */ #define eps 1e-05 /* 容许误差 */ #define MAX_N 100 /* 最大迭代次数 */ 163 4310 232 A 5 10 第6章矩阵特征值及特征向量计算 int fmifa(double aNN)int fmifa(double aNN) double lNN,uNN; double lNN,uNN; double

16、 xN,nxN; double xN,nxN; double xmax,oxmax; double xmax,oxmax; int i,j,k,flag; int i,j,k,flag; flag=0; flag=0; for(i=0;iN;i+) / for(i=0;iN;i+) /* * u u矩阵对交线元素为矩阵对交线元素为1 1 * */ / uii=1.0;uii=1.0; for(k=0;kN;k+) for(k=0;kN;k+) for(i=k;iN;i+) / for(i=k;iN;i+) /* * 计算计算l l矩阵矩阵 * * lik=aik;lik=aik; for(j=

17、0;j=k-1;j+) for(j=0;j=k-1;j+) lik=lik-(lij lik=lik-(lij* *ujk);ujk); for(j=k+1;jN;j+) / for(j=k+1;jN;j+) /* * 计算计算u u矩阵矩阵 * */ / ukj=akj;ukj=akj; for(i=0;i=k-1;i+) for(i=0;i=k-1;i+) ukj=ukj-(lki ukj=ukj-(lki* *uij);uij); ukj/=lkk; ukj/=lkk; 第6章矩阵特征值及特征向量计算 for(i=0;iN;i+) xi=1.0;for(i=0;iN;i+) xi=1.0

18、; oxmax=0.0;oxmax=0.0; for(i=0;iMAX_N;i+)for(i=0;iMAX_N;i+) for(j=0;jN;j+) / for(j=0;jN;j+) /* * 反幂乘反幂乘 * */ / nxj=xj;nxj=xj; for(k=0;k=j-1;k+) nxj=nxj-ljk for(k=0;k=0;j-) for(j=N-1;j=0;j-) xj=nxj; xj=nxj; for(k=j+1;kN;k+) xj=xj-ujk for(k=j+1;kN;k+) xj=xj-ujk* *xk;xk; xmax=0.0; xmax=0.0; for(j=0;jN;

19、j+) / for(j=0;jxmax) xmax=fabs(xj);if(fabs(xj)xmax) xmax=fabs(xj); for(j=0;jN;j+) xj/=xmax; for(j=0;jN;j+) xj/=xmax; if(fabs(xmax-oxmax)eps) / if(fabs(xmax-oxmax)eps) /* *输出按模最小的特征值和特征向量输出按模最小的特征值和特征向量 * */ / printf(n Min EigenValuen%11lfn,xmax);printf(n Min EigenValuen%11lfn,xmax); printf(n Min Eig

20、enVectorn); printf(n Min EigenVectorn); 第6章矩阵特征值及特征向量计算 for(i=0;iN;i+) for(i=0;iN;i+) printf(%11lfn,xi); printf(%11lfn,xi); flag=1; flag=1; break; break; oxmax=xmax; oxmax=xmax; return(flag); return(flag); 第6章矩阵特征值及特征向量计算 main()main() double aNN=2.0,3.0,2.0,10.0,3.0,4.0,3.0,6.0,1.0; double aNN=2.0,3

21、.0,2.0,10.0,3.0,4.0,3.0,6.0,1.0; int i,j,flag; int i,j,flag; clrscr(); clrscr(); printf( Matrixn); printf( Matrixn); for(i=0;iN;i+) / for(i=0;iN;i+) /* * 输出矩阵输出矩阵 * */ / printf( );printf( ); for(j=0;jN;j+) for(j=0;jN;j+) printf(%11lf,aij); printf(%11lf,aij); printf(n); printf(n); printf(n); printf(n

22、); flag=fmifa(a); flag=fmifa(a); if(flag=0) / if(flag=0) /* * 输出无解信息输出无解信息 * */ / printf(Return no solovedn);printf(Return no solovedn); 第6章矩阵特征值及特征向量计算 程序运行结果：程序运行结果： MatrixMatrix 2.000000 3.000000 2.000000 2.000000 3.000000 2.000000 10.000000 3.000000 4.000000 10.000000 3.000000 4.000000 3.000000

23、6.000000 1.000000 3.000000 6.000000 1.000000 Min EigenValueMin EigenValue 0.500014 0.500014 Min EigenVectorMin EigenVector 0.200012 0.200012 0.399985 0.399985 -1.000000 -1.000000 第6章矩阵特征值及特征向量计算 6.3 6.3 雅可比（雅可比（Jacobi）法法 6.3.1 6.3.1 雅可比法的基本思想雅可比法的基本思想 雅可比法是求实对称矩阵雅可比法是求实对称矩阵 的特征值和的特征值和 特征向量的常用方法。其基本思

24、想是通过一特征向量的常用方法。其基本思想是通过一 组平面旋转变换（正交相似变换），将实对组平面旋转变换（正交相似变换），将实对 称矩阵称矩阵 化为对角阵，从而得到其全部特征化为对角阵，从而得到其全部特征 值和相应的特征向量。值和相应的特征向量。 A A 第6章矩阵特征值及特征向量计算 6.3.2 6.3.2 实现雅可比法的基本步骤实现雅可比法的基本步骤 雅可比法求雅可比法求 阶实对称矩阵阶实对称矩阵 的全部特征值和相应的特征向量的的全部特征值和相应的特征向量的 基本步骤如下基本步骤如下: : （1 1） 输入矩阵的阶数输入矩阵的阶数 、矩阵、矩阵 、预先给定的精度、预先给定的精度 和最大迭代次

25、数和最大迭代次数 KMAX KMAX； （2 2） , , 为为 阶单位矩阵；阶单位矩阵； （3 3）在）在 中选取非对角线元素中绝对值最大的元素中选取非对角线元素中绝对值最大的元素 ； （4 4）若）若 ，则迭代结束。此时，对角线元素，则迭代结束。此时，对角线元素 即为特征值即为特征值 ，矩阵，矩阵 的第的第 列为与列为与 对应的特征向量；否则对应的特征向量；否则 继续做下一步；继续做下一步； nA n IV I n pq a pq a) 121(n ,iaij V i i A A i 第6章矩阵特征值及特征向量计算 pq ax )( 2 1 ppqq aay ),2, 1,(nkkji 2

26、2 )sign( yx x yw wsin2 )-1(12 sin 2 w w 2 sin1cos cos2sincos 221 qpqqpppp aaaa sin2cossin 221 pqqqppqq aaaa （5 5） 按下列公式计算矩阵按下列公式计算矩阵 的新元素：的新元素： 第6章矩阵特征值及特征向量计算 0 11 qppq aa qp,j,aaa jqjpjp sincos 1 qp,j,aaa jqjpjq cossin 1 qp,i,aaa iqipip sincos 1 qp,i,aaa iqipip cossin 1 qp,ji,aa ijij 1 第6章矩阵特征值及特征

27、向量计算 10 cossin sincos 01 q p q,p,R qp （6 6） ，转第（，转第（2 2）步计算，其中）步计算，其中 为如下所示的正交矩阵。为如下所示的正交矩阵。 )(q,p,RVV )(,q,pR 第6章矩阵特征值及特征向量计算 例例6.36.3 使用雅可比法计算下列实对称矩 阵的特征值和相应的一组特征向量，特 征值精度为 。 210 121 012 A 6 10 第6章矩阵特征值及特征向量计算 # #include include #include #include #define N 3 /#define N 3 /* * N N为矩阵的阶数为矩阵的阶数 * */

28、/ # #define eps 1e-06 /define eps 1e-06 /* * 给定精度要求给定精度要求 * */ / # #define KMAX 60 /define KMAX 60 /* * 最大迭代次数最大迭代次数 * */ / 第6章矩阵特征值及特征向量计算 int jacobieigen(double aN,double uN)int jacobieigen(double aN,double uN) int i,j,k; int i,j,k; int p,q; int p,q; double b,d,m,x,y; double b,d,m,x,y; double sn,c

29、n,w; double sn,cn,w; for(i=0;iN;i+) / for(i=0;iN;i+) /* * 产生单位矩阵产生单位矩阵 * */ / uii=1.0;uii=1.0; for(j=0;jN;j+) for(j=0;jN;j+) if(i!=j) if(i!=j) uij=0.0; uij=0.0; k=1; k=1; while(1) while(1) m=0.0; m=0.0; for(i=1;i=N-1;i+) / for(i=1;i=N-1;i+) /* * 选取绝对值最大的对角线元素选取绝对值最大的对角线元素 * */ / for(j=0;j=i-1;j+)for

30、(j=0;jm) if(i!=j) p=i; q=j; m=d; p=i; q=j; if(meps) / if(mKMAX) / if(kKMAX) /* * 超过最大迭代次数返回超过最大迭代次数返回 * */ / return(-1);return(-1); k=k+1; k=k+1; x=-apq; x=-apq; y=(aqq-app)/2.0; y=(aqq-app)/2.0; w=x/sqrt(x w=x/sqrt(x* *x+yx+y* *y);y); if(y0.0) w=-w; if(y0.0) w=-w; sn=1+sqrt(1.0-w sn=1+sqrt(1.0-w* *

31、w);w); sn=w/sqrt(2.0 sn=w/sqrt(2.0* *sn);sn); cn=sqrt(1.0-sn cn=sqrt(1.0-sn* *sn);sn); m=app; / m=app; /* * 计算矩阵的新元素计算矩阵的新元素 * */ / 第6章矩阵特征值及特征向量计算 app=m app=m* *cncn* *cn+aqqcn+aqq* *snsn* *sn+apqsn+apq* *w;w; aqq=m aqq=m* *snsn* *sn+aqqsn+aqq* *cncn* *cn-apqcn-apq* *w;w; apq=0.0; apq=0.0; aqp=0.0;

32、 aqp=0.0; for(j=0;jN;j+) for(j=0;jN;j+) if(j!=p) apj=m m=apj;apj=m* *cn+aqjcn+aqj* *sn; aqj=-msn; aqj=-m* *sn+aqjsn+aqj* *cn; cn; for(i=0;iN;i+) for(i=0;iN;i+) if(i!=p) m=aip; aip=m aip=m* *cn+aiqcn+aiq* *sn;sn; aiq=-m aiq=-m* *sn+aiqsn+aiq* *cn; cn; for(i=0;iN;i+) for(i=0;iN;i+) m=uip; m=uip; uip=m

33、 uip=m* *cn+uiqcn+uiq* *sn;sn; uiq=-m uiq=-m* *sn+uiqsn+uiq* *cn;cn; 第6章矩阵特征值及特征向量计算 main()main() double aNN=2.0,-1.0,0.0,-1.0,2.0,-1.0,0.0,-1.0,2.0; double aNN=2.0,-1.0,0.0,-1.0,2.0,-1.0,0.0,-1.0,2.0; double bNN,vNN; double bNN,vNN; int i,j,k; int i,j,k; printf(Matrixn); printf(Matrixn); for(i=0;iN

34、;i+) for(i=0;iN;i+) printf( ); printf( ); for(j=0;jN;j+) printf(%11f,aij); for(j=0;jN;j+) printf(%11f,aij); printf(n); printf(n); k=jacobieigen(a,v); k=jacobieigen(a,v); if(k=1) if(k=1) printf(nEigenValuen); printf(nEigenValuen); for(i=0;iN;i+) for(i=0;iN;i+) printf( %d:,i+1); printf( %d:,i+1); for(

35、j=0;jN;j+) for(j=0;jN;j+) if(i=j) printf(%11f,aij); if(i=j) printf(%11f,aij); printf(n); printf(n); printf(n); printf(n); 第6章矩阵特征值及特征向量计算 for(i=0;iN;i+) for(i=0;iN;i+) for(j=0;jN;j+) for(j=0;jN;j+) bij=vji; bij=vji; printf(EigenVectorn); printf(EigenVectorn); for(i=0;iN;i+) for(i=0;iN;i+) printf( %d

36、:,i+1); printf( %d:,i+1); for(j=0;jN;j+) for(j=0;jN;j+) printf(%11f,bij); printf(%11f,bij); printf(n); printf(n); else else printf(no EigenValue); printf(no EigenValue); 第6章矩阵特征值及特征向量计算 程序运行结果：程序运行结果： MatrixMatrix 2.000000 -1.000000 0.0000002.000000 -1.000000 0.000000 -1.000000 2.000000 -1.000000 -1

37、.000000 2.000000 -1.000000 0.000000 -1.000000 2.000000 0.000000 -1.000000 2.000000 EigenValueEigenValue 1 1： 3.414214 3.414214 2: 0.585786 2: 0.585786 3: 2.000000 3: 2.000000 EigenVectorEigenVector 1: 0.500000 -0.707107 0.500000 1: 0.500000 -0.707107 0.500000 2: 0.500000 0.707107 0.500000 2: 0.50000

38、0 0.707107 0.500000 3: -0.707107 0.000000 0.707107 3: -0.707107 0.000000 0.707107 第6章矩阵特征值及特征向量计算 6.4 6.4 雅可比过关法雅可比过关法 6.4.1 6.4.1 雅可比过关法的基本思想雅可比过关法的基本思想 在雅可比法中，每进行一次平面旋转变换 前都需要在非对角的元素中选取绝对值最大的 元素，这是很费时间的。雅可比过关法对此进 行了改进，具体做法如下。 （1）输入矩阵的阶数 、矩阵 以及容许误 差 ； （2）计算 阶实对称矩阵 的所有非对角线元 素平方之和的平方根： n A nA 1 11 2

39、0 2 n i n ij ij as 第6章矩阵特征值及特征向量计算 （3 3）设置关口）设置关口 , ,在非对角线元素中逐行扫在非对角线元素中逐行扫 描选取第一个绝对值大于或等于描选取第一个绝对值大于或等于 的元素的元素 进行进行 平面旋转变换，直到所有非对角线元素的绝对值都平面旋转变换，直到所有非对角线元素的绝对值都 小于小于 为止；为止； （4 4）再设置关口）再设置关口 ，重复以上过程；，重复以上过程； （5 5）以此类推，经过一系列的关口）以此类推，经过一系列的关口 ，直，直 到对于某一个关口满足条件到对于某一个关口满足条件 为止。为止。 nss 01 1 s pq a 1 s 2

40、012 nsnss ,s ,s ,s 321 k s 第6章矩阵特征值及特征向量计算 6.4.2 6.4.2 实现雅可比过关法的基本步骤实现雅可比过关法的基本步骤 雅可比过关法求雅可比过关法求 阶实对称矩阵阶实对称矩阵 的的 特征值和相应的特征向量的基本步骤，与特征值和相应的特征向量的基本步骤，与 雅可比法基本相同，只是在选取非对角线雅可比法基本相同，只是在选取非对角线 上的绝对值最大元素时采用了前述方法，上的绝对值最大元素时采用了前述方法， 因此，在此不再重复。因此，在此不再重复。 A 第6章矩阵特征值及特征向量计算 例例6.4 6.4 使用雅可比过关法计算下列实对称矩阵的特征值和使用雅可比

41、过关法计算下列实对称矩阵的特征值和 相应的一组特征向量，特征值精度为相应的一组特征向量，特征值精度为 。 #include #include #define N 3 /* N为矩阵的阶数 */ #define eps 1e-06 /* 给定精度要求 */ 151534 112323 53712 32191 432110 A 6 10 第6章矩阵特征值及特征向量计算 void jacobieigeng(double aN,double uN)void jacobieigeng(double aN,double uN) int i,j,k; int i,j,k; int p,q; int p,q;

42、 double b,d,m,x,y; double b,d,m,x,y; double sn,cn,s,w; double sn,cn,s,w; for(i=0;iN;i+) / for(i=0;iN;i+) /* * 产生单位矩阵产生单位矩阵 * */ / uii=1.0;uii=1.0; for(j=0;jN;j+) for(j=0;jN;j+) if(i!=j) if(i!=j) uij=0.0; uij=0.0; s=0.0; s=0.0; for(i=1;i=N-1;i+) for(i=1;i=N-1;i+) for(j=0;j=i-1;j+) for(j=0;j=i-1;j+) d

43、=fabs(aij); d=fabs(aij); s=s+d s=s+d* *d;d; 第6章矩阵特征值及特征向量计算 s=sqrt(2.0 s=sqrt(2.0* *s);s); loap0: / loap0: /* * 设置关口设置关口 * */ / s=s/(1.0s=s/(1.0* *(N-1);(N-1); loap1: loap1: for(i=1;i=N-1;i+) / for(i=1;i=N-1;i+) /* * 选取对角线元素中的绝对值最大的元素选取对角线元素中的绝对值最大的元素* */ / for(j=0;j=i-1;j+)for(j=0;js) / if(ds) /* *

44、 未过关口，做旋转变换未过关口，做旋转变换 * */ / p=i; q=j; goto loap2; p=i; q=j; goto loap2; if(seps) / if(seps) /* * 满足精度要求，返回满足精度要求，返回 * */ / return;return; goto loap0; / goto loap0; /* * 设置下一道关口设置下一道关口 * */ / loap2:loap2: x=-apq; x=-apq; y=(aqq-app)/2.0; y=(aqq-app)/2.0; w=x/sqrt(x w=x/sqrt(x* *x+yx+y* *y);y); if(y0

45、.0) if(y0.0) w=-w; w=-w; 第6章矩阵特征值及特征向量计算 sn=1+sqrt(1.0-w sn=1+sqrt(1.0-w* *w);w); sn=w/sqrt(2.0 sn=w/sqrt(2.0* *sn);sn); cn=sqrt(1.0-sn cn=sqrt(1.0-sn* *sn);sn); m=app; / m=app; /* * 计算矩阵计算矩阵A A的新元素的新元素 * */ / app=mapp=m* *cncn* *cn+aqqcn+aqq* *snsn* *sn+apqsn+apq* *w;w; aqq=m aqq=m* *snsn* *sn+aqqs

46、n+aqq* *cncn* *cn-apqcn-apq* *w;w; apq=0.0; apq=0.0; aqp=0.0; aqp=0.0; for(j=0;jN;j+) for(j=0;jN;j+) if(j!=p) m=apj; apj=m apj=m* *cn+aqjcn+aqj* *sn;sn; aqj=-m aqj=-m* *sn+aqjsn+aqj* *cn;cn; for(i=0;iN;i+) for(i=0;iN;i+) if(i!=p) m=aip; aip=m aip=m* *cn+aiqcn+aiq* *sn;sn; aiq=-m aiq=-m* *sn+aiqsn+ai

47、q* *cn;cn; 第6章矩阵特征值及特征向量计算 for(i=0;iN;i+) for(i=0;iN;i+) m=uip; m=uip; uip=m uip=m* *cn+uiqcn+uiq* *sn;sn; uiq=-m uiq=-m* *sn+uiqsn+uiq* *cn;cn; goto loap1; / goto loap1; /* *继续逐行扫描继续逐行扫描, ,让每一个非对角线元素过关让每一个非对角线元素过关* */ / 第6章矩阵特征值及特征向量计算 main()main() static double aN=10.0,1.0,2.0,3.0,4.0, static doub

48、le aN=10.0,1.0,2.0,3.0,4.0, 1.0,9.0,-1.0,2.0,-3.0, 1.0,9.0,-1.0,2.0,-3.0, 2.0,-1.0,7.0,3.0,-5.0, 2.0,-1.0,7.0,3.0,-5.0, 3.0,2.0,3.0,12.0,-1.0, 3.0,2.0,3.0,12.0,-1.0, 4,-3,-5,-1,15; 4,-3,-5,-1,15; double bNN,vNN; double bNN,vNN; int i,j; int i,j; clrscr(); clrscr(); printf(Matrixn); printf(Matrixn);

49、for(i=0;iN;i+) for(i=0;iN;i+) printf( ); printf( ); for(j=0;jN;j+) for(j=0;jN;j+) printf(%11f,aij); printf(%11f,aij); printf(n); printf(n); 第6章矩阵特征值及特征向量计算 jacobieigeng(a,v); jacobieigeng(a,v); printf(nEigenValuen); printf(nEigenValuen); for(i=0;iN;i+) for(i=0;iN;i+) printf( %d:,i+1); printf( %d:,i+

50、1); for(j=0;jN;j+) for(j=0;jN;j+) if(i=j) printf(%11f,aij); if(i=j) printf(%11f,aij); printf(n); printf(n); printf(n); printf(n); for(i=0;iN;i+) for(i=0;iN;i+) for(j=0;jN;j+) bij=vji; for(j=0;jN;j+) bij=vji; printf(EigenVectorn); printf(EigenVectorn); for(i=0;iN;i+) for(i=0;iN;i+) printf( %d:,i+1);

51、printf( %d:,i+1); for(j=0;jN;j+) printf(%11f,bij); for(j=0;jN;j+) printf(%11f,bij); printf(n); printf(n); 第6章矩阵特征值及特征向量计算 程序运行结果：程序运行结果： MatrixMatrix 10.000000 1.000000 2.000000 3.000000 4.000000 10.000000 1.000000 2.000000 3.000000 4.000000 1.000000 9.000000 -1.000000 2.000000 -3.000000 1.000000 9.000000 -1.000000 2.000000 -3.000000 2.000000 -1.000000 7.000000 3.000000 -5.000000 2.000000 -1.000000 7.000000 3.000000 -5.00000