T_NAKAの阿房ブログ

アクセスカウンタ

zoom RSS 高次モーメントを考える

<<   作成日時 : 2008/08/10 00:56   >>

ブログ気持玉 0 / トラックバック 0 / コメント 24

「基底状態の平均と分散を求めてみる」 http://teenaka.at.webry.info/200808/article_7.html で、分散レベルでは周期的境界条件で考えたモデルと同じであることを述べました。そこで「高次モーメントでは異なってくる」というご指摘を受けましたが、どうもそう思えないのです。これについてMaximaで確かめてみました。

まず、高次モーメントについての確認です。 積分範囲-∞→+∞ で、

1次モーメント; これは平均を表わします。(力学での重心に相当)(違う表現もありますが、本質的な違いは無いでしょう。)

μ1=∫yf(y)dy

2次モーメント; これは分散を表わします。(力学での重心周りの回転モーメントに相当)

μ2=∫(y-μ1)2f(y)dy

3次モーメント; これは分布の左右非対称の度合いを表わします。

μ3=∫(y-μ1)3f(y)dy

4次モーメント; これは分布の峰の鋭さ(裾野の広さ)を表わします。

μ4=∫(y-μ1)4f(y)dy

というところが、実用的な統計学で使う高次モーメントです。
もちろん、理屈ではもっと上の次数のモーメントも考えられ、数理統計学で使用されます。
つまり、一般的に m次モーメントは

μm=∫(y-μ1)mf(y)dy

で表わされます。

ここで分かることは、確率密度関数 f(y) は y=μ1 を中心に左右対称であるとすると、m が奇数ならば

(y-μ1)mf(y)

が奇関数になるので、μm はゼロですね。

さて、本題に入ると、

f(p)=(πα/2)[cos2(αp)/{(αp)2-(π/2)2}2]

ですが、θ≡αp 、さらに

g(θ)=(π/2)[cos2θ/{θ2-(π/2)2}2]

とすると、f(p)=αg(θ) なので、g(θ)の様相を見てみると、次のようになります。

画像


つまり、f(p)はp=0で左右対称であり、pの平均はゼロであることが分かります。
これから、必然的に μ1=0 から、

μ2n+1=∫p2n+1f(p)dp=0 (n=1,2,…)

つまり、奇数次モーメントはゼロであることが分かります。
では偶数次モーメントはどうなるのでしょう。 n=1,2,… として

μ2n=∫p2nf(p)dp=(1/α2n+1)∫(αp)2nf(p)d(αp)=(1/α2n+1)∫θ2nαg(θ)dθ=(1/α2n)∫θ2ng(θ)dθ

と計算されます。
とすると、

I2n=∫θ2ng(θ)dθ

という定積分(積分範囲-∞→+∞)の値が分かれば良いことになりますね。
Maximaを使って確かめましょう。(真ん中がコマンドで、右が結果です。)

n=1 → (%pi/2)*(integrate((x^2)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)); → %pi^2/4
n=2 → (%pi/2)*(integrate((x^4)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)); → %pi^4/16
n=3 → (%pi/2)*(integrate((x^6)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)); → %pi^6/64
n=4 → (%pi/2)*(integrate((x^8)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)); → %pi^8/256
n=5 → (%pi/2)*(integrate((x^10)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)); →%pi^10/1024

つまり、

I2=(π/2)2 ; I4=(π/2)4 ; I6=(π/2)6 ; I8=(π/2)8 ; I10=(π/2)10

なので、

I2n=(π/2)2n

ということになり、

μ2n=(1/α2n)(π/2)2n=(π/2α)2n

となります。ここで α=a/2ħ から、最終的に

μ2n=(πħ/a)2n

という結果が出ます。

これは、周期的境界条件で考えたモデルで計算したものと同じです。

テーマ

関連テーマ 一覧


月別リンク

ブログ気持玉

