スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

マクスウェル分布 - 平均速さ編

今度はちゃんと統計力学でいうところの、マクスウェルの速度分布のグラフを描いてみました。横軸が速さ、縦軸が確率密度(分子の相対数)になっています。

速度分布式

見やすさだけを考えればやっぱり前回の確率密度関数のほうがいいかもしれませんね。でも物理的意味は断然こちらのほうがわかりやすいです。以下にRのスクリプトを示します。摂氏温度の値を色々変えてグラフを描いてみました。

# 摂氏温度(温度は任意に変えて下さい)
temp <- 25

# 横軸範囲
v = seq(0, 2500, .1)

# 希ガス原子
name = c("He =4", "Ne =20", "Ar  =40", "Kr  =84", "Xe =132")
weight = c(4,20,40,84,132)

######## 調整はここまで ########

#filename <- paste(temp,"_deg",".jpg", sep="")  # ファイル名の拡張子(jpg,pdfなど)
#jpeg(filename,width=600, height=650)           # px指定
#pdf(filename,width=6, height=6)                # in指定

# 確率密度関数
maxdis_true <- function(a,v){
  sqrt(2/pi)*(1/a^3)*v^2*exp(-v^2/(2*a^2))
}

abst = 273.15 + temp                # 絶対温度
k = 1.381e-23                       # ボルツマン定数
Na = 6.022e+23                      # アボガドロ定数

weight = weight/1000                # kgに直す
N = length(weight)
mm = max(weight)/Na                 # 重い分子が一番背が高い
ma = sqrt(k*abst/mm)
mp = max(maxdis_true(a=ma,v))*1000  # 縦軸の限界値

par(tck=.01, bty="l", lwd=2, mai = c(0.85, 0.8, 1.5, 0.35))
for (i in 1:N) {
 m = weight[i]/Na
 a = sqrt(k*abst/m)
 dist = maxdis_true(a,v) * 1000
 plot(v,dist,type="l", xlim=range(v),ylim=c(0,mp),
 xlab="Speed: v [m/s]", ylab="Probability Density: f(v) x1000 [s/m]", col = i)
 par(new=T)
}

title(main="Maxwell-Boltzmann Molecular Speed \n
Distribution for Noble Gases", line=3)
mtext(paste("at", temp, "deg."), line = 1)

legend(x= max(v)*0.6, y= mp*0.75, paste(name),
col = 1:N,lwd=2, merge = TRUE, bg='gray90')

#dev.off()

摂氏温度は変数tempにテキトウに代入するだけです。ファイル名はJPEGで出力すれば「(摂氏温度)_deg.jpg」のようになります。PDFが欲しければそのようにコメントアウトしてください。

Rの場合、JPEG出力はあまりキレイではないのが残念です。本当はPNGが欲しいところですが、上手くいきません。

100℃での分布 500℃での分布 1000℃での分布

温度が高くなるほど「山が潰れてくる」ことがわかります。確率密度関数に対して、全区間での積分は常に1です。縦軸だけ見ると全然1に達しないように見えますが、横軸が大きい値なのでスケール的には1になるわけです。

全区間の積分は1

確率密度関数を知らなくても、「全事象の各々の確率を全て足したら100%になる」ということが直感的には分かると思います。

確率・統計をある程度勉強した人ならば、絶対に知っているのが次の統計量です。

そこで、この速度分布に対しても平均を定義することが出来ます。確率密度関数に対して平均または期待値を次のように定義します。

平均の定義

要するに変数xを賭けてから積分すればいいわけです。マクスウェルの速度分布も確率密度関数ですから、平均の速さを次のように定義します。

平均速さ

この積分を計算するには、結構面倒な部分積分が必要になります。

絶対に途中で計算ミスを起こすので、いつまで経っても正解に到達しないこともあります。こっちはお遊びでマクスウェル分布を見ているので、そんなに時間を費やせるはずもありません。そこでMaximaを使って積分を計算してしまいましょう。

まず確率密度関数であることを確認するため、無限積分を求めてみます。パラメータaは正であることに注意してください。積分が1であることを確かめたら、今度は平均を求めることにします。Maximaに次のように打ち込みます。

f(v):=sqrt(2/%pi)*(1/a^3)*v^2*exp(-v^2/(2*a^2));
integrate(f(v), v, 0, inf);
positive;

g(v):=v*f(v)
integrate(g(v), v, 0, inf);

計算結果はつぎの画像のようになると思います。

Maxima計算結果

したがって求める平均速さはつぎのようになります。

平均速さの公式

この結果を用いて、もう少し色々な議論ができるのですが、今回はここまでにしておきましょう。

みねさん。大丈夫ですよね?


コメント

マクスウェル分布の真髄が分かったような気がします。
教科書ではやたら難しく書いてあるので取っ付きにくいですが、ブログ形式で書いてもらえると分かりやすくて良いです。

しかし平均速度の数式を導くのはかなり大変ですね。平均や分散は現在勉強しているところですが、積分の知識まで必要だとは私にとっては痛いです。

また、Rの導入についても詳しく書かれていたので、簡単に実行まで至ることができました。

最後にセクシー動画についてですが、あくまでセクシーであることを忘れないでください。

マクスウェル分布の真髄

このようなコメントを戴き、恐縮しております。数式を導くのが大変な場合でもコンピュータでどうにか解決してしまうのが真の理系です。

マクスウェル分布については
http://en.wikipedia.org/wiki/Maxwell_Speed_Distribution
が一番よくわかる解説ですかね。

学校でも先生自身がマクスウェル分布を軽く見てますからね。こうやって実際にグラフを描いてみたりして、直感的に学ぶことも大切です。

まいどっ!コメントありがとうございます。

いえいえ、おせっかいじゃあありませんよ。尊敬のまなざしで読ませてもらっています。(あなたを見ることはできないがi-229

マクスウェル分布 の記事を読ませていただきました。
・・・当然のことながらさっぱり意味が分かりませんが。(汗

話は変わりまして、Maximaって、便利なのですね。
一応私は数Ⅲの微積はやっていますが、こういう、実際に作ってみるのはやったことがありません(まぁ、当然ですが)
こういう授業が学校でほしいものです(ちなみに、プログラミング等の授業はありません。情報の授業は3時間ほどしましたが。そのときは、プレゼンのような発表のためのスライド作りをしました。)

そうそう、YouTubeの大きさを変えて、私の紹介文の下に動画をくっつけてみました。
本当に、こういういじるのがおもしろいです。

分布

>yuukaさん
結局、軽い分子ほどスピードが出て、温度が高いほどスピードが出るということです。それを確率密度という形で表現したのがマクスウェル分布ということです。
高校物理では一様な速さで運動していると考えますが、大学ではさまざまな速さが分布していると考えるわけです。

こんなプログラム初めてです!
PDFをいただきたいです。

コメントの投稿















管理者にだけ表示を許可

トラックバック

この記事のトラックバックURI
http://wavelet.blog58.fc2.com/tb.php/119-d1f3e93a

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。