投稿日:

【第４話】強制振り子の計算アルゴリズム【ルンゲクッタで行こう！④】

############################################################################
#　強制振り子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0

#支点の初期位置
box_r0 = np.array([0.0, 0.0, 0.0])

#初期位置と速度
r0 = np.array([0.0, 0.0, -1.0])
v0 = np.array([0.0, 0.0, 0.0])

#位置ベクトル
L0 = r0 - box_r0
#ひもの長さ
L0_abs = np.sqrt( np.sum( L0**2 ) )

#強制振動の大きさと角振動数
box_l = 0.2
box_omega = 0.5 * (2.0 * np.pi) 

#支点の運動関数
def boxFunction( t ):
	box_rx = box_r0[0] + box_l * np.sin( box_omega * t)
	box_vx = box_l * box_omega * np.cos( box_omega * t)
	box_ax = - box_l * box_omega**2 * np.sin( box_omega * t)
	box_r = np.array([box_rx, 0.0, 0.0])
	box_v = np.array([box_vx, 0.0, 0.0])
	box_a = np.array([box_ax, 0.0, 0.0])
	return box_r, box_v, box_a


#重力加速度
g = np.array([0.0, 0.0, -10.0])
#補正パラメータ
compensationK = 0.0
compensationGamma = 0.0

#張力ベクトル
def S ( t, r_t, v_t):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)	
	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum(L_t**2) )

	return -m / L_abs ** 2 * ( np.sum(v_bar**2) - np.dot(box_a, L_t) + np.dot(g, L_t) ) * L_t

def F ( t, r_t, v_t ):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)

	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum( L_t**2 ) )
	#おもり速度の大きさ
	v_bar_abs = np.sqrt( np.sum( v_bar**2 ) )
	#単位方向ベクトル
	L_hat = L_t.copy() / L_abs
	#張力ベクトル
	S_ = S( t, r_t, v_t ) 
	#張力の大きさ
	S_coefficient = np.sqrt( np.sum( S_**2 ) )

	#補正力ベクトル
	compensationK = S_coefficient 
	compensationGamma = v_bar_abs * 10
	fc1 = - compensationK * ( L_abs - L0_abs ) * L_hat
	fc2 = - compensationGamma * np.dot(v_bar, L_hat) * L_hat

	#合力ベクトル
	return  S_ + m * g + fc1 + fc2


