熱物理工学 2016（松本担当分）講義資料

機械理工学・マイクロエンジニアリング・航空宇宙工学専攻
熱物理工学 2016（松本担当分）講義資料
松本充弘
第５回 2016/06/27
[email protected]
*レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで．
前回のレポート課題について
課題２値化画像を使って，フラクタル解析をおこなう．
サンプルプログラムの一部
Image[y][x]=idata;
...
// 黒 =0, 白 !=0
for (x = 0; x < Width/resolution; x++) {
for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく
・・・
この間のコードを作成する
・・・
} }
○レポートから（今福壮貴さん）図より，タイルの面積が小さいときは D=1 に近く，タイルの面積が大きくなると
D=2 によく一致している．タイルの面積が小さいときは 1
次元的な発想で解析しており，面積が大きいときは 2 次
元的で解析しているので妥当な結果だと考えられる．ま
た，タイル面積を大きくするとほぼ D=2 と一致しているこ
とから，この図形はフラクタル図形ではないと考えられる．
また，コードの空白部分を以下に示す．
for (x = 0; x < Width/resolution; x++) {
for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく
plus = 0;
multiply = 1;
for (xi = x*resolution; xi < (x + 1)*resolution; xi++){
for (yi = y*resolution; yi < (y + 1)*resolution; yi++){
plus += Image[yi][xi];
multiply *= Image[yi][xi];
}}
if (plus != 0 && multiply == 0){
sum++;
}
}}
松本のコメント：なるほど，Image=0 または ≠0 というのを利用して，積から判断するのは，うまい手ですね．
○レポートから（古田歩さん）
1000000
100000
次元は 2.012 となった．穴を埋めたコードは以下の通りです．
for (xi = x*resolution; xi<(x+1)*resolution; xi++) {
for (yi = y*resolution; yi<(y+1)*resolution; yi++) {
if (Image[yi][xi] !=0) {
sum++;
xi=(x+1)*resolution;
yi=(y+1)*resolution;
}
}
}
10000
系列1
1000
累乗 (系列1)
100
10
y = 714210x-2.012
1
1
10
100
1000
松本のコメント：タイル内に「白」のピクセルが見つかったら，xi と yi をタイルの上限にリセットすることで２重
ループを抜け出すという，なかなかよくできた（技巧的な）コードです．
熱物理工学講義資料
1
○レポートから（柴田悠希さん）
タイルサイズと図を覆うのに必要なタイルの枚数との関係を表
すグラフを fig1 に示す．ここで，レジュメのコードでは resolution
は 2 から始まっていたが，今回は 1 から解析した．また，R で割
り切れず余った領域に関しては無視した．
R が大きすぎない範囲ではおよそ R^(-1.9)と一致しており，フラ
クタル次元は 1.9 程度であると考えられる． R が大きくなると傾
きも大きくなるが，これはタイルのサイズがかなり大きいため，図
を覆うのに全てのタイルが必要となる領域となったからだと考えられる．（R でちょうど割り切れるサイズの
画像なら，次元は 2 に近づくはず．ここにあまりを無視した影響が加わっている．）
松本のコメント：「余った領域」については私も深くは考えていませんでしたが，確かに影響はありそうですね．
本来なら，最初から縦・横とも，２のべき乗（28 = 256 とか 29 = 512 とか 210 = 1024 とか）に揃えて
切り落としておくといいのですが，この出題では，そこをサボりました．
○レポートから（馬見新彩さん）
一部を抜粋
解析の結果，がん組織画像のフラクタル次元は 1.049 となった．
2 次元よりははるかに 1 次元に近く，縮小拡大にほとんど反比
例する形で面積が変化することがわかった．枝分かれした細部
は肉眼では確認することができないが，直感的に見積もったも
のよりもはるかに隙間が多いという結果になった．
松本のコメント：はて，多くの人が D=2 前後の値を報告しているの
に対して，D～1 は少数派ですね．．．
○レポートから（西岡寿朗さん）
一部を抜粋
（前半省略）コンパイラは visualc++ 2008 express edition を用いた．プログラムを実行した
ところ…タイルの総数は常に 0 であった．デバッグを行ったところ，Image[ ][ ] の中身は 0 と
それ以外の二値ではなく，-872415232 と-855638017 の
二値であることが確認された．試しにコンパイラとして
Code::Blocks を用いたところ，Image[ ][ ]は 0 と
16777215 の二値であった．コンパイラにより
Image[ ][ ] の値が変わることが確認された．この理由
を考えることは本課題の範疇を超えると判断した．
-872415232 と-855638017 のどちらが黒(0) に対応し
ているか不明であったのでプログラム中の==0 を
==-872415232，==-855638017 に変更し両方のフラクタ
ル解析を行う．
（以下省略）
熱物理工学講義資料
2
松本のコメント：申し訳ありません．手元のコンパイラ（visual studio 2013）では問題が現れなかったので，全く
気付いていませんでした．手元でこの現象が再現できないので，100%の確信は持てないのですが，おそ
らく，元のプログラムで危険な部分は，データを読み込むところです．
// Read Image
for (y=0;y<Height;y++) {
for (x=0;x<Width;x++) {
if (fread(&idata, 3, 1, fin)<1) {
printf("Error: data end while reading¥n");
return -1;
}
Image[y][x]=idata;
// 黒 =0, 白 !=0
}}
ちょっと専門的な話になりますが，fread(&idata, 3, 1, fin) は，ファイル fin から，3 バイトのデータ１つを
idata に読み込みます．問題は，idata が unsigned int という型をもち，通常は 4 バイトのサイズであるの
に対して，3 バイトしか読み込まないという点です．つまり，残りの 1 バイトは読み込み時に上書きされず
に残ってしまう，ということです．3 バイトというのは普通の変数のサイズではないのですが，ここで扱う
bitmap ファイルが，１ピクセルを 24 bit (=3 バイト) で表現しているのでこうなっています．残りの１バイト
に何が入っているかによって，4 バイト変数の Image[y][x]に代入するときの idata が変わってしまったのだ
ろうというのが私の推測です．たとえば
// Read Image
for (y=0;y<Height;y++) {
for (x=0;x<Width;x++) {
idata=(unsigned int)0;
if (fread(&idata, 3, 1, fin)<1) {
printf("Error: data end while reading¥n");
return -1;
}
Image[y][x]=idata;
// 黒 =0, 白 !=0
}}
のように，読み込む前に idata を初期化してやると，このような現象がなくなるのではないかと思いますが，
手元で実験できないので・・・来年，似たような出題をするときは，注意します．ご指摘をありがとうござ
いました．
＊ここから先，もっと maniac な話：ちなみに，謎の数字 -872415232 と -855638017 ですが，差をとると
−855638017 – (−872415232) = 16777215 = 224 − 1
ですので，多分，もともと idata に入っていた「ゴミ」の下部２４ビットを初期化したのが -872415232 で，
下部ビットをすべて 1 で埋めたのが -855638017 だろうと見当がつきます．
全体的なコメント
（１）「黒」を含むタイルの枚数を数えるのと，
「白」を含むタイルの枚数を数えるのでは，結果が違
ってきます．今回は，多分，「黒」ががん組織だと思うので，
「黒」を数えるのがいいと思います
が…
熱物理工学講義資料
3
（２）私が用意したプログラムは次のようなものです．プログラミングの教科書では，
「goto を使うと
プログラムの流れがわかりにくくなるので使わないように」と書いてあることが多いですが，２
重ループから抜け出すときのように，限定的に使うのは便利だと思います．flag を用意したり，
ループ変数を再設定したりするよりも，むしろわかりやすい
かと思います．
for (x = 0; x < Width/resolution; x++) {
for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく
for (xi = x*resolution; xi< (x+1)*resolution; xi++) {
for (yi = y*resolution; yi< (y+1)*resolution; yi++) {
if (Image[yi][xi]==0) {
sum++;
goto NextTile;
}
}}
NextTile: ;
}}
今回の課題に限っては，白と黒のどちらをカウントしてもほとん
ど同じになるようです．
熱物理工学講義資料
4
第５回レポート課題
出題日：6 月 27 日（月）
提出期限：7 月 2 日（土）
提出方法：電子メールで，[email protected]
へ．
課題粘性支配のブラウン運動の解析練習として，右図のような１次
元対称２重井戸 double well 中の粒子の運動を考える．このポテン
シャル
𝑈(𝑥) = 25𝑥 4 − 100𝑥 2
は，深さが 100 程度の２つの谷をもつように適当に作ったものである．
このポテンシャル中での，粘性支配のブラウン運動
𝑑𝑥
𝑑𝑈
=−
+ 𝐹𝑟
𝑑𝑡
𝑑𝑥
を考える．ただし𝐹𝑟 はこれまでと同じく，white noise の性質をもつ random force,
{
〈𝐹𝑟 (𝑡)〉
〈𝐹𝑟 (𝑡1 )𝐹𝑟 (𝑡2 )〉
=
0
= 𝐶𝛿(𝑡1 − 𝑡2 )
である．なお摩擦係数 γ は１とする．また，第２種搖動散逸定理（１次元）より，𝐶 = 2𝑇 の関係がある．
*講義ノート p. 38, (62)式: 𝐶 = 2𝑛𝑑 𝛾𝑘𝐵 𝑇,において，次元数𝑛𝑑 = 1，摩擦係数γ = 1，Boltzmann 定数
𝑘𝐵 = 1 （つまり温度を𝑘𝐵 単位で測る）としたものである．
講義中に紹介した粘性支配のブラウン運動のプログラム（講義ノート p.82 図 5-40，あるいは以下のプログ
ラムリスト）を改造して，いくつかの温度において計算を実行し，次の２つを調べよ．ただし、初期条件は
𝑥 = −2.0 とせよ．
1.
十分に時間が経過した後の粒子の存在確率分布 𝑝(𝑥) を求め，予想される Boltzmann 分布
𝑝𝑒𝑞 (𝑥) ∝ exp [−
2.
𝑈(𝑥)
𝑇
] と比較せよ．(𝑘𝐵 = 1としたことに注意)
井戸の間を遷移する頻度を求め，その温度依存性を調べよ．
なお，横軸を反応座標 reaction coordinate と考え，２つのポテンシャル極小点をそれぞれ反応物
reactant，生成物 product に対応すると考えると，これは化学反応の反応速度を調べる簡単なモデルになっ
ている．
ヒント１「井戸を移る」というのをどのようにカウントするかが key point です．単純に，座標 x の符号が変わ
った回数（－ ⇔ ＋），と考えることもできますし，頂上付近の揺らぎを除くために，もう一工夫するのもいいで
しょう．
ヒント２温度依存性を考えるキーワード：アレニウスプロット (Arrhenius plot)
本日のレポート課題は以上です．次回は７月４日（月）
熱物理工学講義資料
5
講義ノート p.82 のプログラムを，この課題の１次元対称井戸型に書き
換えたものと２つの温度での計算例．分布を求めたり，遷移頻度を求
めたりできるように改造して欲しい．force が大きい領域でも計算でき
るように，時間刻み delta_t をやや小さくした．
// ２重井戸中の粘性支配ブラウン運動
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define
#define
#define
#define
gamma
1.0
delta_t
0.0002
TOTAL_STEP 200000
SAVE_STEP
10
// 摩擦係数
// 時間刻み
// 計算ステップ数，もっと増やしてもいい
// ファイル出力頻度
//--------------------------------------------------------double fr(double factor)
{
return ( rand( )/(double)RAND_MAX*2 - 1.0 ) * factor;
}
//--------------------------------------------------------double potential(double x) // 計算の中では使わない
{
return (25.0*x*x*x*x - 100.0*x*x);
}
//--------------------------------------------------------double force(double x)
{
return (-100.0*x*x*x + 200.0*x);
}
//--------------------------------------------------------int main( )
{
double x;
double factor, temperature;
int step;
FILE *fout;
fout = fopen("brown-viscous.dat","w");
temperature= 30.0;
// いろいろな「温度」で計算せよ
factor=sqrt((3*2*temperature)/(2*delta_t));
x=-2.0;
// 初期位置は -2.0 に固定せよ
for (step=0; step<=TOTAL_STEP; step++) {
if (step%SAVE_STEP == 0) {
printf("%9d %10.3f¥n",step,x);
fprintf(fout,"%9d %10.3f¥n",step,x);
}
x+= delta_t * (force(x)+fr(factor))/gamma;
}
fclose(fout);
return 0;
}
熱物理工学講義資料
6

Download Report