2018/12/23

"東京1,2月B-T"の14×14ハミルトニアンのエネルギー固有値が縮退するか調べた。









コード

```python
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()

```