#時間区間
t_min = 0
t_max = 40.0
#時間区間数
NT = 4000
skip =  2
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)

		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	box_r, box_v, box_a = boxFunction(t)
	print( "{:.2f}".format(t),  np.sqrt((r[0]-box_r[0])**2 + (r[2]-box_r[2])**2) - 1  )

	#各コマを描画
	img  = plt.plot([box_r[0]], [box_r[2]], colors[0], marker = 's', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([box_r[0], r[0]], [box_r[2], r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv

plt.title( u"強制単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()
投稿日:

【第26話】【結果は渋いです】２次元電子波束を斜めに薄膜へ（一応トンネル効果）

#################################################################
## 薄膜へ照射したガウス波束（トンネル効果）
#################################################################
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mat
import matplotlib.animation as animation
from multiprocessing import Pool

flag_calculate = False
flag_animation = True

#図全体
aspect = 0.9
fig = plt.figure(figsize=(6, 6 * aspect))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#並列数
ParallelNumber = 6
#フォルダ指定
dir_path = "output30/"
#フォルダ生成
os.makedirs(dir_path, exist_ok = True)
#################################################################
## 物理定数
#################################################################
#プランク定数
h = 6.62607015 * 1.0E-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.1093837015 * 1.0E-31
#電子ボルト
eV = 1.602176634 * 1.0E-19
#虚数単位
I = 0.0 + 1.0j

######################################
#　物理系の設定
######################################
#波束の中心エネルギー
E0 = 10.0 * eV
#入射角度
theta0 = 30.0 / 180 * np.pi
#重ね合わせの数
NK = 100
#ガウス分布の幅
sigma = 0.1 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
kx0 = k0 * np.cos(theta0)
ky0 = k0 * np.sin(theta0)

#空間刻み間隔
nm = 1E-9
#空間刻み数
NX = 150
NY = 150
#壁の高さ
V1 = 0.0 * eV
V2 = 10.5 * eV
V3 = 0.0 * eV
#壁の幅
d = 1.0 * nm

#計算時間の幅
ts = -200
te = 200
#時間間隔
dt = 1.0 * 10**-16

#空間幅
x_min = -30.0 * nm
x_max =  30.0 * nm
y_min = x_min
y_max = x_max

######################################
#　反射係数と透過係数を計算する関数定義
######################################
def ReflectionCoefficient( MM ):
	return - MM[1][0] / MM[1][1] 
def TransmissionCoefficient(MM):
	return (MM[0][0] * MM[1][1] - MM[0][1] * MM[1][0] ) / MM[1][1]
def calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d):
	# 転送行列の定義
	M21 = np.zeros((2,2), dtype=np.complex)
	T2  = np.zeros((2,2), dtype=np.complex)
	M32 = np.zeros((2,2), dtype=np.complex)

	# 境界
	M21[0][0] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
	M21[0][1] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
	M21[1][0] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
	M21[1][1] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
	M32[0][0] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
	M32[0][1] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
	M32[1][0] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
	M32[1][1] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
	T2[0][0] = np.exp( I * k2 * cos2 * d )
	T2[0][1] = 0
	T2[1][0] = 0
	T2[1][1] = np.exp( - I * k2 * cos2 * d )
	#転送行列の計算
	M31 = M32 @ T2 @ M21
	#反射係数と透過係数の計算
	rc = ReflectionCoefficient(M31)
	tc = TransmissionCoefficient(M31)
	return rc, tc, M21

#透過係数と反射係数と転送行列の配列
rcs = np.zeros((NK, NK))
tcs = np.zeros((NK, NK))
rcs = np.array(rcs, dtype=complex)
tcs = np.array(tcs, dtype=complex)
M21s = [ 0 ] * NK
for n1 in range(NK):
	M21s[ n1 ] = [ np.zeros((2,2), dtype=np.complex) ] * NK
for kxn in range(NK):
	for kyn in range(NK):
		kx = kx0 + dk * (kxn - NK/2)
		ky = ky0 + dk * (kyn - NK/2)
		k = np.sqrt( kx**2 + ky**2 )
		k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
		k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
		k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 ))
		sin = ky / k
		cos = kx / k
		sin1 = sin * k / k1
		sin2 = sin * k / k2
		sin3 = sin * k / k3
		cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
		cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))
		cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 ))
		rcs[kxn][kyn], tcs[kxn][kyn], M21s[kxn][kyn] = calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d)

######################################
# 波動関数の計算
######################################
#ガウス波束の初期位置
x0 = 0 #-20.0 * nm * np.cos(theta0)
y0 = 0 #-20.0 * nm * np.sin(theta0)

#ガウス波束の値を計算する関数
def Psi( x, y, t ):
	#波動関数値の初期化
	psi = 0 + 0j
	#各波数ごとの寄与を足し合わせる
	for kxn in range(NK):
		for kyn in range(NK):
			#各波数を取得
			kx = kx0 + dk * (kxn - NK/2)
			ky = ky0 + dk * (kyn - NK/2)
			k = np.sqrt( kx**2 + ky**2 )
			k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
			k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
			k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 ))
			sin = ky / k
			cos = kx / k
			sin1 = sin * k / k1
			sin2 = sin * k / k2
			sin3 = sin * k / k3
			cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
			cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))
			cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 ))
			k1x = k1 * cos1
			k1y = k1 * sin1
			k2x = k2 * cos2
			k2y = k2 * sin2
			k3x = k3 * cos3
			k3y = k3 * sin3

			#波数から各振動数を取得
			omega = hbar / (2.0 * me) * ( kx**2  + ky**2 )

			W = np.exp( - I * k1x * x0 - I * k1y * y0) * np.exp( -( (kx - kx0) / (2.0 * sigma) )**2  -( (ky - ky0) / (2.0 * sigma) )**2 )
			if( x < 0 ):
				psi += np.exp( I * ( k1x * x + k1y * y ) - I * omega * t ) * W
				psi += rcs[kxn][kyn] * np.exp( I * ( - k1x * x + k1y * y ) - I * omega * t ) * W
			elif( x < d ):
				A2_p = M21s[kxn][kyn][0][0] + M21s[kxn][kyn][0][1] * rcs[kxn][kyn]
				A2_m = M21s[kxn][kyn][1][0] + M21s[kxn][kyn][1][1] * rcs[kxn][kyn]
				psi += A2_p * np.exp( I * (   k2x * x + k2y * y ) - I * omega * t ) * W
				psi += A2_m * np.exp( I * ( - k2x * x + k2y * y ) - I * omega * t ) * W
			elif( d <= x ):
				psi += tcs[kxn][kyn] * np.exp( I * ( k3x * (x - d) + k3y * y ) - I * omega * t ) * W

	return psi

