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

機械理工学・マイクロエンジニアリング・航空宇宙工学専攻
熱物理工学 2016（松本担当分）講義資料
松本充弘
第３回 2016/06/13
[email protected]
*レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで．
前回のレポート課題について
課題ブラウン運動の速度自己相関関数 VAF(t) を求める：
𝑉𝐴𝐹(|𝑡𝑡 − 𝑡2 |) ≡ 〈𝑣⃗(𝑡1 ) ⃗⃗⃗⃗⃗(𝑡
∙ 𝑣 2 )〉
○レポートから（根岸知和さん）
ブラウン運動について速度自己相関数 VAF(𝑡) を計算し
た結果を図 1 に、 y 軸を対数とした片対数グラフおよび近
似直線を図 2 に示す。𝑡0≫𝑚/𝛾のとき、VAF(𝑡)は今回の単
位系（𝛾=1,𝑚=1,𝐶=1）では
VAF(𝑡) →
𝐶
2𝑚𝛾
𝛾
𝑒 −𝑚𝑡 = 0.5𝑒 −𝑡
(1)
と近似される。図１において、片対数グラフで直線となってい
図 1 二次元ブラウン運動の VAF
ると思われる範囲において近似直線を引くと
VAF(𝑡) = 0.479𝑒 −1.08𝑡
(2)
となった。これは式(1)で表されるものとほんど等しい。
図 1 を見ると、 t が 3 となる頃からグラフが直線から外れて
いることがわかる。この原因としては、ランダム力が完全に
ランダムではなく偏りを持った力が出されてしまうことが考え
られる。
松本のコメント：長時間での振る舞いが指数関数からはずれ
ていくのは，もちろんランダム力が完全でないという理
由もあり得ますが，多分，統計平均が不足していること
が最大の原因だろうと思います．それでも２桁近く小さ
くなる（つまり t=0 での値の 1/100）まで指数関数的振
図 2 二次元ブラウン運動の VAF 近似直線
る舞いが観測されたのなら，とりあえずは十分でしょう．
○レポートから（勝浦知也さん）
m = 1,2，gamma = 1,2 の計 4 通りについて計算を行っ
た．時定数は予想通りおおむね gamma/m となっている
ことが読み取れる．また，VAF の初期値は
1/(gamma*m) に比例していることも読み取れ，これら
の結果は式(2-59)のとおりになっていることが分かる．
プログラムの構造が非効率であることについては，速
度自己相関関数の計算のところで，for 文の２つ目が
熱物理工学講義資料
1
total_step まで計算するようになっているが，実際には i+out_step で十分であるため，そこで止めることで計
算量を減らすことができると考えられる．
○レポートから（白水仁さん）（前半省略）
・プログラムの非効率性： outstep を 1000 ステップまでしか計算しないのに対し、totalstep ばかり大きくする
と、無駄に記録する量が増えてしまいます。よって、今回の場合では 1001 ステップ以降の計算を省略する
と、プログラムが改善されると思います。
松本のコメント：計算を省略すると積算回数が減って統計平均の精度が悪くなりますね．計算はする
けど記録はしない，というのは，やり方によっては可能です．ただ，メモリへの記録にかかる「コ
スト」は，よほど totalstep が大きくない限り，一般にはそれほど大きくありません．
○レポートから（山本洋介さん）
一部を抜粋
求めた VAF は急速に減衰しているので，粒子の運動は気体的ではないといえる． VAF は振動している
ようにも見えることから，粒子の状態は固体に近いと予測できるが，ブラウン運動は液体中の微粒子が示
す運動なので，事実と異なる． VAF の振動は，シミュレーションの精度が原因かもしれない．
松本のコメント：そうですね，
「振動」は「復元力」の存在を示唆しますが，ここで扱っている Langevin
方程式には復元力はありませんから，
「振動」は統計平均
の不足によるものでしょう．
○レポートから（森村啄守さん）
自己相関数は，時刻がおよそ 2.8 および 7.8 のときに極
大となる。
松本のコメント：
「自己相関」について，どこかに誤解がありま
すね．
○レポートから（角木亮介さん）
一部を抜粋
プログラムの非効率性については，
「速度データの蓄積」と「VAF の計算」を別のループで行って
いる点が非効率であると考えた．そこで 2 つの計算を同じループで計算するプログラムに変更し
た．具体的には，
「カウンタ初期化」と「VAF の計算」を消去し，「速度データの蓄積」に下に示
すコードを追加した．効率化の結果をわかりやすくするために計算ステップを 1000000 回にして
検証したところ，計算時間は，元のプログラムでは約 17 秒，新しいプログラムでは約 13 秒とな
り，効率化できていることが確認できた．
追加したコード
for (i = step; i > 0 && step-i < out_step ; i--) {
//速度自己相関関数の計算
step2=step-i;
vaf[step2]=velx[step]*velx[i] + vely[step]*vely[i];
count[step2]++;
}
熱物理工学講義資料
2
○レポートから（伊藤圭人さん）
一部を抜粋
・コード改善： VAF を計算するループを以下のようにした。繰り返し回数を大きく減らすことが
できる。
for (step=0; step<out_step; step++) {
for (j=step; j<total_step; j++) {
vaf[step] += velx[j-step]*velx[j]+vely[j-step]*vely[j];
count[step]++;
} }
total_step=10^6 として clock 関数を用いてループ前後の経過時間を計測したところ、
改善前：1842.65 秒
改善後：5.77 秒
となった。
○レポートから（新井希さん）
一部を抜粋
dx, dy を定義した意味が分かりませんでした。
松本のコメント：おっと，多分，前回に変位を求めたときのコードの断片が残ってしまったのでしょ
うね．今回のコードでは明らかに不要でした．済みません．
補足：単位系について
質量：
基礎方程式に現れるパラメタの次元解析をすると
[𝑚] = M
摩擦係数：
[γ] = 力 / 速度 = M T-1
ランダム力の強さ：
[𝐶] = 力 2 / 時間 = M2 L2 T-5
となります．これから逆に，各基本次元を３つの量で表すと
M = [𝑚]
T = [𝑚⁄𝛾]
L = [𝑚3/2 𝛾 −5/2 𝐶 1/2 ]
が得られます．つまり計算プログラムの中で「質量，摩擦係数，ランダム力の大きさをすべて１にす
る」ということは，例えば「時間の単位を 𝑚⁄𝛾 とする」ということに相当します．
試しに，現実系の値を適当に見積もって代入してみましょう．
「花粉からの微粒子」の質量は詳しく
はわかりませんが，サイズとしてはきっと 10μm 程度のものでしょう．で，これが水に浮かぶのだか
ら，比重もだいたい１と考えると
𝑚 ～ 1 [g/cm3] ×4π/3×(10-5)3 [m3] ～ 4×10-12 [kg]
一方，Stokes 式 [講義ノートの(67)式] を利用すると，水の粘性係数(室温) η= 0.009 Pa・s より
Γ ～ 6π𝑎η ～ 20×10-5 [m] ×0.009 [Pa s] = 2×10-6 [kg/s]
よって，この計算系の実単位時間は 𝑚⁄𝛾 ～2×10-6 s = 2 マイクロ秒，ということになります．マイ
クロ秒程度で緩和する（昔の記憶を失う）からこそ，我々が観測すると微粒子はふらふらとさまよう
（＝酔歩）ように見えるのですね．
熱物理工学講義資料
3
第３回レポート課題
出題日：6 月 13 日（月）
提出期限：6 月 18 日（土）
提出方法：電子メールで，[email protected]
へ．
課題今日の授業で取り上げたスペクトル解析の練習です．講義ノート p.55 に例をあげたような「楽曲の
スペクトル解析」の真似をしてみよう．本格的におこなおうとすると，対象となるデータを準備するのが面倒な
ので，ここでは，次のようにごく簡単な「遊び」をやってみる：
(1) 音の強弱は無視して，音程のみを考える．
(2) 長い音は，同じ音程の音符が続いていると見なす．
(3) 休符も無視する．音の強弱を考える場合は，休符は「強さゼロ」と見なせばいいのだが，今回は強弱を
考えないため，休符の扱いに困るので…
ごく簡単な例として，「きらきら星 Twinkle, twinkle, little star」（もとはフランス民謡だそうです）を考え
てみよう：
メロディを数値化するのにも，いろいろな考え方がある．例えば
（方法１）ドレミファ…を，1, 2, 3, 4, … に置き換える．この場合，上の例なら
1 1 5 5 6 6 5 5 4 4 3 3 2 2 1 1 5 5 4 4 3 3 2 2 5 5 4 4 3 3 2 2 1…
となる．この場合，シャープ ♯ やフラット ♭ が出てきたら，1.5 などとすることになる．
（方法２）１２音階主義で，ドレミファ…を，1, 3, 5, 6 …に置き換える．上の例なら
1 1 8 8 10 10 8 8 6 6 5 5 3 3 1 1 8 8 6 6 5 5 3 3 8 8 6 6 5 5 3 3 …
となる．
ひんぱん
ここでは単なる遊びなのでどちらでもよいだろう．（シャープやフラットが頻繁に出てくるメロディーを選んだ場
合は，方法２のほうがいいだろう．）以下の計算例は（方法２）で行った．
このような数値データファイルを，Fourier 変換してパワースペクトルを求める．この程度のデータ量なら
FFT (Fast Fourier Transform) などを持ち出すまでもなく，次のような単純なプログラムで十分である：
#include <stdio.h>
#include <math.h>
//三角関数を使うのに数学ライブラリのリンクが必要
#define MAX_DATA 10000
//データ個数の上限
#define PI (4.0*atan(1.0))
//円周率 3.141592…を精度よく与える１つの方法， = arctan(1.0)
π
4
int main()
{
int ndata,i,j;
熱物理工学講義資料
4
double indata[MAX_DATA];
double omega,iomega,cdata,sdata,power;
FILE *fin, *fout;
fin =fopen("star.dat","r");
//読み込むデータファイルを指定
if (fin==NULL) {
printf("No INPUT File¥n");
return -1;}
fout=fopen("power.dat","w");
//結果を出力するデータファイルを指定
ndata=0;
//
↓ 読み込むデータがなくなるまで読む
while ((ndata<MAX_DATA) && (fscanf(fin,"%lf",&indata[ndata])==1)) ndata++;
fclose(fin);
printf("Number of data =%5d¥n",ndata);
omega=2.0*PI/ndata;
//離散 Fourier 変換の基本角振動数の定義
for (i=0;i<ndata;i++) {
//角振動数についてのループ
cdata=0.0; sdata=0.0;
//積算する変数をゼロに初期化
iomega=omega*i;
for (j=0;j<ndata;j++) {
cdata += indata[j]*cos(iomega*j);
//入力データについてのループ
//素朴に sin, cos を掛けて足し合わせる
sdata += indata[j]*sin(iomega*j);
}
power=cdata*cdata+sdata*sdata;
//パワー (= 絶対値の２乗) を求める
fprintf(fout,"%5d %10.4f¥n",i,power);
}
fclose(fout);
return 0;
}
このようにして得られたパワースペクトルの例が右図
である．「きらきら星」のスペクトルはおおよそ 1⁄𝑓 2 に
近いように見える．
では，以上をまねて，何か適当な音楽（メロディの一部でよい）について，パワースペクトルを求めてグラ
フを描いてみてください．
本日のレポート課題は以上です．次回は６月２０日（月）．
熱物理工学講義資料
5

Download Report