プログラミング言語「Python」には, 数値計算ライブラリ NumPy のためのグラフ描画ライブラリ Matplotlib が提供されているので, データの可視化を容易に行うことが出来るようだ.そこで, Python によって, 調和振動子の波動関数などのグラフ化やアニメーション化をしてみよう.しかし実際に Python でプログラミングする前に, まずは Bernd Thaller:「Visual Quantum Mechanics」の第7章から必要な事柄を抜粋して, 波動関数をグラフ化する際に必要な手続きを確認しておくことにする.
7.1.3. ハミルトニアンのスケーリング変換
グラフの表示や数学的な検討には, 長さ, 時間, エネルギーの単位を変えて, 式を簡略化することが有効である.適切な「スケーリング変換」(scaling transformation)を用いれば [1]「スケーリング」(scaling):(1) 系の振る舞いが物理量の絶対値だけでなく,相対的な大きさのみで決まること.特に2変数関数
2次の微分 (運動エネルギー演算子) の場合は次となる:
スケーリング変換に於いて, 位置演算子
従って, スケーリングされた関数に対するシュレディンガー演算子の作用は, 次の様に記述される:
7.1.4. 無次元の単位系(Orders of magnitude)
式 (7.15) 中の2つの加数が相反するスケーリングの振舞いを行うことから, 次式が成り立つようなスケーリングパラメータ
すなわち,
ただし
この式が意味するのは,「変位は長さ
式 (7.7) のシュレディンガー方程式は, 時間もスケーリングして次式
が成り立つように定義する, すなわち
と定義するならば, さらに簡略することが出来る.すると新たな単位による最終的なシュレディンガー方程式は, (
更には, エネルギーのスケーリングを再定義し
とするならば, 新たな単位で測定された調和振動子のエネルギーに対する演算子, すなわちハミルトニアン演算子は次となる:
新たな変数
7.1.5. 大きさのオーダー
ここで出会う物理現象に於いて期待される桁数を推定するには, 物理定数に実際の値を入れてみることが有効である.質量
すると, 時間, 長さ, そしてエネルギーは次のように設定されることになる:
今後は, この無次元の単位を用いる.特に図を見るときには, このことを忘れないようにすべきである.バネ定数が
水素原子のボーア半径(
光の速さ (
7.2.2. 固有関数への展開
以上を踏まえると, 調和振動子のエネルギー固有関数
は,
ただし,
全ての2乗可積分の波動関数は, 調和振動子の固有関数
7.3.1. 時間発展
次に系の時間発展を考える.そのためには, 与えられた初期関数
任意状態
初期状態
任意状態
調和振動子波動関数の Python によるグラフ化とアニメーション化
以上の事柄を踏まえて, Python により式 (7.34) の調和振動子の固有関数
このグラフ化の Python ソースコードを以下に示しておく [2]ただし, Google Colaboratory で実行する場合, このままではフォントエラーが発生し日本語がきちんと表示されない.その場合は, … Continue reading .
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 | import numpy as np import math import matplotlib.pyplot as plt # Google Colaboratory で実行する場合は, 下行先頭の#を外すこと. #import japanize_matplotlib % matplotlib inline amp = math. pow (math.pi, - 1 / 4 ) x = np.arange( - 5 , 5 , 0.01 ) y = amp * np.exp( - 0.5 * x * x) y1 = y y2 = math.sqrt( 2 ) * x * y y3 = ( 2 * x * x - 1 ) / math.sqrt( 2 ) * y y4 = ( 2 * x * * 3 - 3 * x) / math.sqrt( 3 ) * y y5 = ( 4 * x * * 5 - 20 * x * * 3 + 15 * x) / math.sqrt( 60 ) * y y6 = ( 1024 * x * * 10 - 23040 * x * * 8 + 161280 * x * * 6 - 403200 * x * * 4 + 302400 * x * * 2 - 30240 ) / math.sqrt(math. pow ( 2 , 10 ) * math.factorial( 10 )) * y vx = 1 / 2 * x * x fig,axes = plt.subplots( 1 , 2 , figsize = ( 13 , 7 )) axes[ 0 ].plot(x, vx, linewidth = 1.5 , color = "grey" , linestyle = "dashed" ) axes[ 0 ].plot(x, y1 + 1 / 2 , linewidth = 1.5 , color = "red" , label = "n = 0" ) axes[ 0 ].plot(x, y2 + 3 / 2 , linewidth = 1.5 , color = "blue" , label = "n = 1" ) axes[ 0 ].plot(x, y3 + 5 / 2 , linewidth = 1.5 , color = "green" , label = "n = 2" ) axes[ 0 ].plot(x, y4 + 7 / 2 , linewidth = 1.5 , color = "cyan" , label = "n = 3" ) axes[ 0 ].plot(x, y5 + 9 / 2 , linewidth = 1.5 , color = "magenta" , label = "n = 4" ) axes[ 1 ].plot(x, vx, linewidth = 1.5 , color = "grey" , linestyle = "dashed" ) axes[ 1 ].plot(x, y1 * y1 + 1 / 2 , linewidth = 1.5 , color = "red" , label = "n = 0" ) axes[ 1 ].plot(x, y2 * y2 + 3 / 2 , linewidth = 1.5 , color = "blue" , label = "n = 1" ) axes[ 1 ].plot(x, y3 * y3 + 5 / 2 , linewidth = 1.5 , color = "green" , label = "n = 2" ) axes[ 1 ].plot(x, y4 * y4 + 7 / 2 , linewidth = 1.5 , color = "cyan" , label = "n = 3" ) axes[ 1 ].plot(x, y5 * y5 + 9 / 2 , linewidth = 1.5 , color = "magenta" , label = "n = 4" ) # グラフのタイトル axes[ 0 ].set_title(r "調和振動子のエネルギー固有関数 $\psi_n(x)$" ) axes[ 1 ].set_title(r "エネルギー固有関数の絶対値2乗 $|\psi_n(x)|^{2}$" ) # X軸の範囲の設定 axes[ 0 ].set_xlim( - 5 , 5 ) axes[ 1 ].set_xlim( - 5 , 5 ) # X軸のラベル axes[ 0 ].set_xlabel( "x" ) axes[ 0 ].set_ylabel( "$E_n$ , $\psi(x)$" ) axes[ 0 ].set_ylim( - 0.5 , 6.0 ) axes[ 1 ].set_xlabel( "x" ) axes[ 1 ].set_ylabel( "$E_n$ , $|\psi(x)|^{2}$" ) axes[ 1 ].set_ylim( - 0.5 , 6.0 ) # 補助目盛りの描画 axes[ 0 ].minorticks_on() axes[ 1 ].minorticks_on() # 主目盛の設定 axes[ 0 ].set_xticks(np.linspace( - 5 , 5 , 5 )) axes[ 0 ].set_yticks(np.linspace( 0.5 , 4.5 , 5 )) axes[ 0 ].set_yticklabels([ "E0" , "E1" , "E2" , "E3" , "E4" ]) axes[ 1 ].set_xticks(np.linspace( - 5 , 5 , 5 )) axes[ 1 ].set_yticks(np.linspace( 0.5 , 4.5 , 5 )) axes[ 1 ].set_yticklabels([ "E0" , "E1" , "E2" , "E3" , "E4" ]) # 目盛り線の描画 axes[ 0 ].grid() axes[ 1 ].grid() # 凡例表示するには #axes[0].legend(loc="upper right", bbox_to_anchor=(1,1)) #axes[1].legend(loc="upper right", bbox_to_anchor=(1,1)) # サブプロット間スペースの調整 fig.subplots_adjust(wspace = 0.3 ) plt.show() |
また, 式 (7.49) の調和振動子固有状態, 及び式 (7.51) の任意関数として
の時間発展の様子をシミュレーションする Python アニメーションは例えば次のようである:
このアニメーションを見ると次のことが分かると思う:
「調和振動子のエネルギー固有関数は, 定常状態に相当するので, その時間変化は定常波となる.しかし, 幾つかのエネルギー固有関数を重ね合わせたものは, もはや定常状態の波動関数ではなく, その時間変化は進行波となる」.
このアニメーション化の Python ソースコードを以下に示しておく.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | import numpy as np import math import matplotlib.pyplot as plt # Google Colaboratory で実行する場合は, 次の行の先頭の # を外すこと. #import japanize_matplotlib from matplotlib import animation, rc from IPython.display import HTML fig, axes = plt.subplots( 1 , 2 , figsize = ( 10 , 5 ), constrained_layout = True ) # 各軸のラベル # --> rを入れないとLaTeX出力されないので注意する. axes[ 0 ].set_title(r "調和振動子のエネルギー固有関数 $\psi_n(x)$ の時間変化" ) axes[ 1 ].set_title(r "$\psi(x)=(\psi_0(x)+\psi_1(x)+\psi_2(x))/\sqrt{3}$と$|\psi(x)|^{2}$の時間変化" ) axes[ 0 ].set_xlabel(r "$x$" , fontsize = 15 ) axes[ 0 ].set_ylabel(r "$\psi_0(x),\quad\psi_1(x),\quad\psi_2(x),\quad E_n$" , fontsize = 15 ) axes[ 1 ].set_xlabel(r "$x$" , fontsize = 15 ) axes[ 1 ].set_ylabel(r "$\psi(x),\quad |\psi(x)|^{2}$" , fontsize = 15 ) # グラフの範囲を設定 axes[ 0 ].set_xlim([ - 4 , 4 ]) axes[ 0 ].set_ylim([ - 0.5 , 3.5 ]) axes[ 1 ].set_xlim([ - 4 , 4 ]) axes[ 1 ].set_ylim([ 1 , 5 ]) axes[ 0 ].set_xticks(np.linspace( - 4 , 4 , 9 )) axes[ 1 ].set_xticks(np.linspace( - 4 , 4 , 9 )) axes[ 0 ].set_yticks(np.linspace( 0.5 , 2.5 , 3 )) axes[ 0 ].set_yticklabels([r "$\hbar\omega/2$" , r "$3\hbar\omega/2$" , r "$5\hbar\omega/2$" ]) axes[ 1 ].set_yticks(np.linspace( 2.5 , 4.5 , 3 )) axes[ 1 ].set_yticklabels([r "$5\hbar\omega/2$" , r "$7\hbar\omega/2$" , r "$9\hbar\omega/2$" ]) axes[ 0 ].grid() axes[ 1 ].grid() #axes[1].legend() # x点の計算範囲とその点数 x = np.linspace( - 5 , 5 , 80 ) amp = math. pow (math.pi, - 1 / 4 ) T = 2 * math.pi dt = 1 / T vx = 1 / 2 * x * x axes[ 0 ].plot(x, vx, color = 'gray' , linestyle = "dashed" ) axes[ 1 ].plot(x, vx, color = 'gray' , linestyle = "dashed" ) ims = [] for a in range ( 200 ): y1 = amp * np.exp( - 0.5 * x * x) * np.cos( 1 / 2 * dt * a) y2 = amp * np.exp( - 0.5 * x * x) * math.sqrt( 2 ) * x * np.cos( 3 / 2 * dt * a) y3 = amp * np.exp( - 0.5 * x * x) * ( 2 * x * x - 1 ) / math.sqrt( 2 ) * np.cos( 5 / 2 * dt * a) ytotal = ( 1 / math.sqrt( 3 )) * amp * np.exp( - 0.5 * x * x) * (np.cos( 1 / 2 * dt * a)\ + math.sqrt( 2 ) * x * np.cos( 3 / 2 * dt * a) + ( 2 * x * x - 1 ) / math.sqrt( 2 ) * np.cos( 5 / 2 * dt * a)) im1 = axes[ 0 ].plot(x, y1 + 0.5 , color = 'red' ) im2 = axes[ 0 ].plot(x, y2 + 1.5 , color = 'green' ) im3 = axes[ 0 ].plot(x, y3 + 2.5 , color = 'blue' ) im4 = axes[ 1 ].plot(x, ytotal + 3.5 , alpha = 1.0 , color = 'brown' , label = r "\psi(x)=\psi_0(x)+\psi_1(x)+\psi_2(x)" ) im5 = axes[ 1 ].plot(x, ytotal * ytotal + 3.5 , linewidth = 8 , alpha = 0.5 , color = 'magenta' , label = r "|\psi(x)|^{2}" ) ims.append(im1 + im2 + im3 + im4 + im5) # ArtistAnimationにfigオブジェクトとimsを代入してアニメーションを作成 anim = animation.ArtistAnimation(fig, ims, interval = 70 ) # interval値でグラフの描画速度が決まる. #anim.save('anim1.mp4') # Jupyter notebookの場合は #plt.show() # Google Colaboratoryの場合必要 rc( 'animation' , html = 'jshtml' ) plt.close() anim |
References
↑1 | 「スケーリング」(scaling):(1) 系の振る舞いが物理量の絶対値だけでなく,相対的な大きさのみで決まること.特に2変数関数 | ||
---|---|---|---|
↑2 | ただし, Google Colaboratory で実行する場合, このままではフォントエラーが発生し日本語がきちんと表示されない.その場合は, プログラム実行前にまず, 次を実行する:
その後, プログラム中の5行目の先頭の#を消してその行を付加すること. アニメーションプログラムの方も同じ作業が必要である. |