#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0, 0 ))

######################################
# 空間分布の計算
######################################
#実時間
t = 0
#x軸・y軸のメッシュ生成
xs = np.linspace(x_min, x_max, NX)
ys = np.linspace(y_min, y_max, NY)
X, Y = np.meshgrid(xs, ys)

#並列用ラッパー関数
def wrap_timeEvolution( data ):
	return timeEvolution( *data ) 

#時間発展の計算
def timeEvolution( filepath , tn ):
	#実時間
	t = ( ts + tn ) * dt
	Z = np.sin( np.complex128(X/nm) )
	for ix in range(NX):
		for iy in range(NY):
			x = xs[ix]
			y = ys[iy]
			Z[iy][ix] = Psi( x, y, t ) / psi_abs

	#ファイルへ保存
	np.savetxt( filepath, abs(Z))

#アニメーションフレーム作成用
def update( tn ):
	#グラフ消去
	fig.clear()
	#両軸位置の設定
	axes = fig.add_axes(    [0.07, 0.1/aspect, 0.78, 0.78/aspect])
	#カラーバー位置の設定
	bar_axes = fig.add_axes([0.87, 0.1/aspect, 0.05, 0.78/aspect])
	#データ読み込み	
	filepath = dir_path + str(tn + ts) + ".txt"
	Z = np.loadtxt( filepath )
	#透過領域の振幅を定数倍
	for ny in range(len(X)):
		for nx in range(len(X[ny])):
			if(X[ny][nx] > 1.0 * nm ):
				Z[ny][nx] = Z[ny][nx] * 1000

	Z = Z[:-1, :-1]

	#２次元カラーグラフ描画
	img = axes.pcolormesh( X/nm, Y/nm, Z, vmin=0.0, vmax=1.0, cmap=cmap)
	#カラーバーの追加
	fig.colorbar(img, cax=bar_axes)



####### メイン関数として実行 ########
if __name__ == "__main__":

	#計算開始
	if(flag_calculate):
		list = []
		for tn in range( te - ts + 1 ):
			filepath = dir_path + str(tn + ts) + ".txt"
			list.append( (filepath, tn) )

		with Pool(processes=ParallelNumber) as p:
			p.map(wrap_timeEvolution, list)

	#アニメーション作成
	if( flag_animation ):
		#カラーマップの配色配列
		color_list = []
		color_list.append( [ 0,   "black" ] )
		color_list.append( [ 1,   "white" ] )

		#カラーマップ用オブジェクトの生成
		cmap = mat.colors.LinearSegmentedColormap.from_list('cmap', color_list)

		#アニメーションの設定
		ani = animation.FuncAnimation( fig, update, range(te - ts + 1), interval=10)
		#アニメーションの保存
		#ani.save("output.html", writer=animation.HTMLWriter())
		#ani.save('anim.gif', writer='pillow')
		ani.save("output.mp4", writer="ffmpeg", dpi=300)

		#グラフの表示
		#plt.show()
投稿日:

【第３話】拘束条件の改良！【ルンゲクッタで行こう！③】

20210.08.27 アルゴリズムを改良しました。破綻しません。

############################################################################
#　単振り子運動シミュレーション（補正力あり）
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#ひもの長さ
r_abs0 = 1.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])
#補正パラメータ
compensationK = 100.0
compensationGamma = 0.0