クリックして気持ちを伝えよう!
ログインしてクリックすれば、自分のブログへのリンクが付きます。
→ログインへ

トラックバック(0件)

タイトル (本文) ブログ名/日時

トラックバック用URL help


自分のブログにトラックバック記事作成(会員用) help

タイトル
本 文

コメント(24件)

内 容 ニックネーム/日時
T_NAKAさま 先にコメントしましたようにp→無限大での減衰の具合から発散と存じます。それらしい数値がでてくるのはプログラム中でフーリエ変換がなされ、つまり座標空間で、pは微分にかわり井戸の壁のところでの2階以上の微分のギャップ(または超関数)が考慮されていないのかもしれないと思いました。
参考 folo:fphys/290/topic/1/39 議論が乱歩のようで恐縮ですが...
甘泉法師
2008/08/10 11:46
甘泉法師さん

ご自分の説が正しいと主張されるのは自由ですので、主張されること自体に対して私は何も申し上げません。ただ、Maximaというツールの力を借りていますが、私はそれなりの時間を掛けて考察しているのです。ですから、私の行った方法に則った論理の何処かに誤ったところがあるなら、そこを指摘して下さい。「p→無限大での減衰の具合から発散と存じます」という言説は私の論理の前では、意味を成しません。これによって、私を折伏・説得することは無駄です。
この記事ではフーリエ変換後のp表示の確率密度関数について議論していることはお分かりで居られると思っておりましたが、「それらしい数値がでてくるのはプログラム中でフーリエ変換がなされ」などという言説が出てくることで、理解されていらっしゃらないことが分かります。"integrate"という命令はただ単に積分を表わしているだけですよ。「プログラム中でフーリエ変換」は実施されておりません。 「ランダウ=リフシッツ物理学小教程_量子力学」で示された答えを確率密度関数として高次モーメントを求めているだけです。
T_NAKA
2008/08/10 13:52
 どもTOSHIです。

 ここでの話は計算方法の話なので関係はないのでしょうが,folomyでの甘泉法師さんの参加していた議論でも感じたのと同じ感想を述べさせてもらいます。私自身も元の職業はプログラマーで積分はシンプソンやガウス積分で数値的にいくらでも正しく得られることを知っていますのでT_NAKAさんの方に理があると思います。

 そこでのはんどるさんと甘泉法師さんの議論でも計算技術の微細に終始していたのでおそらく「検算するとか,理論を確かめるだけ」という意味だろうと静観してりましたが,フーリエ変換という数学の話とは別に物理の定式化として座標運動量の両表示は等価なのでいかに高次の分散であろうと表示の違いで期待値は変わるはずがなく,変わっていたなら「計算間違い」かあるいは「考え違い」だろうと思っていました。
 今回もそうではないでしょうか。。   
                          TOSHI
TOSHI
2008/08/11 07:46
単純にt→∞での振る舞いが悪いのでこの積分が定義できているはずはないと思いますが、なぜそこを指摘されているのに調べないのでしょう。
マキシマの答えを見る限り、これは何らかの解析接続の式を使ったと思われます。例えば
< (θ^2-(Π/2)^2)^2>
= ∫dθ g(θ) (θ^2-(Π/2)^2)^2
= ∫dθ cos(θ)^2
を考えるとTANAKAさんの計算では、之は4次、2次、0次のモーメント線形結合ですが、これも収束することになりますがどうなのかな・・・。こういった式だと解析接続が使えないのでプログラムは発散を返してくると思われますがどうでしょう。

もっと直接的にはn次のモーメントを数値積分すればよいと思いますが、被積分の分母にポールがありますが、そこは多分マキシマが処理してくれると期待します。(見かけのポールの除去はスタンダード)例えば有限区間で積分してみると答えが返ってきて、∞区間だとエラーがでると期待しますが、どうでしょうか? マキシマは使ったことがないのでなんとも言えませんが、結果を是非教えてください。

