匯入標頭檔案
import numpy as np
import numpy.matlib
from scipy.linalg import block_diag
from numpy import linalg as la
import matplotlib.pyplot as plt
from math import sqrt,pi
# 生成塊三對角矩陣函式
def tridiag(c, u, d, n):
# c, u, d are center, upper and lower blocks, repeat n times
cc = block_diag(*([c]*n))
shift = c.shape[1]
uu = block_diag(*([u]*n))
uu = np.hstack((np.zeros((uu.shape[0], shift)), uu[:,:-shift]))
dd = block_diag(*([d]*n))
dd = np.hstack((dd[:,shift:],np.zeros((uu.shape[0], shift))))
return cc+uu+dd
def hamiltonian_h0_su4(k,n,t=-1):
"""t 是最近鄰hopping係數
k 是bloch波向量
n 是寬度
"""a = np.matrix([[sqrt(3),0,1,np.exp(-1j*k)],[0, sqrt(3),1,1],[1,1,sqrt(3),0],[np.exp(1j*k),1,0,sqrt(3)]])
b = np.matrix([[0,0,0,0],[0,0,0,0],[1,0,0,0],[0,-1,0,0]])
h = tridiag(a,b,b.h,n)
h[0,0] = 1000
return h
nk = 256
n = 256
k = np.linspace(0,2*pi,nk)
band = np.zeros((4*n, nk))
ak = np.zeros((4*n,nk),dtype="complex")
for i in range(nk):
hk0 = hamiltonian_h0_su4(k[i],n)
e, a = la.eigh(hk0)
band[:,i] = e
ak[:,i] = a[:,n-1]
for i in range(4*n):
plt.plot(k,band[i,:])
plt.ylim((-1,1))
ak = ak[:16,:]def hamiltonian_heff(nk,iq,ak,u,n):
h = np.zeros((nk,nk),dtype='complex')
for j in range(nk):
j_q = (j - iq + nk)%nk
for i in range(nk):
i_q = (i-iq + nk)%nk
h[j,i] -= np.sum(ak[:,i_q].conj()*ak[:,i]*ak[:,j_q]*ak[:,j].conj())
m = h.copy()
for i in range(nk):
i_q = (i - iq + nk)%nk
h[i,i] += np.sum(np.power(np.absolute(ak),2).sum(axis=1)*np.power(np.absolute(ak[:,i]),2))
return h/n, m/n
u = 1.0
band_mag = np.zeros((nk, nk))
band_m = np.zeros((nk,nk))
psi = np.zeros((nk,nk,nk),dtype='complex')
psi_m = np.zeros((nk,nk,nk),dtype='complex')
for i in range(nk):
heff,m = hamiltonian_heff(nk,i,ak,u,n)
e, a = la.eigh(heff)
em, am = la.eigh(m)
psi[:,:,i] = a
band_mag[:,i] = e
band_m[:,i] = em
psi_m[:,:,i] = am
for i in range(nk):
plt.plot(k,band_mag[i,:],color='gray')
#plt.ylim((0,0.1))
for i in range(2):
plt.plot(k,band_m[i,:],color='gray')
def circshift(matrix,shiftnum1,shiftnum2):
""""""
h,w=matrix.shape
matrix=np.vstack((matrix[(h-shiftnum1):,:],matrix[:(h-shiftnum1),:]))
matrix=np.hstack((matrix[:,(w-shiftnum2):],matrix[:,:(w-shiftnum2)]))
return matrix
delta = 0.002
nw = 128
w = np.linspace(0,0.05,nw)
spectral_a = np.zeros((nw,nk),dtype='complex')
for i in range(nk):
a_kq = circshift(ak,0,i)
s_plus = np.sum(ak.conj()*a_kq,axis=0).reshape(nk,1)
for j in range(nw):
for p in range(nk):
spectral_a[j,i] += np.absolute(np.dot(psi_m[:,0,i].t,psi[:,p,i]))**2/(w[j]-band_mag[p,i]+1j*delta)
x,y = np.meshgrid(k,w)
plt.figure(figsize=(5,5))
plt.pcolormesh(x, y, -spectral_a.imag/pi)
plt.colorbar()
#plt.yscale("log")
#plt.ylim(0.02, 0.05)
python矩陣對角化 大矩陣對角化python
我使用scipy中的linalg來得到155x156矩陣的e值和特徵向量。然而，與矩陣相比，特徵值的階數似乎是隨機的。我希望第乙個特徵值對應於矩陣中的第乙個數。請看下面我的例行程式。我首先讀取的是乙個包含所有浮點數的檔案 1u1o.dat 2533297.650278 2373859.531153 ...
對角化求可逆矩陣 矩陣對角化方法
矩陣對角化方法 摘要 本文給出了一種不同於傳統方法的矩陣對角化方法，利用矩陣的初等變換，先求出矩陣的特徵根與特徵向 量，接著再判斷矩陣是否可對角化。矩陣特徵根 特徵向量 對角化the methods of the diagonalization of the matrix gabstract in ...
相合以及對角化
介紹相合的定義以及其引出的標準型.定義 1 設給定 a,b in m n 如果存在乙個非奇異的矩陣 s 使得 a b sas 那麼就說 b 與 a 是 相合的 也稱為星相合的或者共軛相合的 b b sas t 那麼就說 b 與 a 是相合的，或者 t 相合的 也稱為 t 相合的 乙個矩陣乘以非奇異矩...