厄密矩陣對角化中的特徵向量不連續性

Question

我需要對一個2x2厄密矩陣進行對角化，該矩陣取決於連續變化的參數x。爲了對角化我使用EISPACK。當我繪製特徵向量的實部和虛部作爲x的函數時，我注意到它們具有不連續性。特徵值計算是可以的。當我繪製Maxima中的特徵向量時，解決方案看起來是連續的。我需要連續的特徵向量，因爲在下一步我需要計算它們的導數。

在f77代碼下面，我用作測試（在gingw上編譯gfortran）。

program Eigenvalue 

    implicit none 

    integer n, m 
    parameter (n=2) 
    integer ierr, matz, i, j 
    double precision x, dx, xf, amp, xin 
    double precision w(n) 
    double precision Ar(n,n), Ai(n,n) 
    double precision xr(n,n), xi(n,n) 
    double precision fm1(2,n) ! f77 
    double precision fv1(n)  ! f77 
    double precision fv2(n)  ! f77  
    double complex psi1a, psi1b, psi2a, psi2b 

    m = 51    

    xf = 10.d0 
    xin = 0.0d0  

    amp = 2.d0 
    dx = (xf - xin)/(m-1) 

    do i = 1, m  
    x = dx*(i-m) + xf 

    Ar(1,1) = dsin(x)**2 
    Ar(1,2) = amp*dcos(x) 
    Ar(2,1) = amp*dcos(x) 
    Ar(2,2) = dcos(x)**2 

    Ai(1,1) = 0.0d0 
    Ai(1,2) = amp*dsin(x) 
    Ai(2,1) = -amp*dsin(x) 
    Ai(2,2) = 0.0d0 

    matz = 1 

    call ch (n, n, ar, ai, w, matz, xr, xi, fv1, fv2, fm1, ierr) !f77 

    write(20,*) x, w(1), w(2) 

    write(21,*) x, xr(1,1), xi(1,1) 
    write(22,*) x, xr(2,1), xi(2,1) 

    write(23,*) x, xr(1,2), xi(1,2) 
    write(24,*) x, xr(2,2), xi(2,2)  

!  autovetor 1 

    psi1a = cmplx(xr(1,1),xi(1,1)) 
    psi1b = cmplx(xr(1,2),xi(1,2)) 

!  autovetor 2 

    psi2a = cmplx(xr(2,1),xi(2,1)) 
    psi2b = cmplx(xr(2,2),xi(2,2))   

    end do 
    end 

Answer 1

雖然不是一個真正的答案，接下來是我與LAPACK使用的代碼。 我用LAPACK和BLAS的最新版本，具有下列編譯器選項：

gfortran -oG -std = F2008 -Wall -Wextra {} location_of_lapack {/liblapack.a} location_of_blas /blas_LINUX.a main.f90時 - 主

我在Mac OS X上編譯gfortran 6.3.0從自制軟件。

正如伊恩所提到的，像dcos這樣的東西被替換爲cos，並且我使用了KIND =公式來確保相同的精度。

伊恩上面提到的關於任意階段。 回答此問題here;我已將此解決方案轉換爲下面的代碼。 「魔術」發生在致電ZHEEV之後。

有了此修復程序，我看不到任何不連續性。

program Eigenvalue 
    !> This can be used with the f2008 call 
    use, intrinsic :: iso_fortran_env 
    implicit none 

    !> dp contains the kind value for double precision. 
    !> Use below if compiling to f2008 
    integer, parameter :: dp = REAL64 
    !> Use below if compiling with f95 up and comment out iso_fortran_env 
    !>integer, parameter :: dp = SELECTED_REAL_KIND(15, 300) 
    !> Set wp to the desired precision. 
    integer, parameter :: wp = dp 

    integer, parameter :: n = 2 
    integer :: i, j, k, m 
    real(kind=wp) x, dx, xf, amp, xin 
    real(kind=wp), dimension(n) :: w 
    real(kind=wp), dimension(n, n) :: Ar, Ai 
    complex(kind=wp), dimension(n, n) :: A 
    complex(kind=wp), dimension(max(1,2*n-1)) :: WORK 
    integer, parameter :: lwork = max(1,2*n-1) 
    real(kind=wp), dimension(max(1, 3*n-2)) :: RWORK 
    integer :: info 
    complex(kind=wp) :: psi1a, psi1b, psi2a, psi2b 
    real(kind=wp) :: mag 
    m = 51    

    xf = 10.0_wp 
    xin = 0.0_wp 

    amp = 2.0_wp 

    if (m .eq. 1) then 
    dx = 0.0_wp 
    else 
    dx = (xf - xin)/(m-1) 
    end if 

    do i = 1, m 
    x = dx*(i-m) + xf 

    Ar(1,1) = sin(x)**2 
    Ar(1,2) = amp*cos(x) 
    Ar(2,1) = amp*cos(x) 
    Ar(2,2) = cos(x)**2 

    Ai(1,1) = 0.0_wp 
    Ai(1,2) = amp*sin(x) 
    Ai(2,1) = -amp*sin(x) 
    Ai(2,2) = 0.0_wp 

    do j = 1, n 
     do k = 1, n 
      A(j, k) = cmplx(Ar(j, k), Ai(j, k), kind=wp) 
     end do 
    end do 
    call ZHEEV('V', 'U', N, A, N, W, WORK, LWORK, RWORK, INFO) 
    do j = 1, n 
     A(:, j) = A(:, j)/A(1, j) 
     mag = sqrt(real(A(1, j)*conjg(A(1, j)))+ real(A(2, j)*conjg(A(2, j)))) 
     A(:, j) = A(:, j)/mag 
    end do 

    psi1a = A(1, 1) 
    psi1b = A(1, 2) 
    psi2a = A(2, 1) 
    psi2b = A(2, 2) 

    end do 
end program Eigenvalue