"東京1,2月B-T"の14×14ハミルトニアンのエネルギー固有値が縮退するか調べた。
コード
import xlrd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy
from scipy import linalg
import openpyxl
# 変数
tau = 1 * 10 ** (3) # 時間スケール
taup = tau + 1
dt = 1 / tau
time = np.zeros(taup)
h_bar = 1 * 10 ** (-2)
dim = 14 # ハミルトニアンの次元
E = np.zeros((taup, dim), dtype=np.complex)
EI = np.zeros(taup, dtype=np.complex)
# グラフのセッティング
mpl.rcParams['font.family'] = 'Times New Roman' # 使用するフォント名
mpl.rcParams['pdf.fonttype'] = 42 # フォントが埋め込まれるようになる
params = {'backend': 'ps', # バックエンド設定
'axes.labelsize': 24, # 軸ラベルのフォントサイズ
'font.size': 8, # テキストサイズ
'legend.fontsize': 8, # 凡例の文字の大きさ
'xtick.labelsize': 40, # /J, # x軸の数値の文字の大きさ
'ytick.labelsize': 40, # /J, # y軸の数値の文字の大きさ
'lines.linewidth': 3, # グラフデータの線幅
'axes.linewidth': 4, # 軸の線幅
'xtick.color': 'black', # x軸の目盛の色
'ytick.color': 'black', # y軸の目盛の色
'xtick.major.width': 3, # x軸主目盛り線の線幅
'ytick.major.width': 3, # y軸主目盛り線の線幅
'xtick.major.size': 24, # /J, #x軸主目盛り線の長さ
'ytick.major.size': 24, # /J, #y軸主目盛り線の長さ
'xtick.direction': 'in', # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
'ytick.direction': 'in', # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
'axes.facecolor': 'white', # グラフ背景の色
'text.usetex': False, # 使用するフォントをtex用（Type1）に変更
'figure.figsize': [30 / 2.54, 30 / 2.54]
}
mpl.rcParams.update(params)
# ハミルトニアン
H0 = np.zeros( ( dim, dim ), dtype = np.float )
for i in range( 0, dim ):
for j in range(0,dim):
if i == j:
H0[i][j] = 1
else:
H0[i][j] = -1
data = xlrd.open_workbook('/Users/mizuguchitaishi/Documents/B-T/東京1,2月B-T.xlsx')
sheets = data.sheet_by_name('Sheet2')
def get_rowlist_2d(sheets,start_row,end_row,start_col,end_col):
return [sheets.row_values(row,start_col,end_col) for row in range(start_row,end_row+1)]
Htau = -1*np.array(np.reshape(get_rowlist_2d(sheets,30,43,38,52),(14,14))) #気象データによる予測
def H(t):
H = t * Htau + (1 - t) * H0
return H
# 固有値、固有ベクトルを普通に求める
eig_val, eig_vec = scipy.linalg.eigh(H(0))
for i in range(len(eig_vec)):
eig_vec[i] = eig_vec[i] / scipy.linalg.norm(eig_vec[i])
min0 = np.where(eig_val == min(eig_val)) # 固有値が最小となるリストの場所を格納
psi = eig_vec[:, min0] # 最小固有値の固有ベクトルをリストから取り出す
psi = np.reshape(psi, [dim, 1]) # 取り出した固有ベクトルをリストの形から縦ベクトルに直す
# H0の基底状態（量子アニーリングとは関係なく固有値問題を解く）
l, P = scipy.linalg.eigh(H(0))
for i in range(0, dim):
E[0][i] = l[i]
# エネルギー期待値
EI[0] = (np.conj(psi).T @ H(0) @ psi)
# 横軸（時間）
time[0] = 0
# 量子アニーリング
def R(i, psi):
R = 1 / (1j * h_bar) * (H(i) @ psi)
return R
for t in range(1, tau + 1):
# 波動関数を1つ進める(ルンゲクッタ法)
k1 = R((t - 1) * dt, psi)
k2 = R((t - 1) * dt + dt / 2, psi + dt / 2 * k1)
k3 = R((t - 1) * dt + dt / 2, psi + dt / 2 * k2)
k4 = R((t - 1) * dt + dt, psi + dt * k3)
psi = psi + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
#波動関数を1つ進める(expかける)
#psi = scipy.linalg.expm( H( dt * (t-1) ) * dt / (1j *h_bar) ) @ psi
psi = psi/scipy.linalg.norm(psi)
print(psi,"\n")
# 横軸（時間）
time[t] = t * dt
# エネルギー固有値の答え（量子アニーリングとは関係なく固有値問題を解く）
l, P = scipy.linalg.eigh(H(t * dt))
for i in range(0, dim):
E[t][i] = l[i].real
# エネルギー期待値
EI[t] = (np.conj(psi).T @ H(t * dt) @ psi).real
# ターゲットの固有値、固有ベクトルを普通に求める
eig_val, eig_vec = scipy.linalg.eigh(Htau)
for i in range(len(eig_vec)):
eig_vec[i] = eig_vec[i] / scipy.linalg.norm(eig_vec[i])
min_g = np.where(eig_val == min(eig_val)) # 固有値が最小となるリストの場所を格納
psi_tau = eig_vec[:, min_g] # 最小固有値の固有ベクトルをリストから取り出す
psi_tau = np.reshape(psi_tau, [dim, 1]) # 取り出した固有ベクトルをリストの形から縦ベクトルに直す
print(eig_val,"\n")
def v(i):
v = np.reshape(eig_vec[:,i],[dim,1])
return v
for i in range(dim):
print(v(i),"\n")
def E_(i):
E_ = np.conj(v(i)).T @ Htau @ v(i)
return E_
print("E_0",E_(0),"\n","E_1",E_(1))
#print(psi_tau)
# 描画
# ハミルトニアンの固有値と量子アニの期待値の時間発展グラフ
for i in range(0,2):
plt.plot(time, E[:, i], ".", markersize=12)
plt.plot(time, EI, ".", markersize=12)
plt.tick_params(labelsize=36) # 目盛りの文字サイズ
plt.tick_params(width=3, length=20) # 目盛りのサイズ（幅、長さ）
plt.show()