100分 用雅可比求实对称矩阵的特征值和特征向量

[shenxy0596]

public int jcbi(Double[,] A, int n, Double[,] v, Double eps,int jt)//A存放n阶实对称矩阵,返回时对角线上存放n个特征值，n实对称矩阵A的阶数,V[n,n]返回特征向量，eps控制精度要求，jt控制最大迭代次数， int jcbi()方法返回标志值
{
int i, j, p, q, u, w, t, s, l;
p = 0; q = 0;
Double fm, cn, sn, omega, x, y, d;
l = 1;
for (i = 0; i <= n - 1; i++)
{
v[i,i] = 1.0;
for (j = 0; j <= n - 1; j++)
{
if (i != j) v[i ,j] = 0.0;
}
}
while (l == 1)
{
fm = 0.0;
for (i = 1; i <= n - 1; i++)
{
for (j = 0; j <= i - 1; j++)
{
d = Math.Abs(A[i, j]);
if ((i != j) && (d > fm))
{
fm = d; p = i; q = j;
}

}
}
if (fm < eps) return 1;
if (l > jt) return -1;
l = l + 1;
u = p * n + q;
w = p * n + p;
t = q * n + p;
s = q * n + q;
x =0-A[p,q];
y = (A[q,q] - A[p,p])/2.0;
omega = x / Math.Sqrt(x * x + y * y);
if (y < 0.0) omega = 0-omega;
sn = 1.0 + Math.Sqrt(1.0 - omega * omega);
sn = omega / Math.Sqrt(2.0 * sn);
cn = Math.Sqrt(1.0 - sn * sn);
fm = A[p,p];
A[p,p] = fm * cn * cn + A[q,q] * sn * sn + A[p,q] * omega;
A[q,q] = fm * sn * sn + A[q,q] * cn * cn - A[p,q] * omega;
A[p,q] = 0.0;
A[q,p] = 0.0;
for (j = 0; j <= n - 1; j++)
{
if ((j != p) && (j != q))
{
u = p * n + j;
w = q * n + j;
fm = A[p,j];
A[p,j] = fm * cn + A[q,j] * sn;
A[q,j] = -fm * sn + A[q,j] * cn;
}
}
for (i = 0; i <= n - 1; i++)
{
if ((i != p) && (i != q))
{
u = i * n + p;
w = i * n + q;
fm = A[i,p];
A[i,p] = fm * cn + A[i,q] * sn;
A[i,q] = -fm * sn + A[i,q] * cn;
}
}
for(i=0;i<=n-1;i++)
{
u=i*n+p;
w=i*n+q;
fm=v[i,p];
v[i,p]=fm*cn+v[i,q]*sn;
v[i,q]=-fm*sn+v[i,q]*cn;
}



}
return 1;

}

这是从徐士良书上改编成C#的，不知道哪里有错，哪位高手帮我看一下

[shenxy0596]

不知道徐士良那本书上的代码是否有误
或谁有用C写的代码也可以，拜托!谢谢

[dai_weitao]

书上代码应该没错误吧
C#不懂, 帮你顶了

[spofmy]

只能顶了

[星羽]

同上

[jixingzhong]

用C语言来实现的，雅可法可用于法实对称矩阵的全部特征值与对应的特征值向量：
http://www.programsalon.com/downloads31/sourcecode/math/detail100765.html

[jixingzhong]

http://www.dearbook.com.cn/book/viewbook.aspx?pno=TS00124175

