仮想計算機構

IT業界と無縁な派遣社員のブログ

デルタ関数と階段関数のシミュレーション【Sage】【Python】

デルタ関数

小出昭一郎『量子力学(Ⅰ)』の3章7節より


f(x)=\sqrt{\frac{\alpha}{2\pi}}\exp \{-\frac{\alpha}{2}(x-x_0)^2\}

で定義される関数の極限をとるとデルタ関数になる。


\delta (x-x_0)=\lim_{\alpha \to +\infty} f(x)

SageShellを立ち上げて、関数 f を下記で定義する。

x,a = var('x,a')
x0 = 1
f = lambda x,a:sqrt(a/2/pi)*exp(-a*(x-x0)*(x-x0)/2)

a を 1 から 50 まで 1 刻みで変化させ、長さ 50 のリストを作る。

sage: l2 = [plot(f(x,a),(x,-1,3),ymin=0,ymax=3)+line([(x0,0),(x0,3)],color="gray") for a in sxrange(1,51,1)]

a=1 の場合のグラフを見る場合は以下のようにする。

sage: l2[0]
Launched png viewer for Graphics object consisting of 2 graphics primitives

f:id:riverta1992:20201230155200p:plain

灰色部分がデルタ関数を表しています。a=1 なので関数 f とデルタ関数がまったく違う関数に見えます。

階段関数

デルタ関数、関数 f の積分を下記で定義します。


\eta(x-x_0)=\int_{-\infty}^{x}\delta (x^{\prime}-x_0)dx^{\prime}\\
\eta^{\prime}(x-x_0)=\int_{-1}^{x}f(x^{\prime})dx^{\prime}

a を 1 から50 まで変化させたときの関数 \eta^{\prime} の様子を可視化するためにサイズ50のリストを作ります。

sage: l = [plot(lambda r:numerical_integral(f(x,a),-1,r)[0],(r,-1,3),ymin=0,ymax=1.2)+line([(x0,0),(x0,1)],color="gray")+line([(x0,1),(3,1)],color="gray")+text("a={}".format(a),(2,0.5),fontsize=20,bounding_box={'boxstyle':'round','fc':'w'}) for a in sxrange(1,51,1)]

a=1 のときの関数 \eta^{\prime} のグラフを見る場合は下記のようにします。

sage: l[0]
Launched png viewer for Graphics object consisting of 4 graphics primitives

f:id:riverta1992:20201230155100p:plain

a=1 のときの関数 f と関数 \eta^{\prime} を並べて可視化するには

sage: graphics_array([[l2[0]],[l[0]]])
Launched png viewer for Graphics Array of size 2 x 1

とします。

f:id:riverta1992:20201230155315p:plain

最後に今まで作ってきたリストlとl2を使ってアニメーションを作ります。

sage: l3 = [graphics_array([[l2[i]],[l[i]]]) for i in range(len(l))]
sage: animate(l3)
Launched gif viewer for Animation with 50 frames

f:id:riverta1992:20201230154732g:plain