基于雅可比方法的特征值特征向量求解函数(原创)
#include<iostream.h>#include<math.h>
#include "Matrix.h"
void Eig(int n,Matrix a,Matrix &v,Matrix &r)
//n is the dimension of square matrix a;
//a is the matrix be dealed;
//v is the arrray for returning eigenvalues
//r is the matrix for returning eigenvectors(columns)
{
//preparing
int i,j;
Matrix b(n,n);
Matrix s(n,n);
//yacobi method
int p,q,I=1;
double v0=0,n1,eps=1e-6,n2,c,tt,sine,cosine,temp=0;
for(i=0;i<n;i++){
for(j=0;j<n;j++){
if(i==j)
r(i,j)=1;
else
r(i,j)=0;
}
}
for(i=0;i<n;i++){
for(j=0;j<n;j++){
if(i!=j)
temp+=a(i,j)*a(i,j);
}
}
n1=sqrt(temp);
n2=(eps/n)*n1;
n1=n1/n;
while(n1>n2){
while(I){
I=0;
for(i=0;i<n-1;i++){
for(j=i+1;j<n;j++){
if(fabs(a(i,j))>=n1){
I=1;
c=(a(i,i)-a(j,j))/(2*a(i,j));
if(c>=0)
tt=1/(fabs(c)+sqrt(c*c+1));
else
tt=(-1)/(fabs(c)+sqrt(c*c+1));
cosine=1/sqrt(1+tt*tt);
sine=tt/sqrt(1+tt*tt);
for(p=0;p<n;p++){
for(q=0;q<n;q++){
b(p,q)=a(p,q);
s(p,q)=r(p,q);
}
}
for(p=0;p<n;p++){
for(q=0;q<n;q++){
if(q==i && p!=i && p!=j)
b(p,i)=b(i,p)=a(p,i)*cosine+a(p,j)*sine;
if(q==j && p!=i && p!=j)
b(p,j)=b(j,p)=a(p,j)*cosine-a(p,i)*sine;
if(q==i)
s(p,i)=r(p,i)*cosine+r(p,j)*sine;
if(q==j)
s(p,j)=r(p,j)*cosine-r(p,i)*sine;
}
}
b(i,i)=a(i,i)*cosine*cosine+a(j,j)*sine*sine+2*a(i,j)*sine*cosine;
b(j,j)=a(i,i)*sine*sine+a(j,j)*cosine*cosine-2*a(i,j)*sine*cosine;
b(i,j)=b(j,i)=0;
for(p=0;p<n;p++){
for(q=0;q<n;q++){
a(p,q)=b(p,q);
r(p,q)=s(p,q);
}
}
}
}
}
}
n1/=n;
I=1;
}
//get eigenvalues from matrix a and set them to vector v as the output variable
for(i=0;i<n;i++){
for(j=0;j<n;j++){
if(i==j)
v(i,0)=a(i,j);
}
}
}
具体算法是参考浙大版《数值分析》,#include "Matrix.h" 是包含自定义的矩阵类,与matlab对比调试无误!
[ 本帖最后由 风花雪月 于 2007-12-24 11:07 编辑 ] 本人所编程序中过程最痛苦的一个!!!花了一个星期。谨以此纪念我死去的脑细胞。 谢谢共享 楼主好人啊 雅可比方法是比较成熟的方法,适用于自由度比较少的系统。
楼主把程序调试出来是需要费些幸苦。 原帖由 hyl2323 于 2007-12-20 19:28 发表 http://www.chinavib.com/forum/images/common/back.gif
本人所编程序中过程最痛苦的一个!!!花了一个星期。谨以此纪念我死去的脑细胞。
我想请问一下,有些矩阵用matlab算出来的特征值是复数,这个程序能表示出来吗?谢谢 原帖由 mulan 于 2008-2-26 17:02 发表 http://www.chinavib.com/forum/images/common/back.gif
我想请问一下,有些矩阵用matlab算出来的特征值是复数,这个程序能表示出来吗?谢谢
这个程序不是所有矩阵都能计算的,建议楼主注明适用条件 楼主能把Matrix.h的代码发上来吗??谢谢 原帖由 snowswallowhe 于 2008-4-9 17:13 发表 http://www.chinavib.com/forum/images/common/back.gif
楼主能把Matrix.h的代码发上来吗??谢谢
http://ggt.sourceforge.net/html/Matrix_8h-source.html
看看是不是这个,或者直接和搂主短消息联系 一堆垃圾 ..雅可比几乎没有实际价值.. 这个,为什么没主函数的,请指教
页:
[1]