回 帖 发 新 帖 刷新版面

主题:[讨论]关于矩阵特征向量(fortran与matlab对照)

最近一直很疑惑,为什么对于复数矩阵,得到的特征向量会有差距呢?尽管从数学上讲,对于一个矩阵,特征向量有无穷多个,但是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个回复)

31 楼

[quote]
[em35][em35][em35][em35][em35][em35][em35][em35][em46][em46][em46][em46][em46][em46][em52][em52][em52][em52]
今年第一次编程[/quote]

所以需要花大力气修改你的程序,不要到头来对着一堆错误的计算结果,得出一通错误的结论![em54][em55][em56][em58][em75][em71][em52][em22][em21][em23][em24][em19][em18][em20][em17][em7][em15][em14][em13][em11]

32 楼

[quote][quote]
[em35][em35][em35][em35][em35][em35][em35][em35][em46][em46][em46][em46][em46][em46][em52][em52][em52][em52]
今年第一次编程[/quote]

所以需要花大力气修改你的程序,不要到头来对着一堆错误的计算结果,得出一通错误的结论![em54][em55][em56][em58][em75][em71][em52][em22][em21][em23][em24][em19][em18][em20][em17][em7][em15][em14][em13][em11][/quote]
静听教诲!!我发现程序中有个地方不完善。由于矩阵h有很多特征向量,我选择能满足条件的特征向量,要求就是特征向量与它的转置的乘积为1。你看我这样修改好不好:

eps=1.0e-8
do ii=1,4
     if(dot_product(eigenvector(:,j)*transpose(eigenvector(:,j)))<=1.0-eps)  then
       kk=1
    else
      kk=0
       if(kk==1) then
!求矩阵特征值和特征向量
 eigenvalue=eig(c,w=eigenvector) 
 !特征值特征向量同时排序
  do i=1,4
   do j=i+1,4
    if(eigenvalue(i)>eigenvalue(j)) then
      temp=eigenvalue(j)
      eigenvalue(j)=eigenvalue(i)
      eigenvalue(i)=temp

      eigenvectortmp(:,j)=eigenvector(:,j)
      eigenvector(:,j)=eigenvector(:,i)
      eigenvector(:,i)=eigenvectortmp(:,j)
    endif
   end if
  end if
   end do
  end do
end do




!将特征值和特征向量排序,便于计算
先找满足条件的特征向量,找到之后再对 特征值和特征向量排序。

我来回复

您尚未登录,请登录后再回复。点此登录或注册