（
１次元井戸型ポテンシャル（１）からの続き）
テレビは朝から即位の礼のニュースで埋めつくされているのでなかなか気持ちが悪い。ラグビーワールドカップが終わった（実はまだ終わっていないの）と思ったらこれだ。オリンピックかＩＲカジノまでこの調子なのだろうか。
Mathematicaによる１次元井戸型ポテンシャルの解法をjuliaに移植してみた。Mathematicaプログラミングは土地勘があるので，簡単なガイドがあれば大丈夫だ。juliaプログラミングはそこまで熟達していないので，地図とガイドとネットでの評判を駆使して歩き回ることになる。ポイントは２つ。非線形方程式を解くFindRootや代数方程式を解くNSolveをどうするか。グラフをどうするか。それさえクリアすればよいのだが，なかなか難しかった。
非線形方程式を解くパッケージ NLsolve，１次元の数値積分を実行するパッケージQuadGKを導入した。図形描画のためのPlotとGRは既に導入済みである。こういうときに助けになるのが阪大のサイバーメディアセンターの
降旗大介さんのページ（Applied Mathematics 9）。NLSolveは生で使うと非常にわかりにくい仕様になっているので，降旗さんが
シンタックスシュガーを作ってくれている。おかげで比較的簡単に使うことができるが，MathematicaのFindRootの方がわかりやすいと思うのは気のせいか。規格化条件から波動関数の振幅を求める連立方程式も，MathematicaのNSolveに対応するものが見当たらなかったので，NLsolveを使うことにした。
問題はグラフだ。MathematicaのPlotルーチンにはなじんでいるので，およその様子はわかるが，juliaの方はさっぱりで難渋した。データを離散的なリストの形にするところまでは問題なかったが，そうすると横軸がデータ数でプロットされる。これをもとの変数に変換するためには，Plotの引数にxのリストを与える必要があることに気付くまで半日要した。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using NLsolve
using QuadGK
gr()
function nls(func, params...; ini = [0.0])
#
#スカラー変数 x スカラー関数 f(x,params)=0
# nls( f, params, ini = xの初期値 )
#ベクトル変数 x ベクトル関数 f(x,params)=0
# nls( f, params, ini = xの初期ベクトル )
#
if typeof(ini) <: Number
r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
v = r.zero[1]
else
r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
v = r.zero
end
# return v, r.f_converged
return v
end
function heaviside(t)
0.5 * (sign(t) + 1)
end
function r2(v0,a)
r=10^6*v0*a^2/(2000)^2
return r
end
function h1(x, p)
a,b = x # x[1] とか x[2] と書くのは面倒なので，a,b で代用
c,d = p # p[1] とか p[2] と書くのは面倒なので，c,d で代用
return [b+a/tan(a)+c, a^2+b^2-d]
end
function h2(y, q)
a,b = y # y[1] とか y[2] と書くのは面倒なので，a,b で代用
c,d = q # q[1] とか q[2] と書くのは面倒なので，c,d で代用
(f,hf) =quadgk(x -> sin(c*x)^2, 0, 1)
(g,hg) =quadgk(x -> exp(-2*d*x),1,10)
return [a*sin(c)-b*exp(-d), a^2*f+b^2*g-1]
end
function wf(x,s,t)
(pa,qa)=s
(a,b)=t
return [heaviside(1-t)*a*sin(pa*t)+heaviside(t-1)*b*exp(-qa*t) for t in x]
end
r = [0, r2(50, 2)]
ini_v = [2.0, 1.0]
s = nls(h1, r, ini = ini_v)
ini_u = [1.0, 1.0]
t = nls(h2, s, ini = ini_u)
x = 0:0.01:3
plot(x,wf(x,s,t))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
図 １次元井戸型ポテンシャルの波動関数
