3

我在 fortran 中编写了以下程序，该程序使用名为ZGEEV的 lapack 子例程。这个想法是看矩阵的特征值如何矩阵对角化 随着 k 从实数变为复数而变化。分析地，答案应该是 2 和 0，无论 k 是否复杂。但是我得到了一个显示很多变化的图。
特别是对于真正的 k，情节看起来像这样 -
绘制特征值 vs k 这是我写的代码 - 

           program main
             implicit none
       !**********************************************                 
             complex(8) :: k,mat(2,2)
             complex(8) :: eigenvals(2)
             real(8), parameter    :: kmax = 2.d0 
             real(8), parameter    :: dk = 1.d-1 
             real(8)               :: kr,ki
       !**********************************************
             kr=-kmax
             do while (kr.le.kmax)
                ki= -1.d-3
                do while (ki.le.1.d-3)
                   k=cmplx(kr,ki)
                   call init_mat(k,mat)
                   call diagonalize(mat,eigenvals)
                   print*, real(k), real(eigenvals(2)),aimag(eigenvals(2))
                   ki=ki+1.d-4
                end do
                kr=kr+dk
             end do
           end program main





           subroutine init_mat(k,mat)
             implicit none
             complex(8),intent(in) :: k
             complex(8),intent(out):: mat(2,2)
             complex(8),parameter  :: di=(0.d0,1.d0)
             complex(8),parameter  :: d1=(1.d0,0.d0)
        !**********************************************
             mat(1,1) = d1 
             mat(1,2) = exp(di*k)
             mat(2,1) = exp(di*k) 
             mat(2,2) = d1 
             return
           end subroutine init_mat

           subroutine diagonalize(mat,eigenvals)
             implicit none
             complex(8),intent(in) :: mat(2,2)
             complex(8),intent(out):: eigenvals(2)
             complex(8)            :: vl(2,2),vr(2,2)
             complex(8),allocatable:: work(:)
             integer(4)            :: lwork
             complex(8)            :: rwork(4)
             complex(8)            :: mat2(2,2)
             integer(4)            :: info
        !**********************************************
             mat2(:,:) = mat(:,:)
             allocate(work(6))
             call zgeev('N', 'N', 2, mat2, 2, eigenvals, vl, 2, vr, 2, work, -1, rwork, info)
             lwork = work(1)
             deallocate(work)
             allocate(work(lwork))
             call zgeev('V', 'V', 2, mat2, 2, eigenvals, vl, 2, vr, 2, work, lwork, rwork, info)
             if (info.ne.0) print*, info
             stop 'diagonalize failed'
           end subroutine diagonalize

欢迎在评论中提出任何关于这种异常原因的懒惰理论！

PS：我在python中写了一个类似的代码，特征值是y = 2和y = 0的两条恒定线。

4

1 回答 1

3

在子程序 init_mat(k,mat)

mat(1,2) = exp(di*k) 和 mat(2,1) = exp(di*k)

但是其中之一，例如 mat(2,1) 应该 = exp(-di*k)

尽管您的数学项目需要一个在非对角线上带有 e^ik 和 e^-ik 的矩阵，但显示的代码是在两个非对角线上创建一个带有 e^ik 的矩阵。实际编码的矩阵具有复杂的特征值，因此用于查找特征值的子程序可能工作正常，并且所示的输入有错误的规范。

那么 [[1, e^ik], [e^ik, 1]] 的特征值是什么？

好吧，迹仍然是 2，所以特征值总和为 2。

并且行列式是 1-e^(2ik)，所以乘积是复数。

这表明实际输入的矩阵的特征值是复共轭，总和为 2。通过检查，特征值似乎是 1 +/- e^ik

于 2013-06-20T08:36:00.047 回答