機械理工学・マイクロエンジニアリング・航空宇宙工学専攻 熱物理工学 2015（松本担当分）講義資料 松本充弘 第５回 2015/05/11 [email protected] *レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで． 前回のレポート課題について 課題 雲の形状のフラクタル解析． ＊対象図形が同じで，考え方も同じ（はず）なので，皆さん同じ結果が得られ ると思っていたのですが，「図形境界をタイルで覆う」にも，さまざまな考 え方（明らかに誤っているものも…）があったようです． ○レポートから（小西保彰さん） 雲のデジタル画像を２値化したもの Image[y*resolution + yi][x*resolution + xi] とすることでタイルごとの検 索領域にかぶりがないようにしている． データ点列に対して最小二乗法フィッティングしたグラフを青線で 示した．その結果，ハウスドルフ次元 D の値は 1.53 となった．よっ て，「雲の周囲の長さ」のフラクタル次元は 1.53 と見積もることがで きた． 松本のコメント：多分，最も素直な方法の１つ． ○レポートから（池田垣平さん） (sum0,sum1 を変数に追加) for (x = 0; x < Width / resolution; x++) { for (y = 0; y < Height / resolution; y++) { // タイルを１枚ずつチェックしていく for (xi = 0; xi < resolution; xi++){ for (yi = 0; yi < resolution; yi++){ sum0 = sum0 + (Image[(x*resolution) + xi][(y*resolution) + yi] *Image[(x*resolution) + xi][(y*resolution) + yi]); //二乗和：sum0 が 0 なら all 黒、sum0 が!=0 なら黒白 or all 白 sum1 = sum1 * Image[(x*resolution) + xi][(y*resolution) + yi]; //掛け合わせ:sum1 が 0 なら all 黒 or 黒白、sum1 が!=0 なら all 白 熱物理工学 講義資料 1 } } if (sum0 !=0 && sum1 == 0){ //この 2 つの条件が検査タイルの中に黒白タイルが含まれる条件 sum += 1; //(黒白タイル)カウントするタイル一枚追加 } else sum += 0; //( All 黒 or All 白 )カウントするタイル追加ならず } } 松本のコメント：ピクセル値の積と和を利用して，白黒両方が含ま れていることを検出しようというのは面白い方法で，原理的に は最初の人のグラフと同じものが得られるはずです．ただ残念 ながら，その意図と実際のコードが違っている気がします．あ るいは，隣のタイルに移った時に sum0, sum1 の初期化がされ ていないのが悪いのかも． ○レポートから（岸本章裕さん） 空欄部に for (xi = x*resolution ; xi <= (x+1)*resolution-1; xi++) { for (yi = y*resolution ; yi <= (y+1)*resolution-1; yi++) { if (Image[yi][xi]!=0) {sum++;break;} } if (Image[yi][xi]!=0) {break;} } というコードを追加して計算を行いました．計算結果とその他の線と の比較により，フラクタル次元はおよそ 1.9 だと推測できます. 松本のコメント：ひとつでもゼロ（黒）でないピクセルがあればカ ウントするという意図はよくわかります．でも，それだと境界 ではなく図形全体を覆っていることになり，大きな図形だと （境界の複雑さとは関係なく）ほぼ２次元ということになってしまうでしょう． ○レポートから（松下睦生さん） for (x = 0; x < Width/resolution; x++) { for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく count=0; for (xi = 0;xi<resolution;xi++) { for (yi = 0;yi<resolution;yi++) { if (Image[yi+y*resolution][xi+x*resolution]==0) count++; } } if ((count==0)||(count==resolution*resolution)) continue; sum++; } } 図より、resolution が小さいときは D～1 になり、大きいときは D ～2 になっていることがわかる。すなわち画像の分割が細かいとき にはフラクタル次元は 1 次元になり、分割が粗いときには２次元に なると考えられる。この理由について考える。分割が細かいときほ ど sum は境界線の長さに近くなる。すなわち講義ノートの海岸線 熱物理工学 講義資料 2 のように１次元的な発想で解析していることになる。一方、分割が粗いほど２次元で解析している効果が大 きくなるため D=2 に近くなると考えられる。したがって、フラクタル図形でなければ resolution が大きく なるにつれて D が大きくなるのは妥当なことであると考えられる。さらに、この図をみると D=1 からなだらか に D=2 に移行しているというよりかはある値を境に急速に D の値が変化し 1 から 2 に変わっているととるこ とができる。もし講義ノートにある海岸線のように非常に入り組んだ形をしているもののフラクタル次元を考 えられるとするとコッホ曲線のように１と２の値をとるはずである。しかし今回の雲のフラクタル次元は分割が 細かいと D=1 で直線に等しく、荒いと D=2 で平面に等しい。したがって今回の例では全体としてみると次 元が変わっていてグラフが直線的でないのでフラクタル的図形とは言えないと考えられる。 松本のコメント： 「全部が黒(count==0) でもなく 全部が白(count==resokution*resokution)でもな い」場合に continue を使ってループの初めに戻す，というのはなかなか高度ですが，完全に正 しい書き方です．初心者には if ((count!=0)&&(count!=resolution*resolution)) sum++; のほうが読みやすいとは思いますが． ○レポートから（陳愷唯さん） for (x = 0; x < Width/resolution; x++) { for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく if ( Image[x][y] != 0) //count the white area { sum += 1; } } } 松本のコメント：これだと，図形の左下のタイル１枚だけを調べていること になりますね．当然，「２次元」が得られるでしょう． ○レポートから（塩見勇祐さん） for (x = 0; x < Width / resolution; x++) { for (y = 0; y < Height / resolution; y++) { // タイルを１枚ずつチェックしていく int flag; flag = 0; if ((Image[y][x] != 0) && (Image[y + 1][x] == 0)) flag = 1; else if ((Image[y][x] != 0) && (Image[y - 1][x] == 0)) flag = 1; else if ((Image[y][x] != 0) && (Image[y][x + 1] == 0)) flag = 1; else if ((Image[y][x] != 0) && (Image[y][x - 1] == 0)) flag = 1; else flag = 0; if (flag == 1) sum++; } } 松本のコメント：これも，一枚のタイルの中だけを調べていますね． ○レポートから（篠塚尚明さん） タイル数が 10 程度のときはフラクタル次元は 0.86 程度であるが，タイ ル数の増加につれてフラクタル次元は増加し，1.05 に漸近することが分かる． 松本のコメント：多分，N-R のグラフ（左）は意図通りの結果だと思いますが，そこからどのように「フラクタル次 元」を評価したのでしょうか？ 左の両対数 グラフの傾きはおよそ -2 ぐらいと読み取れ るのですが… 熱物理工学 講義資料 3 ○レポートから（内藤一馬さん） 一部を抜粋 配布プログラムでは画像の端が解析できな い（画像サイズが 256 の倍数ではないので、 それ以上の R では端が余る）ため、雲の形全 ノイズ除外前 ノイズ除外後 体を評価できるよう改編した。（中略） 元画像をよく見ると、人間の目で見て「境界」と思われる領域以外にも、 黒い領域に白い点があったり、白い領域に黒い点が散在したりする のがわかる。今回のプログラムでは、そのような散在するノイズのよう な情報も雲の境界をなすものとしてカウントしている。（中略）よって雲 の境界線は Resolution が大きい条件ではノイズのせいで面積のよう な情報を持つってしまい、次元が不当に大きくなっていると考えられ る。そこで、より正確に解析を行うため、このようなノイズを除く方法を 考えた。穴埋め部の条件に、「白が 7.5%以下、または 92.5%以上なら ばその領域は境界線とみなさない」という条件を追加した（この数字は 試行錯誤で与えた）。このようにすると、ノイズと思われる点をきれいに 除いてカウントすることができた。全域にわたってハウスドルフ次元は 1.5 となっている。 ○レポートから（山下遼人さん） タイルが小さいとき変な結果になってしまいました。どういう風に間違 ったかよくわかりませんでした。大きい時の結果を使ったらフラクタル 次元は２と求められました。 for (xi = 0; xi < Width/resolution; xi++) { for (yi = 0; yi < Height/resolution; yi++) { // タイルを１枚ずつチェックしていく count = 0; for(x = xi; x < resolution; x++){ for(y = yi; y < resolution; y++){ if(Image[y][x] == 0) count += 1; } } if( (count != 0) && (count != resolution * resolution)) sum += 1; } } こんな風にプログラム書き換えましたがどこで間違ったのかわかりません 松本のコメント：Image[y][x] の指し示しているピクセルの場所が間違っているようですね．xi, yi は「タイ ルの背番号」なので，内側のループは，たとえば for (x=xi*resolution; xi<(xi+1)*resolution; xi++) { for (y=yi*resolution; yi<(yi+1)*resolution; yi++) { とでもすべきでしょうね． 熱物理工学 講義資料 4 松本の補足： テスト系として，「円 circle」の白黒画像も作って同じ解析をしてみました．タイルサイズが小さい 領域で，１次元よりも小さな傾きを持つように見えるのが少し不思議です．さらに簡単に，上半分黒， という完全な１次元図形の解析もやってみました．一応，きちんと１次元になっているようです． ということで，完全な理解には至っていませんが，タイルサイズが小さい場合，その中を斜めに境 界が通過しているときの処理の問題かも知れないと思っています． 実際問題としては，内藤さんの指摘のように，ノイズをどう処理するかのほうが重要だろうと思い ます． 熱物理工学 講義資料 5 第５回レポート課題 出題日：5 月 11 日（月） 提出期限：5 月 23 日（土）次回の準備の都合上，できるだけ期限内に送って下さい． 提出方法：電子メールで，[email protected] へ． 課題 粘性支配のブラウン運動の解析練習として，右図のような １次元非対称２重井戸 double well 中の粒子の運動を考える． このポテンシャル 𝑈(𝑥) = 30𝑥 4 − 90𝑥 2 − 20𝑥 は，深さが 50～100 程度の２つの谷をもつように適当に作ったも のである． このポテンシャル中での，粘性支配のブラウン運動 𝑑𝑥 𝑑𝑈 =− + 𝐹𝑟 𝑑𝑡 𝑑𝑥 を考える．ただし𝐹𝑟 は white noise の性質をもつ random force, { 〈𝐹𝑟 (𝑡)〉 〈𝐹𝑟 (𝑡1 )𝐹𝑟 (𝑡2 )〉 = 0 = 𝐶𝛿(𝑡1 − 𝑡2 ) である．なお摩擦係数 γ は１とした．また，第２種搖動散逸定理より，𝐶 = 2𝑇 の関係がある． *講義ノート p. 35, (62)式: 𝐶 = 2𝑛𝑑 𝛾𝑘𝐵 𝑇,において，次元数𝑛𝑑 = 1，摩擦係数γ = 1，Boltzmann 定数 𝑘𝐵 = 1（つまり温度を𝑘𝐵 単位で測る）としたものである． 講義中に紹介した粘性支配のブラウン運動のプログラム（講義ノート p.76 図 5-39，あるいは以下のプログ ラムリスト）を改造して，幾つかの温度において計算を実行し，次の２つを調べよ．ただし、初期条件は 𝑥 = −2.0 とせよ． 1. 十 分 に 時 間 が 経 過 し た 後 の 粒 子 の 存 在 確 率 分 布 𝑝(𝑥) を 求 め ， 予 想 さ れ る Boltzmann 分 布 𝑝𝑒𝑞 (𝑥) ∝ exp [− 2. 𝑈(𝑥) 𝑇 ] と比較せよ． 左の井戸から右の井戸へ移る頻度，および右から左への頻度の温度依存性． なお，横軸を 反応座標 reaction coordinate と考え，２つのポテンシャル極小点をそれぞれ 反応物 reactant，生成物 product に対応すると考えると，これは化学反応の反応速度を調べる簡単なモデルになっ ている． ヒント１ 「井戸を移る」というのをどのようにカウントするかが key point です．単純に，座標 x の符号が変わ った回数（－⇔＋），と考えることもできますし，頂上付近の揺らぎを除くために，もう一工夫するのもいいでし ょう． ヒント２ 温度依存性を考えるキーワード：アレニウスプロット (Arrhenius plot) 本日のレポート課題は以上です． 来週は外国出張のため休講にします（当初の予定通り）．次回は 5 月 25 日． 熱物理工学 講義資料 6 講義ノート p.76 のプログラムを，この課題の１次元非対称井戸型 に書き換えたものと２つの温度での計算例．分布を求めたり，遷 移頻度を求めたりできるように改造して欲しい．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 (30.0*x*x*x*x - 90.0*x*x - 20.0*x); } //--------------------------------------------------------double force(double x) { return (-120.0*x*x*x + 180.0*x + 20.0); } //--------------------------------------------------------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; } 熱物理工学 講義資料 7
© Copyright 2024 ExpyDoc
