1次元調和振動子の可視化 その2

はじめに

前回のつづきです。

今回も数値計算っぽいことをやります。

小出昭一郎『量子力学』のp.53より の期待値は

となります。前回と同様 とし、これを数値計算で確かめます。

方法

から をサンプリングして、 の平均値をとれば が計算できます。 からサンプリングする方法については下記サイトを参考にしました。

任意の確率密度分布に従う乱数の生成（von Neumannの棄却法） | Pacocat's Life

プログラム再掲します。

def X(n,M,w):
    A = sqrt(sqrt(w*M/pi)/2^n/factorial(n))

    return lambda x : A*hermite(n,x)*exp(-x^2*M*w/2)

def X2(n,M,w):
    A = sqrt(sqrt(w*M/pi)/2^n/factorial(n))

    return lambda x : (A*hermite(n,x)*exp(-x^2*M*w/2))^2

サンプリング用の関数も作っておきます。

import numpy as np

def sample(f,xmax,xmin,ymax,n):
    xs = []
    ys = []
    x2s = []
    while True:
        x = xmin + (xmax - xmin) * np.random.random()
        y = ymax * np.random.random()
        if f(x)>=y:
            xs.append(x)
            ys.append(y)
            x2s.append(x*x)
        if len(xs)>=n:
            break
    return {"x":xs,"y":ys,"x2":x2s}

今回は の場合を考えます。 のグラフは以下のとおりです。

sage: plot(X2(12,1,1),(-10,10))
Launched png viewer for Graphics object consisting of 1 graphics primitive

実際にサンプリングできるのか試してみます。

sage: histogram(sample(X2(12,1,1),10,-10,0.6,2000)["x"], bins=100,density=True)+plot(X2(12,1,1),(-10,10))
Launched png viewer for Graphics object consisting of 2 graphics primitives

グラフとヒストグラムがちゃんと重なっているので大丈夫そうです。

の場合、 になるはずなので、計算してみましょう。

sage: np.mean(sample(X2(12,1,1),10,-10,0.6,2000)["x2"])
12.347786052776227

いい感じの値が出ました。

参考

LaTeXコマンド集 - 上線・下線 (overline,underline)
Histograms — Sage 9.2 Reference Manual: 2D Graphics