主题：[讨论]关于矩阵特征向量（fortran与matlab对照）
tianhy2010
[专家分：60] 发布于 2010-08-31 11:06:00
最近一直很疑惑，为什么对于复数矩阵，得到的特征向量会有差距呢？尽管从数学上讲，对于一个矩阵，特征向量有无穷多个，但是yeg001说，matlab是用lapack求解矩阵的, IMSL也是用lapack求解矩阵的. 也就是说二者应该得到相同的特征向量，为了检查到底哪里出了问题，我做了以下操作：
1。matlab和fortran解同一个复数矩阵c(4,4)，程序如下：
matlab程序:
clear;
clc;
format long;
t=1;
f=1/4;
w=0.0;
n=3;
C=[w,t*exp(-i*(-pi+2*pi/3*f*(1.5*n+1/4))),0,0;t*exp(i*(-pi+2*pi/3*f*(1.5*n+1/4))),w,t,0;0,t,w,t*exp(i*(-pi+2*pi/3*f*(1.5*(n+1)+1/4)));0,0,t*exp(-i*(-pi+2*pi/3*f*(1.5*(n+1)+1/4))),w];
vvvv=eig(C)
[u,v]=eig(C);
fortran程序：
program main
use IMSL
implicit none
real*8, parameter::pi=3.14159265
complex*16, parameter:: fi=(0.0, 1.0)
real*8 t,f,w,temp
complex*16 ::c(4,4)
complex*16::eigenvector(4,4)
real*8:: eigenvalue(4) !特征值的矩阵
integer ::i,n,j,p
w=0.0
t=1.0
f=0.25
call hmatr(c)
!求矩阵特征值和特征向量
eigenvalue=eig(c,w=eigenvector)
open(1,file='dmatr.txt')
do i=1,4
do j=1,4
write(1,*) c(j,i)
end do
end do
close(1)
!特征值排序
do j=1,4-1
p=j
do i=j+1,4
if(eigenvalue(i)<eigenvalue(p)) p=i
end do
temp=eigenvalue(j)
eigenvalue(j)=eigenvalue(p)
eigenvalue(p)=temp
end do
!特征值求值，写入文件
open(10,file='eeig.txt')
do i=1,4
write(10,*) eigenvalue(i)
end do
close(10)
!将特征向量写入文件
open(20,file='jing2.txt')
do i=1,4
do j=1,4
write(20,*) eigenvector(j,i)
end do
end do
close(20)
stop
contains
subroutine hmatr(c)
integer n0
complex*16:: c(4,4)
c=0.0
n0=3
c(1,1)=w
c(2, 1) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(1, 2) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(2, 2) = w
c(2, 3) = t
c(3, 2) = t
c(3, 3) = w
c(3, 4) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 3) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 4) = w
end subroutine hmatr
end program main
两个程序我都求了，理论上讲应该是相同的，但是不知道为什么，得到的矩阵相同，得到的特征值相同，但是特征向量还是不同啊。是因为复数矩阵的原因吗？
回复列表 （共32个回复）
沙发
tianhy2010 [专家分：60] 发布于 2010-08-31 11:21:00
2种编辑器得到的特征向量不同是有可能的，只要它们满足Hr=Ar,H--矩阵，A--矩阵特征值，就说明r是矩阵H的特征向量。
看下它们得到的结果怎样吧：
matlab:几乎完全相同：
C0*u(:,4)
ans =
0.425325404176020 + 0.425325404176020i
0.964922709445035 + 0.127034484677236i
0.964922709445035 + 0.127034484677236i
0.601500955007546 + 0.000000000000000i
>> v(4,4)*u(:,4)
ans =
0.425325404176020 + 0.425325404176020i
0.964922709445035 + 0.127034484677236i
0.964922709445035 + 0.127034484677236i
0.601500955007545
再看下fortran得到的结果吧：在程序中加入语句：
！记录矩阵与特征向量的乘积
open(12,file='ceshi1.txt')
write(12,*) matmul(c,eigenvector(:,4))
close(12)
！记录特征值与特征向量的乘积
open(13,file='ceshi2.txt')
write(13,*) eigenvalue(4)*eigenvector(:,4)
close(13)
得到的结果大跌眼镜啊，是不是我的语句有问题？
板凳
tianhy2010 [专家分：60] 发布于 2010-08-31 21:16:00
大家好，你们感觉哪里不明白？fortran程序太长吗？看完有什么感想说下哈，水帖也要
3 楼
aliouying [专家分：1150] 发布于 2010-08-31 23:20:00
嗯，我看看~ 纯属水贴~
4 楼
aliouying [专家分：1150] 发布于 2010-08-31 23:38:00
IMSL求Hermitian matrix的函数是EVCHF~
5 楼
tianhy2010 [专家分：60] 发布于 2010-09-01 09:28:00
[quote]IMSL求Hermitian matrix的函数是EVCHF~[/quote]
[em21][em21][em21][em21][em21][em21][em21][em21]
能否讲的具体点啊？这个函数怎么百度搜不出来呢？
6 楼
yeg001 [专家分：14390] 发布于 2010-09-01 09:30:00
我觉得先搞清楚几个东西.
1. 在网上我找 cvf IMSL eig 关键词没有找到这个函数的说明, 不过cvf里面确实可以运行. 楼主你可以把你看到这个函数的相关资料帖出来.
2. 如aliouying所讲, 我查看了手头上v3.0的文档和官网http://www.vni.com上提供的v5.0文档, 没有提到(至少没有直接提到)eig这个函数, 所以对它不了解.
这里帖一页关于IMSL求解本征值本证矢的目录清单
Routines
2.1. Eigenvalue Decomposition
2.1.1 Computes the eigenvalues of a self-adjoint
matrix, A.................................................................LIN_EIG_SELF 432
2.1.2 Computes the eigenvalues of an n n matrix, A.... LIN_EIG_GEN 439
2.1.3 Computes the generalized eigenvalues of an n n
matrix pencil, Av = Bv .........................................LIN_GEIG_GEN 448
2.2. Eigenvalues and (Optionally) Eigenvectors of Ax = x
2.2.1 Real General Problem Ax = x
All eigenvalues.................................................................... EVLRG 455
All eigenvalues and eigenvectors .......................................EVCRG 457
Performance index................................................................EPIRG 460
2.2.2 Complex General Problem Ax = x
All eigenvalues.................................................................... EVLCG 462
All eigenvalues and eigenvectors .......................................EVCCG 464
Performance index................................................................EPICG 467
2.2.3 Real Symmetric Problem Ax = x
All eigenvalues..................................................................... EVLSF 469
All eigenvalues and eigenvectors ........................................EVCSF 471
Extreme eigenvalues ...........................................................EVASF 473
Extreme eigenvalues and their eigenvectors.......................EVESF 475
Eigenvalues in an interval ....................................................EVBSF 478
Eigenvalues in an interval and their eigenvectors ...............EVFSF 480
Performance index................................................................ EPISF 483
2.2.4 Real Band Symmetric Matrices in Band Storage Mode
All eigenvalues.....................................................................EVLSB 485
All eigenvalues and eigenvectors ....................................... EVCSB 487
Extreme eigenvalues ...........................................................EVASB 490
Extreme eigenvalues and their eigenvectors.......................EVESB 492
Eigenvalues in an interval ....................................................EVBSB 495
Eigenvalues in an interval and their eigenvectors ...............EVFSB 498
2.2.5 Complex Hermitian Matrices
All eigenvalues .....................................................................EVLHF 502
All eigenvalues and eigenvectors........................................ EVCHF 505
Extreme eigenvalues........................................................... EVAHF 508
Extreme eigenvalues and their eigenvectors ...................... EVEHF 510
Eigenvalues in an interval ................................................... EVBHF 513
Eigenvalues in an interval and their eigenvectors................EVFHF 515
Performance index ................................................................EPIHF 518
7 楼
tianhy2010 [专家分：60] 发布于 2010-09-01 09:38:00
[quote]我觉得先搞清楚几个东西.
1. 在网上我找 cvf IMSL eig 关键词没有找到这个函数的说明, 不过cvf里面确实可以运行. 楼主你可以把你看到这个函数的相关资料帖出来.
2. 如aliouying所讲, 我查看了手头上v3.0的文档和官网http://www.vni.com上提供的v5.0文档, 没有提到(至少没有直接提到)eig这个函数, 所以对它不了解.
这里帖一页关于IMSL求解本征值本证矢的目录清单
Routines
2.1. Eigenvalue Decomposition
2.1.1 Computes the eigenvalues of a self-adjoint
matrix, A.................................................................LIN_EIG_SELF 432
2.1.2 Computes the eigenvalues of an n n matrix, A.... LIN_EIG_GEN 439
2.1.3 Computes the generalized eigenvalues of an n n
matrix pencil, Av = Bv .........................................LIN_GEIG_GEN 448
2.2. Eigenvalues and (Optionally) Eigenvectors of Ax = x
2.2.1 Real General Problem Ax = x
All eigenvalues.................................................................... EVLRG 455
All eigenvalues and eigenvectors .......................................EVCRG 457
Performance index................................................................EPIRG 460
2.2.2 Complex General Problem Ax = x
All eigenvalues.................................................................... EVLCG 462
All eigenvalues and eigenvectors .......................................EVCCG 464
Performance index................................................................EPICG 467
2.2.3 Real Symmetric Problem Ax = x
All eigenvalues..................................................................... EVLSF 469
All eigenvalues and eigenvectors ........................................EVCSF 471
Extreme eigenvalues ...........................................................EVASF 473
Extreme eigenvalues and their eigenvectors.......................EVESF 475
Eigenvalues in an interval ....................................................EVBSF 478
Eigenvalues in an interval and their eigenvectors ...............EVFSF 480
Performance index................................................................ EPISF 483
2.2.4 Real Band Symmetric Matrices in Band Storage Mode
All eigenvalues.....................................................................EVLSB 485
All eigenvalues and eigenvectors ....................................... EVCSB 487
Extreme eigenvalues ...........................................................EVASB 490
Extreme eigenvalues and their eigenvectors.......................EVESB 492
Eigenvalues in an interval ....................................................EVBSB 495
Eigenvalues in an interval and their eigenvectors ...............EVFSB 498
2.2.5 Complex Hermitian Matrices
All eigenvalues .....................................................................EVLHF 502
All eigenvalues and eigenvectors........................................ EVCHF 505
Extreme eigenvalues........................................................... EVAHF 508
Extreme eigenvalues and their eigenvectors ...................... EVEHF 510
Eigenvalues in an interval ................................................... EVBHF 513
Eigenvalues in an interval and their eigenvectors................EVFHF 515
Performance index ................................................................EPIHF 518[/quote]
很重要，谢谢关系，有没有这本书啊？我看下怎么 函数用
8 楼
yeg001 [专家分：14390] 发布于 2010-09-01 09:43:00
另外,
1.
对于你的函数情况的不清楚, 我对语句
eigenvalue=eig(c,w=eigenvector) 存疑
这里能保证eig()不对c进行修改吗?
2.
还有一个疑问的地方.楼主你对本征值进行了你喜好的排序.
原函数对本征值跟本证矢输出有什么规律? 你重新排序是否也应该对本证矢进行重新排序?
我进行了简单的代码修改
加入
...(楼主变量定义)
[color=FF0000]real(kind=8) :: d(4, 4)[/color]
...(楼主代码)
[color=FF0000]d = c[/color]
eigenvalue=eig(c,w=eigenvector)
[color=FF0000]write(*, *) matmul(d,eigenvector(:,4)) - eigenvalue(4)*eigenvector(:, 4)[/color]
...(楼主代码)
结果:
(-3.330669073875470E-016,0.000000000000000E+000)
(-2.220446049250313E-016,1.387778780781446E-016)
(3.885780586188048E-016,-2.498001805406602E-016)
(3.330669073875470E-016,-3.885780586188048E-016)
Press any key to continue
可以认为得到的结果为零, 也就是第4个本征值跟本证矢满足本证方程.
[color=FF0000]我顺手也把前面的1, 2, 3都测试了, 结果也是零!![/color]
(手头上的本本没有装matlab, 没有弄matlab代码)
9 楼
yeg001 [专家分：14390] 发布于 2010-09-01 09:44:00
回7楼
http://www.vni.com/products/imsl/documentation/#ts
下面就有关于fortran的IMSL Fortran Numerical Library Version 5.0
楼主有eig()函数的说明也请帖出来吧, 是在没有搜索到.
10 楼
tianhy2010 [专家分：60] 发布于 2010-09-01 10:24:00
关于eig的用法还是从书上抄来的，台湾人写的那本书，书上只是给出下面的语句：
对于复型矩阵：eigenvalue=eig(c,w=eigenvector)
对于实型矩阵：eigenvalue=eig(c,v=eigenvector)
其他的我也不清楚了
我来回复
您所在位置：