１次元井戸型ポテンシャル（１）

明後日の量子物理学の授業は，１次元井戸型ポテンシャル（復習）の続きなので，数値計算で結果を確認するためのコードをMathematicaで書いてみた。

１次元ポテンシャル $V(x)$に質量$m$の電子が束縛されているとする。$x \le 0$ で$V(x)=\infty$ ，$0 \lt x \lt a $で，$V(x)=0$，$a \le x$で$V(x)=V_0$として [eV]単位で与える。ポテンシャルのレンジ $a$は [Å]単位とする。$2mc^2=10^6$[eV]，$\hbar c = 2000$ [eV・Å] と近似した。$r2= \frac{2mc^2 V_0 a^2}{(\hbar c)^2} =(\pi/2)^2$以上で束縛状態が存在する。

In[1]:= r2[V0_, a_] := 10^6*V0*a^2/(2000)^2
In[2]:= Clear[pa, qa]
In[3]:= sol1 = FindRoot[{p^2 + q^2 == r2[50, 2], q == -p/Tan[p]}, {p, 5}, {q, 3}]
Out[3]= {p -> 5.41164, q -> 4.55128}
In[4]:= {pa, qa} = {p, q} /. sol1
Out[4]= {5.41164, 4.55128}
In[5]:= Clear[A, B]
In[6]:= sol2 =NSolve[{A Sin[pa] == B Exp[-qa], A^2 Integrate[Sin[pa x]^2 , {x, 0, 1}] +
                                     B^2 Integrate[Exp[-2 qa x], {x, 1, Infinity}] == 1}, {A, B}]
Out[6]= {{A -> -1.28052, B -> 92.8589}, {A -> 1.28052, B -> -92.8589}}
In[7]:= {A, B} = {A, B} /. sol2[[1]]
Out[7]= {-1.28052, 92.8589}
In[8]:= Plot[A Sin[pa x] HeavisideTheta[-x + 1] +  B Exp[-qa x] HeavisideTheta[x - 1],
{x, 0, 3}, PlotRange -> {-2, 2}]

図　１次元井戸型ポテンシャルの波動関数

（１次元井戸型ポテンシャル（２）に続く）