φ男
2008/08/11 08:09
>単純にt→∞での振る舞いが悪いのでこの積分が定義できているはずはないと思いますが、なぜそこを指摘されているのに調べないのでしょう。

というのは言いすぎではないですか?私は誰彼から急かされて、何かを実行しなければならない立場の人間ではありません。
私を批判されるより、別の方法で調べて、結果を示して下さるのが礼儀ではないでしょうかねぇ。。
T_NAKA
2008/08/11 08:51
T_NAKAさん 甘泉法師です。
∫ x^2n cos^2 x / (x^2 - a^2)^2 dx > ∫ x^2(n-2) cos^2 x dx
 xを区間[nπ,(n+1)π]にわけcos^2 xの区間平均1/2を区間で最小のnπでの値に掛け 
> π/2 Σ(nπ)^2(n-2) = ∞ for n=4,6,8,...
甘泉法師
2008/08/11 10:49
訂正:1 / (x^2 - a^2)^2 > 1/ x^4 がいえる 区間|x|>a の積分と見てください。区間|x|<a の積分は有限値です。
甘泉法師
2008/08/11 11:02
訂正 ふたつの意味でnをつかってしまいました。区間のほうをmにかえて
> π/2 Σ(mπ)^2(n-2) = ∞ for n=4,6,8,...
    m
甘泉法師
2008/08/11 11:14
Maximaで (%pi/2)*(integrate((x^4)*(cos(x))^2,x,minf,inf)); と入力すると、結果は 0 となる事実があります。そうすると、

∫ x^4 cos^2 x / (x^2 - a^2)^2 dx > 0

というつまらない結果になります。よって、∫x^4 cos^2 x≠0 であり、Maximaが間違っているという理由を指摘して頂かないと水掛け論です。
これ以上議論をしても価値が無いと思います。
T_NAKA
2008/08/11 16:02
T_NAKAさん 甘泉法師です。
>(%pi/2)*(integrate((x^4)*(cos(x))^2,x,minf,inf)); と入力すると、結果は 0 となる事実があります
行の読み方に自信がないのですがいたるところ正の関数 x^4 cos^2 x を積分して0になるのということですか? わたしならコンピュータのほうを疑いますが… =甘泉法師= 
甘泉法師
2008/08/11 16:58
訂正 正の→非負の
甘泉法師
2008/08/11 17:19
高次モーメントについて、ただ予測しているだけで、何方も正確に計算して発散を証明していません。Maximaは答えを出しています。ちゃんと積分を解いてみせて下さい。分散が確定値になることも、正確に計算して示していただかないと私は納得しませんので、これでこの件に対するRESは辞めます。無責任な憶測だけでは時間の無駄です。
いい加減無駄な議論は止めましょう。
T_NAKA
2008/08/11 17:29
積分の発散を示すのは収束値を求めるよりずっと楽で、上のコメント(訂正2つを付けたもの)で示したつもりでした。もし興味を持たれてご批判いただければ対応いたします。
コンピュータがもっともらしい(失礼、これは拙考え。T_NAKAさんにとっては真の)収束値を返してくる仕組みに興味があります。
甘泉法師
2008/08/11 17:47
>コンピュータがもっともらしい

言いすぎです。いい加減にして下さい。
T_NAKA
2008/08/11 17:54
分散の計算は、複素積分で求められます。T_NAKAさんの求められた結果と一致します。ご興味があれば folo:fphys/290/topic/14/79 を見ていただくか、ご要望があればさわりを説明いたします。
=甘泉法師=
甘泉法師
2008/08/11 17:57
止めていただかないと、コメントはこちらで消すことになりますので、ご了承願います。
T_NAKA
2008/08/11 18:10
収束性が悪いというφ男さんの指摘は、もっともだと思います。
 数式処理上、定積分命令で∞範囲の場合、収束の良い他の形に式変形を行ってから数値積分している可能性があります。