C#数值计算算法编程
貌似这本书很适合啊，
楼主去 图书大厦 翻出来看看吧：
第3章 矩阵运算
　3.1 矩阵类设计
　3.2 矩阵基础运算
　3.3 实矩阵求逆的全选主元高斯－约当法
　3.4 复矩阵求逆的全球主元高斯－约当法
　3.5 对称正定矩阵的求逆
　3.6 托利兹矩阵求逆的埃兰特方法
　3.7 求行列式值的全选主元高斯消去法
　3.8 求矩阵秩的全选主元高斯消去法
　3.9 对称正定矩阵的乔里斯基分解与行列式的求值
　3.10 短阵的三角三解
　3.11 一般实矩阵的QR分解
　3.12 一般实矩阵的奇异值分解
　3.13 求广义逆的奇异值分解法
　3.14 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
　3.15 实对称三对角阵的全部特征值与特征向量的计算
　3.16 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
　3.17 求赫申伯格矩阵全部特征值的QR方法
　3.18 求实对称矩阵特征值与特征向量的雅可比法
　3.19 求实对称矩阵特征值与特征向量的雅可比过关法

[jixingzhong]

常用算法程序集，
里面有 求实对称矩阵特征值与特征向量的雅可比法，
应该在第五章， 对应 CH5 文件夹下，
具体那个文件就自己看看吧。

C语言常用算法源代码下载：
http://download.csdn.net/source/188851
http://download.chinaitlab.com/soft/7925.htm

[shenxy0596]

我去看看

[shenxy0596]

那本书我也有，可是在哪章找不到呀

[shenxy0596]

这段代码就是从<<常用算法程序集>>里面改编的

[bobnie]

//求实对称矩阵的特征值及特征向量的雅格比法
//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
//返回值小于0表示超过迭代jt次仍未达到精度要求
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值
//n-矩阵的阶数
//u-长度为n*n的数组，返回特征向量(按列存储)
//eps-控制精度要求
//jt-整型变量，控制最大迭代次数
int eejcb(double a[],int n,double v[],double eps,int jt)
{
int i,j,p,q,u,w,t,s,l;
double fm,cn,sn,omega,x,y,d;
l=1;
for (i=0; i<=n-1; i++)
{
v[i*n+i]=1.0;
for (j=0; j<=n-1; j++)
{
if (i!=j)
{
v[i*n+j]=0.0;
}
}
}
while (1==1)
{
fm=0.0;
for (i=0; i<=n-1; i++)
{
for (j=0; j<=n-1; j++)
{
d=fabs(a[i*n+j]);
if ((i!=j)&&(d>fm))
{
fm=d;
p=i;
q=j;
}
}
}
if (fm<eps)
{
return(1);
}
if (l>jt)
{
return(-1);
}
l=l+1;
u=p*n+q;
w=p*n+p;
t=q*n+p;
s=q*n+q;
x=-a[u];
y=(a[s]-a[w])/2.0;
omega=x/sqrt(x*x+y*y);
if (y<0.0)
{
omega=-omega;
}
sn=1.0+sqrt(1.0-omega*omega);
sn=omega/sqrt(2.0*sn);
cn=sqrt(1.0-sn*sn);
fm=a[w];
a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;
a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;
a[u]=0.0;
a[t]=0.0;
for (j=0; j<=n-1; j++)
{
if ((j!=p)&&(j!=q))
{
u=p*n+j;
w=q*n+j;
fm=a[u];
a[u]=fm*cn+a[w]*sn;
a[w]=-fm*sn+a[w]*cn;
}
}
for (i=0; i<=n-1; i++)
{
if ((i!=p)&&(i!=q))
{
u=i*n+p;
w=i*n+q;
fm=a[u];
a[u]=fm*cn+a[w]*sn;
a[w]=-fm*sn+a[w]*cn;
}
}
for (i=0; i<=n-1; i++)
{
u=i*n+p;
w=i*n+q;
fm=v[u];
v[u]=fm*cn+v[w]*sn;
v[w]=-fm*sn+v[w]*cn;
}
}
return(1);
}

[krocwang]

mark

[青葱岁月不可追]

来学习

[wudeshou82666]

[Quote=引用 13 楼 jillnicky 的回复:]
来学习
[/Quote]
来学习

[miracle159]

mark

[super_hover]

学习学习

[jachanyellow]

学习了

[老王爱上猫]

又见挖坟贴

[洋桔子]

运行结果不对呀，eg3*3的矩阵：-2，1，1，0，2，0，-4，1，3