Document - 松本充弘の個人的ページ

機械理工学・マイクロエンジニアリング・航空宇宙工学専攻
熱物理工学 2015（松本担当分）講義資料
松本充弘
第２回 2015/04/20
[email protected]
*レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで．
前回のレポート課題について
課題ブラウン運動の計算プログラムを実行し，「どのくらい遠くに到達したか」をステップ数の関数としてプロ
ットし，考察する．
○レポートから（東野純平さん）原点からの距離を
r = √𝑥 2 + 𝑦 2 とし，縦軸を r，横軸を時間（=step
数×t）でプロットした．（中略）時間を長くしてい
くにつれて原点から離れていくのが分かる。ただ
し、離れ方に規則性は感じられない。
同じプログラムを再度動かしても同じ結果が得ら
れた。RAND 関数は無作為な数字を毎回出力す
ると思っていたため、予想外の結果であった。
○レポートから（岡田行正さん）計算を N ステップ進めたときの粒子
の移動距離を算出する試行を 5 回行った。N を 1000 から
100000000 まで変化させたときの粒子の移動距離を表 1 に示す。
なお、5 回の試行における平均値も算出した。また、その結果を
両対数グラフとして図 1、2 にプロットした。
図 1 においてステップ数 N の増加に伴って移動距離も増加す
ることが確認できるが、試行ごとにかなり結果にばらつきが出て
くることもわかる。移動距離の N に対する依存性を調べるために、
図 1 ステップ数と移動距離のプロット
五回の試行における平均値を両対数でプロットしたのが図 2 で
ある。図 2 を見てみると、プロットがほぼ直線に乗っていることが
わかる。つまり、移動距離とステップ数は比例の関係にあると予
想できる。ステップ数は時間の経過に対応しているから、時間
経過に比例して粒子は拡散していくことがわかる。
松本のコメント：５回の試行の平均をとるなど，
「ランダムな現象の扱
い」について，正当な考え方だと評価できます．ただ，せっかく
両対数グラフ上に描いたのに「プロットがほぼ直線に乗っている
ことがわかる。つまり、移動距離とステップ数は比例の関係にあ
図 2 ステップ数と平均移動距離のプロット
る」というのはおかしいですね．
○レポートから（蛭谷摩周さん）両対数グラフで表してみると，比例に近い形のグラフとなったため，ステップ関
数を 1,10,100…と大きくすると，原点からの距離は指数関数的に大きくなっていくと考えられる．
熱物理工学講義資料
1
松本のコメント：うーん，
「指数関数」という言葉遣いは間違ってますね．
○レポートから（松岡諒さん） Fig. 1～3 より，計算ステ
ップ数を増やすとブラウン運動する粒子が原点か
ら遠くまで移動することが分かる．また，速度につ
いてステップ数を増やすにつれ，より円に近い軌
跡を描くことから，粒子の速さ√𝑢2 + 𝑣 2 が一定に
近付くと考えられる．
○レポートから（渋野雄一さん） 2 次元ランダムウォー
クでは無限時間後原点に戻るといいますが、案
外遠くへは行きませんでした。100000 回の試行
の場合、グラフの中央部に線が集中していました。
これは。元あった点に点が再び戻ってくること、ひ
いては 2 次元ランダムウォークの無限時間後の回
帰をも示唆しているのではないかと思います。
○レポートから（笠原史禎さん） x 軸にステップ数，y 軸に粒子が最終的に到達した距離をグラフにする（中
略）プロットした点を最小二乗法でフィッティングし，黒実線で示した．実線の式は近似的に以下の式で
表される．
y=
0.21x^(0.52) + 1.21
1.21 は誤差として大きく，グラフは原点を通るとすると，原
点からの粒子の距離はステップ数の約 1/2 に比例すること
がわかる．ステップ数は計算時間と比例するため，これを
言い換えると粒子の距離の２乗は時間に比例していること
がわかる．
全体的なコメントと補足
今日の講義で詳しく扱いますが，このように確率的に決まる問題について「出発点からの距離」を正
しく評価するには，平均２乗変位 mean square displacement (MSD) を求めるのが普通です．ブ
ラウン粒子の場合は，MSD は時間に比例することがちょっとした手計算でわかります．初回のレポー
トは，そのような背景をあえて説明することなしに，
「どのくらい離れたかを調べよ」という漠然とし
た課題にして皆さんに考えてもらいました．
統計平均をとる方法
1.
乱数の「種 seed」を変える：すぐに思いつく，ある意味で「正統な」方法
seed=こちらが指定する非負の整数;
srand(seed);
熱物理工学講義資料
2
2.
時間平均 time average をとる：長い計算（たとえば 108 ステップとか）を１回行った後，そ
れをいくつもの短いステップに分割すれば，平均をとることができる．分割されたそれぞれの動
きが独立と見なせることがその前提となる．
3.
集団平均 ensemble average をとる：たくさんの独立な粒子の動きを計算して，平均する．例
えば下のプログラムのように少し改造すればよい．
でも，よく考えると，実は，どれもほぼ同じことをやっていますね．
// ２次元ブラウン運動の軌跡を追跡する：ver. 2
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define
#define
#define
#define
#define
mass 1.0 // 粒子質量
gamma 1.0 // 摩擦係数
c_rand 1.0 // ランダム力の大きさ
delta_t 0.1 // 時間刻み
total_step 100000 // 計算ステップ数
10 万ステップ目の位置：10 粒子
double fr( )
{
double factor;
factor=sqrt((3*c_rand)/(2*delta_t));
return ( rand( )/(double)RAND_MAX*2 - 1.0 ) * factor;
}
int main( )
{
double x, y;
double u, v;
int step, repeat;
FILE *fout;
10 万ステップ目の位置：102 粒子
fout = fopen("brown-DifferentSeed.dat","w");
if (fout==NULL) {
printf("Error: cannot open file¥n");
return -1;
}
for (repeat=0;repeat<10000;repeat++) {
x=0.0; y=0.0;
u=0.0; v=0.0;
for (step=0; step<total_step; step++) {
x += delta_t * u;
y += delta_t * v;
u += delta_t/mass * (-gamma*u + fr( ));
v += delta_t/mass * (-gamma*v + fr( ));
}
printf("%4d %12.5f %12.5f %12.5f %12.5f¥n",repeat,x,y,u,v);
fprintf(fout,"%12.5f %12.5f %12.5f %12.5f¥n",x,y,u,v);
}
fclose(fout);
return 0;
10 万ステップ目の位置：103 粒子
10 万ステップ目の位置：104 粒子
}
熱物理工学講義資料
3
第２回レポート課題
出題日：4 月 20 日（月）
提出期限：4 月 25 日（土）次週の準備の都合上，できるだけ期限内に送って下さい．
提出方法：電子メールで，[email protected]
へ．
課題今日の講義で，速度自己相関関数 velocity autocorrelation function, VAF から輸送係数（拡散係数）
を求められることを示しました．このように，時系列データの自己相関を求めることは，いろいろな分野のデー
タ解析において基礎的な手法であり，身につけておくとどこかで役に立つかもしれません．そこで，ブラウン運
動について，速度自己相関関数 VAF(𝑡) を求めてみることにします：
VAF(|𝑡1 − 𝑡2 |) ≡ 〈𝑣⃗(𝑡1 ) ∙ 𝑣⃗(𝑡2 )〉
今回は，Ｃプログラミング初心者のために，途中まで作ったプログラムを示しますので，欠けている部分を補っ
て完成させ，VAF を求めてプロットし，簡単な考察と共に送って下さい．
（注意１） total_step 回（ここでは 10 万ステップとした）計算して統計平均をとるが，急速にゼロに減衰
することが予めわかっているので，相関関数は out_step （ここでは 1000 ステップ）までしか計算しな
い．
（注意２）講義中に触れたように，この系の時定数（減衰時間）は 𝛾⁄𝑚 程度と予想される．我々の単位系
（γ=1, m=1 ）では，これは１なので，計算時間刻み delta_t は 0.1 から 0.01 程度に小さくしない
と，詳細がわからない．
（注意３）プログラミングにある程度慣れている人は，このプログラムの構造が少し非効率であることに気付く
8
だろう．total_step が大きいと（たとえば 10 ステップとか）それが顕著に現れる．余力があれば，
その改善方法を考えてみよ．
// ２次元ブラウン運動の速度自己相関関数を求める
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define
#define
#define
#define
#define
#define
mass
1.0
gamma
1.0
c_rand 1.0
delta_t 0.01
total_step 100000
out_step 1000
double fr( )
{
double factor;
// 粒子質量
// 摩擦係数
// ランダム力の大きさ
// 時間刻み
// 計算ステップ数
// 出力ステップ数
// 0 以上 1 以下の一様乱数
factor=sqrt((3*c_rand)/(2*delta_t));
return ( rand( )/(double)RAND_MAX*2 - 1.0 ) * factor;
}
熱物理工学講義資料
4
int main( )
{
double x, y;
double u, v;
int step;
FILE *fout;
int i, j;
double dx, dy;
static int count[out_step];
static double vaf[out_step];
static double velx[total_step], vely[total_step];
x=0.0; y=0.0;
// 初期条件：原点から初速度ゼロで出発
u=0.0; v=0.0;
for (step=0; step<total_step; step++) {
x += delta_t * u;
y += delta_t * v;
u += delta_t/mass * (-gamma*u + fr( ));
v += delta_t/mass * (-gamma*v + fr( ));
velx[step]=u;
// 速度データ蓄積
vely[step]=v;
}
for (step=0; step<out_step; step++) {
count[step]=0;
vaf[step]=0;
}
for (i=0; i<total_step; i++) {
for (j=i; j<total_step; j++) {
step=j-i;
・
・
この間のプログラムを作る
・
・
}
}
// カウンタ初期化
// 速度自己相関関数の計算
fout = fopen("vaf.dat","w");
for (step=0; step<out_step; step++)
// 結果をファイルに出力
fprintf(fout, "%10.3f %10.3f¥n", delta_t*step, vaf[step]/count[step]);
}
fclose(fout);
return 0;
}
本日のレポート課題は以上です．次回は 4 月 27 日．
熱物理工学講義資料
5

Download Report