下記ような∞範囲指定を、
integrate((x^10)*(cos(x))^2/((x^2-(%pi/2)^2)^2),x,minf,inf)
例えば、0〜10^10,0〜10^20,0〜10^30,…として、この処理系での数値積分の精度を評価することが必要だと思います。
 バグレポートに分母にポールがある関数での数値積分のことが記載されていたように思います。
はんどる
2008/08/11 21:27
だから、皆さんご自分で解決してくださいよ。指摘するだけでは、単なるシッタカですよ。
はんどるさんは、私が分散を計算したときに批判しましたよね。周期的境界条件で考えたモデルで計算したものと同じにならないはずだと仰いました。結果は同じです。ということはご自分では計算されていなかった訳ですよね。
T_NAKA
2008/08/11 21:42
> ということはご自分では計算されていなかった訳ですよね。
違います。
folo:fphys/290/topic/1/455
>>p表示で求めた運動量期待値の高次モーメントが発散するから矛盾等々の議論をしていたとき、数値計算を試みて確認しました。
(pcos(p))^2/(p^2-k^2)これを、自前スクリプト(float)で数値積分計算し一致しないことを確認しています。(収束性が悪い式そのもので、直接数値積分を強引に実行し、しかもfloat型で計算したのでかなり誤差があります。)

>結果は同じです。
 数式処理で内部的にフーリエ変換し収束性の良い形に直してから計算しているのなら、同じ結果になって当然ですね。

 はっきり言えば、使用した数式処理ソフトが同じ結果を与えるようにプログラムされているということがいえるだけで、この問題についての結論付けはできない。
はんどる
2008/08/12 06:13
抜けがありましたので訂正します。
(pcos(p))^2/(p^2-k^2)^2
 ps.数値計算での評価は難しいですね。有限桁で計算するので本質的に誤差が避けられません。
はんどる
2008/08/12 07:10
はんどるさんもお集まりなので、TOSHIさんの上コメントの問題意識にそって両者の認識の紹介です。
1 運動量表示での計算結果は共有。<p^2n>(n>2)は発散。
2 座標表示での計算と 一致する(甘泉法師)・しない(はんどるさん)
3 なぜならば波動関数の微分p^nψでは井戸境界に超関数 θ、δ、δ’,.. が
  ある(甘泉法師)・ない(はんどるさん)
4 よって運動量表示は座標表示と対等(甘泉法師)、座標表示からのフーリエ変換(運動量表示)は物理として問題(はんどるさん)
以上。
甘泉法師
2008/08/12 09:27
どうも皆さんは勘違いされているようです。たまたま、Maximaで答えを出せたので発表しただけです。それに反論されるなら、皆さん独自のブログなりを立ち上げてそこに書いて頂きたいのです。膨大なfolomy論議のどこそこで示したなどというのは、こちらに探す労力を要求することで、失礼なこととは思いませんか?重箱の隅を突くような内容は当事者以外は興味を持ちません。そこをご理解願います。特に「甘泉法師」さんと「はんどる」さんは1年以上もfolomyで論議されて決着されていません。場所をここに移して同じことをされてはたまりませんよ。ここは掲示板ではありませんので、占有が目立つようであると、こちらの判断で消します。つづきはfolomyでお願いします。
T_NAKA
2008/08/12 13:02
ご迷惑をおかけしました。おかげさまでMaximaというプログラムがあると知ることができました。ありがとうございました。
甘泉法師
2008/08/12 21:19
フリーですので、Maxima 公式サイト( http://www.cymric.jp/maxima/top.html )からDownLoadして試してみてください。
http://teenaka.at.webry.info/200707/article_30.html 辺りで関連記事を書いています。
T_NAKA
2008/08/12 22:19

コメントする help

ニックネーム
本 文
高次モーメントを考える T_NAKAの阿房ブログ/BIGLOBEウェブリブログ
文字サイズ:       閉じる