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

（１次元井戸型ポテンシャル（１）からの続き）

テレビは朝から即位の礼のニュースで埋めつくされているのでなかなか気持ちが悪い。ラグビーワールドカップが終わった（実はまだ終わっていないの）と思ったらこれだ。オリンピックかＩＲカジノまでこの調子なのだろうか。

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))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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