#張力ベクトル
def S ( t, r_t, v_t ):
	r_abs = np.sqrt( np.sum(r_t**2) )
	return -m / r_abs ** 2 * ( np.sum(v_t**2) + np.dot(g, r_t) ) * r_t

def F ( t, r_t, v_t ):
	r = np.sqrt( np.sum( r_t**2 ) )
	v = np.sqrt( np.sum( v_t**2 ) )
	r_hat = r_t.copy() / r
	v_hat = v_t.copy() / v
	fc  = - compensationK * ( r - r_abs0 ) * r_hat
	fc += - compensationGamma * np.dot(v, r_hat) * r_hat
	return S( t, r_t, v_t ) + m * g + fc


#初期位置と速度
r0 = np.array([0.0 , 0.0, -r_abs0])
v0 = np.array([6.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT = 1000
skip =  1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)	
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	print( "{:.2f}".format(t),  np.sqrt(r[0]**2 + r[2]**2) - 1  )
	#各コマを描画
	img  = plt.plot([0], [0], colors[0], marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([0, r[0]], [0, r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv



plt.title( u"単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
#ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
#plt.show()
投稿日:

【第２５話】量子ビットの作り方３：静電場の与えたときの電子状態【量子ドットコンピュータの原理③】

############################################################################
#　無限に深い量子井戸中の電子に静電場を加えたときの固有状態（ステップ３）
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
import numpy as np
import numpy.linalg as LA

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19
#電気素量
e = 1.60217733 * 10**-19
#複素数
I = 0.0 + 1.0j

######################################
#　物理系の設定
######################################
#量子井戸の幅
L = 1 * 10**-9
#計算区間
x_min = -L / 2.0
x_max = L / 2.0
#状態数
n_max = 10
#行列の要素数
DIM = n_max + 1
#空間分割数
NX = 500
#空間刻み間隔
nm = 1.0 * 10**-9
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#電場の強さ
Ex = 1.0 * 10**10
#電場の強さ
Ex_max = 1.0 * 10**10
#電場強度の分割数
NEx = 10

#静電場無しの固有関数
def verphi(n, x):
	kn = np.pi * (n + 1) / L
	return np.sqrt(2.0 / L) * np.sin(kn * (x + L / 2.0))

#ポテンシャル項
def V(x, Ex):
	return(e * Ex * x)
	
#静電場無しの固有エネルギー
def Energy(n):
	kn = np.pi * (n + 1) / L
	return hbar * hbar * kn * kn / (2.0 * me)

#被積分関数
def integral_matrixElement(x, n1, n2, Ex):
	return verphi(n1 ,x) * V(x, Ex) * verphi(n2, x) / eV


######################################
#　展開係数の計算
######################################
#固有値・固有ベクトルの初期化
eigenvalues = [0] * (NEx + 1)
vectors = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	eigenvalues[nEx] = []
	vectors[nEx] = []




#存在確率分布グラフ描画用の配列初期化
xs = []
phi = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	phi[nEx] = [0] * 2


#中心の電場依存性グラフ描画用の配列初期化
averageX = [0] * 2
for n in range(len(averageX)):
	averageX[n] = [0] * (NEx + 1)


######################################
#　静電場を加えた固有関数
######################################
def verphi_Ex(nEx, n, x):
	phi = 0 + 0j
	for m in range(n_max+1):
		phi += vectors[nEx][n][m] * verphi(m, x)	
	return phi


######################################
#　#静電場強度ごとの固有関数の計算
######################################
for nEx in range(NEx + 1):
	print("電場強度：" + str( nEx * 100 / NEx ) + "%")
	#静電場の強度を設定
	Ex = Ex_max / NEx * nEx
	#エルミート行列（リスト）
	matrix=[]

	###行列要素の計算
	for n1 in range(n_max + 1):
		col=[]
		for n2 in range(n_max + 1):
			#ガウス・ルジャンドル積分
			result = integrate.quad(
				integral_matrixElement, #被積分関数
				x_min, x_max,		        #積分区間の下端と上端
				args=(n1, n2, Ex)			  #被積分関数へ渡す引数
			)
			real = result[0]
			imag = 0j
			#無静電場のエネルギー固有値（対角成分）
			En = Energy(n1)/eV if (n1 == n2) else 0
			#行の要素を追加
			col.append( En + real )
		#行を追加
		matrix.append( col )

	#リスト → 行列
	matrix = np.array( matrix )

	###固有値と固有ベクトルの計算
	result = LA.eig( matrix )
	eig = result[0] #固有値
	vec = result[1] #固有ベクトル

	#小さい順に並べるためのインデックス（配列）
	index = np.argsort( eig )
	#固有値を小さい順に並び替え
	eigenvalues[nEx] = eig[ index ]

	#転置行列
	vec = vec.T
	#固有ベクトルの並び替え
	vectors[nEx] = vec[ index ]

	### 検算：MA-EA=0 ?
	sum = 0
	for i in range(DIM):
		v = matrix @ vectors[nEx][i] - eigenvalues[nEx][i] * vectors[nEx][i]
		for j in range(DIM):
			sum += abs(v[j])**2
	print("|MA-EA| =" + str(sum))


#アニメーション作成用
ims = []
for nEx in range(NEx + 1):
	###固有関数の空間分布
	phi[nEx][0] = abs(verphi_Ex(nEx, 0, x))**2 / (1.0 * 10**9)
	phi[nEx][1] = abs(verphi_Ex(nEx, 1, x))**2 / (1.0 * 10**9)
	#各コマを描画
	img  = plt.plot(x/nm, phi[nEx][0], colors[0], linestyle='solid', linewidth = 5 )
	img += plt.plot(x/nm, phi[nEx][1], colors[1], linestyle='solid', linewidth = 5 )
	time = plt.text(0.4, 4.1, "電場強度：" + str( nEx * 100 / NEx ) + "%" , ha='left', va='center', fontsize=12, fontname="Yu Gothic")
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )


#########################################################################
# グラフの描画（固有関数）
#########################################################################
plt.title( u"無限に深い量子井戸に静電場を加えたときの固有状態", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ |\varphi^{(m)}(x)|^2$", fontsize=30)
#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.12, top = 0.95)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-0.6, 0.6])
plt.ylim([0, 4.0])

#井戸の概形
plt.vlines([-0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0],  -0.5, 0.5, "black", linestyles='solid', linewidth = 10 )


#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()
投稿日:

【第２話】単振り子の作り方【ルンゲクッタで行こう！②】

20210.08.27 アルゴリズムを改良しました。破綻しません。

############################################################################
#　単振り子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#ひもの長さ
r = 1.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#張力ベクトル
def S ( t, r_t, v_t ):
	vv = np.sum(v_t**2)
	bb = np.dot(g, r_t)
	r_abs = np.sqrt( np.sum(r_t**2) )
	S_abs =  -m * (vv + bb) / (r_abs ** 2)
	return S_abs * r_t 

def F ( t, r_t, v_t ):
	return S( t, r_t, v_t ) + m * g

#初期位置と速度
r0 = np.array([0.0 , 0.0, -r])
v0 = np.array([6.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT = 1000
skip = 1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	print( "{:.2f}".format(t),  np.sqrt(r[0]**2 + r[2]**2) - 1  )
	#各コマを描画
	img  = plt.plot([0], [0], colors[0], marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([0, r[0]], [0, r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv


plt.title( u"単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
#ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
#plt.show()
投稿日:

【第24話】静電場の与え方【量子ドットコンピュータの原理②】

############################################################################
#　無限に深い量子井戸中の電子に静電場を加えたときの固有状態（ステップ１）
############################################################################
import math
import cmath
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate

#図全体
fig1 = plt.figure(figsize=(10, 6))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 12 #フォントサイズ

#複素数
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * math.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19
#電気素量
e = 1.60217733 * 10**-19

######################################
#　物理系の設定
######################################
#量子井戸の幅
L = 1.0 * 10**-9
#計算区間
x_min = -L / 2.0
x_max = L / 2.0
#状態数
n_max = 10
#電場の強さ
Ex = 1.0*10**10

#固有関数
def verphi(n, x):
    kn = math.pi * (n + 1) / L
    return math.sqrt(2.0 / L) * math.sin(kn * (x + L / 2.0))

#ポテンシャル項
def V(x, Ex):
    return(e * Ex * x)

#被積分関数
def integral_matrixElement(x, n1, n2, Ex):
    return verphi(n1 ,x) * V(x, Ex) * verphi(n2, x) / eV

###行列要素の計算
for n1 in range(n_max + 1):
    for n2 in range(n_max + 1):
        #ガウス・ルジャンドル積分
        result = integrate.quad(
            integral_matrixElement, #被積分関数
            x_min, x_max,           #積分区間の下端と上端
            args=(n1,n2, Ex)            #被積分関数へ渡す引数
        )
        real = result[0]
        imag = 0
        #ターミナルへ出力
        print( "(" + str(n1) + ", " + str(n2) + ")  " + str(real))
投稿日:

【第23話】正規直交系の有難みって？【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 無限に深い量子井戸の固有関数の直交性の確認
#################################################################
import numpy as np
import scipy.integrate as integrate

#################################################################
## 物理定数
#################################################################
#プランク定数
h = 6.62607015 * 1.0E-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.1093837015 * 1.0E-31
#電子ボルト
eV = 1.602176634 * 1.0E-19

#################################################################
## 物理系の設定
#################################################################
#量子井戸の幅（m）
L = 1.0 * 10**-9
#計算区間
x_min = -L / 2.0
x_max = L / 2.0

#固有関数
def verphi(n, x):
    kn = np.pi * (n + 1) / L
    return np.sqrt(2.0 / L) * np.sin(kn * (x + L / 2.0))

#被積分関数
def integral_orthonomality(x, n, m):
    return verphi(m, x) * verphi(n, x)

#状態数
NS = 3
for n in range( NS + 1 ):
    for m in range( NS + 1 ):
        #ガウス・ルジャンドル積分
        result = integrate.quad(
            integral_orthonomality, #被積分関数
            x_min, x_max,           #積分区間の下端と上端
            args=( n, m )           #被積分関数へ渡す引数
        )
        #ターミナルへ出力
        print( "(" + str(n) + ", " + str(m) + ")  " + str( result[0]) )
投稿日:

【第１話】ルンゲクッタで行こう！①【Pythonコピペで古典力学完全攻略マニュアル】

############################################################################
#　単振動運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt

#図全体
fig = plt.figure(figsize=(12, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理定数
#################################################################
#質量
m = 1.0
#ばね定数
k = 1.0

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
    #コンストラクタ
    def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
        self.r = r0.copy()
        self.v = v0.copy()
        self.dr = np.array([0.0, 0.0, 0.0])
        self.dv = np.array([0.0, 0.0, 0.0])
        self.dt = dt
        self.__v1 = np.array([0.0, 0.0, 0.0])
        self.__v2 = np.array([0.0, 0.0, 0.0])
        self.__v3 = np.array([0.0, 0.0, 0.0])
        self.__v4 = np.array([0.0, 0.0, 0.0])
        self.__a1 = np.array([0.0, 0.0, 0.0])
        self.__a2 = np.array([0.0, 0.0, 0.0])
        self.__a3 = np.array([0.0, 0.0, 0.0])
        self.__a4 = np.array([0.0, 0.0, 0.0])

    #速度を与えるメソッド
    def V(self, t, r, v):
        return v.copy()
    ##########################################################
    #加速度を与えるメソッド
    def A(self, t, r, v):
        return -r.copy() * k / m
    ##########################################################

    #時間発展を計算するメソッド
    def timeEvolution(self, t):
        #１段目
        self.__v1 = self.V(
            t, 
            self.r, 
            self.v
        )
        self.__a1 = self.A(
            t, 
            self.r, 
            self.v
        );
        #２段目
        self.__v2 = self.V(
            t + self.dt / 2.0, 
            self.r + self.__v1 * self.dt / 2.0, 
            self.v + self.__a1 * self.dt / 2.0
        )
        self.__a2 = self.A(
            t + self.dt / 2.0, 
            self.r + self.__v1 * self.dt / 2.0, 
            self.v + self.__a1 * self.dt / 2.0            
        )
        #３段目
        self.__v3 = self.V(
            t + self.dt / 2.0, 
            self.r + self.__v2 * self.dt / 2.0, 
            self.v + self.__a2 * self.dt / 2.0                    
        )
        self.__a3 = self.A(
            t + self.dt / 2.0, 
            self.r + self.__v2 * self.dt / 2.0, 
            self.v + self.__a2 * self.dt / 2.0            
        )
        #４段目
        self.__v4 = self.V(
            t + self.dt, 
            self.r + self.__v3 * self.dt, 
            self.v + self.__a3 * self.dt                    
        );
        self.__a4 = self.A(
            t + self.dt, 
            self.r + self.__v3 * self.dt, 
            self.v + self.__a3 * self.dt            
        );        
        #差分の計算
        self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
        self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 物理系の設定
#################################################################

#初期位置と速度
r0 = np.array([0.0, 0.0, 10.0])
v0 = np.array([0.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 20.0

#時間区間数
NT = 1000
skip = 10
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#インスタンス
rk4 = RK4(dt, r0, v0)
#グラフ描画用位置データ
tgs = []
rgs = []
ags = []

#解析解：初期位置のみ指定
def AnalyticalSolution( t ):
    omega = np.sqrt( k / m )
    return r0[2] *np.cos(omega * t)

#各時刻における運動方程式の計算
for tn in range(len(ts)):
    t = ts[tn]
    if( tn == 0 ):
        tgs.append( t ) 
        rgs.append( rk4.r[2] )
        ags.append( AnalyticalSolution( t ) )
        continue
    #ルンゲ・クッタ法による時間発展
    rk4.timeEvolution( t )
    #位置と速度を更新
    rk4.r += rk4.dr
    rk4.v += rk4.dv
    if( tn%skip == 0  ): 
        tgs.append( t ) 
        rgs.append( rk4.r[2] )
        ags.append( AnalyticalSolution( t ) )

plt.title( u"単振動運動（数値解と解析解）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(u"時刻" + r"$\,[{\rm s}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(u"位置" + r"$\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)
#スケールの指定
#ax = plt.gca()
#ax.set_yscale('log')
#ax.set_xscale('log')
#描画範囲を設定
plt.xlim([0, 20])
plt.ylim([-10.5, 10.5])
#プロット
plt.plot(tgs, ags, linestyle='solid', linewidth = 1 )
plt.plot(tgs, rgs, marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
#グラフの表示
plt.show()
投稿日:

【第８話】有限の高さの障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 有限の高さの障壁へ照射した平面波の時間発展（E = 1.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#虚数単位
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19

######################################
#　物理系の設定
######################################
#電子のエネルギー
E = 1.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E / hbar**2 )
#角振動数
omega = E/hbar
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 1000
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -5.0 * dx
x_max =  5.0 * dx
#描画区間数
NX = 500
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, x_max, NX)

#壁の高さ
V = 0.8 * eV
k2 = np.sqrt(np.complex128( 2.0 * me * (E - V) / hbar**2 ))

#アニメーション作成用
ims=[]
    
#各時刻における計算
for tn in range(NT):
    t = dt * tn

    #反射係数
    rc = ( k1 - k2 ) / ( k1 + k2 )
    #透過係数
    tc = 2.0 * k1 / ( k1 + k2 )

    #波動関数の計算
    phi_I = np.exp( I * k1 * x1 - I * omega * t )
    phi_R = rc * np.exp( - I * k1 * x1 - I * omega * t )
    phi_T = tc * np.exp( I * k2 * x2 - I * omega * t )
    phi = phi_I + phi_R

    #各コマを描画
    img  = plt.plot(x1/dx, phi_I.real, colors[0], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/dx, phi_R.real, colors[1], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/dx, phi.real, colors[2], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x2/dx, phi_T.real, colors[3], linestyle='solid', linewidth = 5.0)
    ims.append( img )

#グラフの描画（波動関数）
plt.title( u"有限の高さの障壁へ照射した平面波の波動関数（" + r"$ E = 1.0 [{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ {\rm Re}\{\psi(x, t)\} $", fontsize=30)

#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.12, top = 0.95)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#壁の描画
plt.vlines([0], -0.05, 1.05, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], -5,  0, "black", linestyles='solid', linewidth = 10 )
plt.hlines([1], 0,  5, "black", linestyles='solid', linewidth = 10 )


#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx])
plt.ylim([-2.1, 2.1])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)

#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()
投稿日:

【第７話】無限に高い障壁に向けた電子パルスの照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 無限に高い障壁に向けてガウス波束を照射したときの時間発展（E = 1.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#虚数単位
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19

#################################################################
## 物理系の設定
#################################################################
#波束の中心エネルギー
E0 = 10.0 * eV
#重ね合わせの数
NK = 200
#ガウス分布の幅
sigma = 2 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
#空間分割サイズ
dx = 1.0 * 10**-9
#空間分割数
NX = 500
#描画範囲
x_min = -10.0 * dx
x_max =  0.0 * dx

#ガウス波束の初期位置
x0 = -7.5 * dx

#計算時間の幅
ts = 0
te = 240
#時間間隔
dt = 0.5 * 10**-16

#ガウス波束の値を計算する関数
def Psi( x, t ):
    #波動関数値の初期化
    psi = x * (0.0 + 0.0j)
    #各波数ごとの寄与を足し合わせる
    for kn in range(NK):
        #各波数を取得
        k = k0 + dk * (kn - NK/2)
        #波数から各振動数を取得
        omega = hbar / (2.0 * me) * k**2
        #入射波と反射波
        phi_I = np.exp( I * k * x - I * omega * t )
        phi_R = - np.exp( - I * k * x - I * omega * t )
        #平面波を足し合わせる
        psi += (phi_I + phi_R) * np.exp( - I * k * x0 ) * np.exp( -( (k - k0) / (2.0 * sigma) )**2 )

    return psi

#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0 ))

#座標点配列の生成
x = np.linspace(x_min, x_max, NX)

#アニメーション作成用
ims=[]

### 波束の空間分布の計算
for tn in range(ts, te + 1):
    #実時間の
    t = dt * tn
    #時刻t=0の波動関数を取得
    psi = Psi( x, t )
    #波動関数の規格化
    psi /= psi_abs
    #各コマを描画
    img  = plt.plot(x/dx, psi.real, colors[0], linestyle='solid', linewidth = 3.0)
    img += plt.plot(x/dx, psi.imag, colors[1], linestyle='solid', linewidth = 3.0)
    img += plt.plot(x/dx, abs(psi),  colors[2], linestyle='solid', linewidth = 3.0)
    time = plt.text(8, 1.15, "t = " +  str(round(t/10**-15, 1)) + r"$[{\rm fs}]$" , ha='left', va='center', fontsize=18)
    #テキストをグラフに追加
    img.append(time)
    #アニメーションに追加
    ims.append(img)

#グラフの描画（波動関数）
plt.title( u"無限に高い障壁に向けて照射したガウス波束の波動関数（" + r"$ E = 10.0 [{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel( r"$ x [{\rm nm} ] $", fontsize=30 )
plt.ylabel( r"$ \psi(x,t) $", fontsize=30 )
#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.15, top = 0.95)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#壁の描画
plt.vlines([0], -5, 5, "black", linestyles='solid', linewidth = 10 )

#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx + 1])
plt.ylim([-1.1, 1.5])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
ani.save("output.mp4", writer="ffmpeg", dpi=300)

#################################################################
## ガウス分布の描画
#################################################################
fig_gaussian = plt.figure(figsize=(15, 8))

#波数の計算範囲
k_min = k0 - dk * NK / 2.0
k_max = k0 + dk * NK / 2.0

#波数座標点配列の生成
k = np.linspace(k_min, k_max, NK)
#正規分布
c_n = np.exp( -((k - k0) / (2.0 * sigma))**2 )

#グラフの描画（固有関数）
plt.title( u"ガウス分布（正規分布）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ (k_n - \bar{k})/\sigma $", fontsize=30)
plt.ylabel(r"$ c_n/c_0 $", fontsize=30)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([(k_min- k0)/sigma, (k_max- k0)/sigma])
plt.ylim([0, 1.05])

#波数分布グラフの描画
plt.plot((k-k0)/sigma, c_n, linestyle='solid', linewidth = 5)

#グラフの表示
plt.show()