基于雅可比方法的特征值特征向量求解函数（原创）

hyl2323

1. #include<iostream.h>
   
2. #include<math.h>
   
3. #include "Matrix.h"
   
4. void Eig(int n,Matrix a,Matrix &v,Matrix &r)
   
5. //n is the dimension of square matrix a;
   
6. //a is the matrix be dealed;
   
7. //v is the arrray for returning eigenvalues
   
8. //r is the matrix for returning eigenvectors(columns)
   
9. {
   
10. //preparing
    
11. int i,j;
    
12. Matrix b(n,n);
    
13. Matrix s(n,n);
    
14. //yacobi method
    
15. int p,q,I=1;
    
16. double v0=0,n1,eps=1e-6,n2,c,tt,sine,cosine,temp=0;
    
17. for(i=0;i<n;i++){
    
18. for(j=0;j<n;j++){
    
19.   if(i==j)
    
20.       r(i,j)=1;
    
21.   else
    
22.    r(i,j)=0;
    
23. }
    
24. }
    
25. for(i=0;i<n;i++){
    
26. for(j=0;j<n;j++){
    
27.   if(i!=j)
    
28.    temp+=a(i,j)*a(i,j);
    
29. }
    
30. }
    
31. n1=sqrt(temp);
    
32. n2=(eps/n)*n1;
    
33. n1=n1/n;
    
34. while(n1>n2){
    
35. while(I){
    
36.         I=0;
    
37.         for(i=0;i<n-1;i++){
    
38.             for(j=i+1;j<n;j++){
    
39.           if(fabs(a(i,j))>=n1){
    
40.            I=1;
    
41.      c=(a(i,i)-a(j,j))/(2*a(i,j));
    
42.      if(c>=0)
    
43.       tt=1/(fabs(c)+sqrt(c*c+1));
    
44.      else
    
45.       tt=(-1)/(fabs(c)+sqrt(c*c+1));   
    
46.      cosine=1/sqrt(1+tt*tt);
    
47.      sine=tt/sqrt(1+tt*tt);
    
48.      for(p=0;p<n;p++){
    
49.       for(q=0;q<n;q++){
    
50.        b(p,q)=a(p,q);
    
51.        s(p,q)=r(p,q);
    
52.       }
    
53.      }
    
54.      for(p=0;p<n;p++){
    
55.       for(q=0;q<n;q++){
    
56.        if(q==i && p!=i && p!=j)
    
57.         b(p,i)=b(i,p)=a(p,i)*cosine+a(p,j)*sine;
    
58.        if(q==j && p!=i && p!=j)
    
59.         b(p,j)=b(j,p)=a(p,j)*cosine-a(p,i)*sine;
    
60.        if(q==i)
    
61.         s(p,i)=r(p,i)*cosine+r(p,j)*sine;
    
62.        if(q==j)
    
63.         s(p,j)=r(p,j)*cosine-r(p,i)*sine;
    
64.       }
    
65.      }
    
66.      b(i,i)=a(i,i)*cosine*cosine+a(j,j)*sine*sine+2*a(i,j)*sine*cosine;
    
67.      b(j,j)=a(i,i)*sine*sine+a(j,j)*cosine*cosine-2*a(i,j)*sine*cosine;
    
68.      b(i,j)=b(j,i)=0;
    
69.      for(p=0;p<n;p++){
    
70.       for(q=0;q<n;q++){
    
71.        a(p,q)=b(p,q);
    
72.        r(p,q)=s(p,q);
    
73.       }
    
74.      }
    
75.     }
    
76.    }
    
77.   }
    
78. }
    
79.     n1/=n;
    
80. I=1;
    
81. }
    
82. //get eigenvalues from matrix a and set them to vector v as the output variable
    
83. for(i=0;i<n;i++){
    
84. for(j=0;j<n;j++){
    
85.   if(i==j)
    
86.       v(i,0)=a(i,j);
    
87. }
    
88. }
    
89. }

复制代码

具体算法是参考浙大版《数值分析》，#include "Matrix.h" 是包含自定义的矩阵类，与matlab对比调试无误！

[ 本帖最后由 风花雪月 于 2007-12-24 11:07 编辑 ]

hyl2323

本人所编程序中过程最痛苦的一个！！！花了一个星期。谨以此纪念我死去的脑细胞。

xiaowanzimm

谢谢共享

zz80

楼主好人啊

wanyeqing2003

雅可比方法是比较成熟的方法，适用于自由度比较少的系统。

楼主把程序调试出来是需要费些幸苦。

mulan

原帖由 hyl2323 于 2007-12-20 19:28 发表

                               
登录/注册后可看大图

本人所编程序中过程最痛苦的一个！！！花了一个星期。谨以此纪念我死去的脑细胞。

我想请问一下，有些矩阵用matlab算出来的特征值是复数，这个程序能表示出来吗？谢谢

风花雪月

原帖由 mulan 于 2008-2-26 17:02 发表

                               
登录/注册后可看大图

我想请问一下，有些矩阵用matlab算出来的特征值是复数，这个程序能表示出来吗？谢谢

这个程序不是所有矩阵都能计算的，建议楼主注明适用条件

snowswallowhe

楼主能把Matrix.h的代码发上来吗??谢谢

风花雪月

原帖由 snowswallowhe 于 2008-4-9 17:13 发表

                               
登录/注册后可看大图

楼主能把Matrix.h的代码发上来吗??谢谢

http://ggt.sourceforge.net/html/Matrix_8h-source.html
看看是不是这个，或者直接和搂主短消息联系

lq12131010

一堆垃圾 ..雅可比几乎没有实际价值..

zyhhy_x

这个,为什么没主函数的,请指教