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

機械理工学・マイクロエンジニアリング・航空宇宙工学専攻
熱物理工学 2016（松本担当分）講義資料
松本充弘
最終回 2016/07/11
[email protected]
*レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで．
前回のレポート課題について
課題 Monte Carlo 法の応用として，Travelling Salesman Problem (TSP) を
焼きなまし（simulated annealing）法で扱う．
ひと言でいえば，どれだけ経路長 L が短くなったかが勝負．サンプルプログラ
ムの最急降下法（SD 法）では L=1498.4 だった．
○レポートから（大谷拓也さん）
例：SD 法の結果
L=1034.34519
○レポートから（高橋弘明さん）
温度を一定の割合で減らす場合，局所解に陥っているのか見た目
ではまだ短くなりそうなものでも，解が収束してしまった．
そこで，ある程度解を収束させた時点で，一度温度を上昇させる
ことにより解に幅を持たせ，改めて収束させることでより優れた
解が得られた．
L=1033.3
○レポートから（三上慎司さん）
ちなみに，最近の強い将棋ソフトは全幅探索だけではなく，自然
な指し手を重視する選択探索も用いたハイブリッドである。たと
えば上図を見ても分かるが，「サイト密度が小さい部分に経路を
引かない」「経路は交差しないように引く」という評価関数を作
れば大幅に改善できるかと思われる。これを巡回セールスマン問
題にも適用できれば，さらに発展するのではないだろうかと課題
に取り組みながら思った。
L=964.8
熱物理工学講義資料
1
○レポートから（大仲浩徳さん）
L=989
○レポートから（上原水優さん）
L=1057.2
○レポートから（樋口智哉さん）
より短い経路を求めるため、main 関数の中に for ループを組み込
み 20000step 時の距離が今までの最短距離よりも短ければ、
座標、
距離を更新するというプログラムを作成した。for ループは 1000
回回した。焼きなましについては、初期温度 T=1600、1000step
ごとに T が半減するという条件を採用し、乱数は
srand((unsigned) time(NULL));によって、その都度初期化した。
L=966.72706
○レポートから（西岡寿朗さん）
∆𝐿
松本のコメント：MC 法では，exp⁡(− )の確率で新しい経路を選択するわけですから，
「温度」は，経
𝑇
路長の変化と同程度であるように選ぶのが効率がよさそうですね．今回は，経路長そのものが
1000 のオーダーをもちますので，T も 1000 程度（あるいはそれ以上）から始めて，
「冷却」して
いくのがいいのかなと思います．なお，𝑇 → 0での挙動についてですが，trial の経路変化を乱数
で決定していますので，かならずしも最急降下法と一致するとは限りません．というよりは，最
も経路が短くなる方向に変化させているわけではないので，真の「最急降下法」ではないことに
注意してください．離散的な問題では gradient が定義できませんから．
熱物理工学講義資料
2
○レポートから（今村建太さん）
L=1042.9 (15000step) L=1044.2 (20000step)
○レポートから（原田浩行さん）
L=966.72706
○レポートから（佐藤広野さん）
L=1046.056
松本のコメント：なかなかいい経路が得られたと思うが，１つだけ注意．多分エクセルで経路を描い
たのだと思うが，こんなときはスムージング機能を使わず，線分でつなぐべき．
○レポートから（広田翔さん）
松本のコメント：出発点に戻ってこない経路ですし，
「都市」の配置も異なるようなので，別の問題を
解いたことになりますね．
熱物理工学講義資料
3
○レポートから（馬場量子さん）
L=1096.62757
松本のコメント：ほかの人とは少し違う経路が得られましたね．これも十分に短い経路です．
○レポートから（前田隆馬さん）
L= 1149.19
松本のコメント：これも，やや違うパターンの経路ですが「アリ」ですね．
○レポートから（西田裕悟さん）
20000 ステップをかけて計算した結果、距離は約 1173 となった。まだ自己交差など
が残っており、人の手を加えれば容易に距離を短縮することができるが、再急降下法
の結果よりかなり改善されていることが分かる。2 つの点の順序を交換するだけでは
自己交差を取り除くことはできず、ある程度の高温をしばらく維持することが自己交
差を取り除くのに必要であると考えられるが、この計算では温度の下げ方が急すぎた
ため、自己交差を十分に除去しきれなかったものと考えられる。
（中略）
また、2000000（2.0e+6）ステップをかけて、焼きなまし法を用い計算した結果、距
離は 968 となった。
L= 968
松本のコメント：そう，うまく「交差」を取り除くというのが TSP の key point なのですが，これは
残念ながら，sample program にあるような単純な「都市」の入れ替えではうまく実現できませ
ん．
○レポートから（小山徹さん）
最後に通る点（31 点目）を最初に通れば赤矢印で示したように経路を短絡で
きるように見える。1 点目(32 点目)と 31 点目の交換を許せばこの短絡は容易
に起こると考えられる。しかし、このプログラムでは 1 点目(32 点目)と 31 点
目の交換を許していないので、この短絡が起こるためには、31 点目と 1 点目
と入れ替え、31 点目と 2 点目と入れ替え、31 点目を 3 点目と入れ替え、
・・・
31 点目を 30 点目と入れ替えるという操作を行う必要がある。この操作の間は
distance が増大しているので、この状態からいくらステップ数を増やしても
赤矢印に示す短絡は起きそうにないと思った。
200
100
0
0
100
200
L=983.70
熱物理工学講義資料
4
最終課題について
３つ，お願いがあります．
(1) 授業アンケートに回答してください．多分，教務からも案内があったかと思いますが，KULASIS
からリンクされた KULIQS という web system から，スマートフォンやＰＣから回答ができます．
回答期限：９月３０日ですが，早めに（できれば７月中に）お願いします．
科目：10G005
熱物理工学
吉田英生／松本充弘
(2) (1)の授業アンケートとは別に，大学のサーバ上に，吉田先生が「熱物理工学」の最終アンケート
フォームを作成されました．
https://www.t.kyoto-u.ac.jp/fs/s-es/netsubutsurikogaku_enquete
こちらへの回答もお願いします．授業で取り上げた
内容，レベル，基幹科目としての妥当性，改善の提
案など，何でも結構です．もちろん成績評価とは関
係しません．
回答期限：７月３１日（日）とします．
(3) レポート課題
次の課題についてのレポートをいつものように，
電子メールで，
[email protected] へ
提出期限：同じく７月３１日（日）
なお，未提出のレポート課題も，この期限までは受け付けますのでどうぞ．
また，この課題について解説をする機会がありませんので，提出締切後の適当な時期（８月中旬
ごろまでに）に，私のホームページ
http://www.mitsuhiromatsumoto.mech.kyoto-u.ac.jp/
に解説を掲載します．興味のある方は，後ほどアクセスしてみてください．
課題 Fokker-Planck (FP) 方程式を解いて，非平衡エントロピーの時間
変化を体験してみよう．第５回のレポート課題で用いたのと同じ，対称２
重井戸型ポテンシャル
𝑈(𝑥) = 25𝑥 4 − 100𝑥 2
を使うことにする．このポテンシャル中の粘性支配 FP 方程式
𝜕𝑃 1 𝜕 𝑑𝑈
𝜕2
= [ ( 𝑃) + 𝑘𝐵 𝑇 2 𝑃]
𝜕𝑡 𝛾 𝜕𝑥 𝑑𝑥
𝜕𝑥
の数値計算プログラムは次のページの通りである．非平衡エントロピーも
求められるように改造し，いくつかの温度において，確率分布やエントロピーの時間変化を調べよ．
これで本年度の熱物理工学は終了です．ご苦労さまでした．
熱物理工学講義資料
5
1 次元 Fokker-Planck 方程式の数値計算プログラム
私のホームページ http://www.mitsuhiromatsumoto.mech.kyoto-u.ac.jp/ からダウンロード可．
注意
微分方程式を数値積分するのに素朴な陽解法を用いているため，特に低温
（convection term >> diffusion term）において，数値不安定性が起こりやすい．振
動的な分布があらわれたら，時間刻み DT を小さくしてみるとうまく計算できるかも
しれない．
// Numerical Integration of Fokker-Planck Equation: 1-D Diffusion in Potential Field
#include <stdio.h>
#include <math.h>
#define MAX_STEP 100000
#define SAVE_STEP 1000
double
double
double
DT =
0.000002;
// 時間刻み．計算精度が悪いときは小さくしてみよう
GAMMA = 1.0;
TEMPERATURE = 20.0; // 温度：いろいろと変えてみよう
#define NX_CELL 200
double XMIN = -4.0;
double XMAX = 4.0;
double
prob[NX_CELL];
double
dprob[NX_CELL];
double prob_eq[NX_CELL];
double dx;
// 領域の分割数
// 領域の下限
// 領域の上限
// ヒストグラムの幅
//------------------------------------------------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;
}
計算例：T=20.0, 50,000step 目の分布関数
//------------------------------------------------void initial( )
{
int ix;
double sum, factor, x;
dx=(XMAX-XMIN)/NX_CELL;
for (ix=0;ix<NX_CELL;ix++) {
prob[ix]=0.0;
}
ix=(int)((-2.0-XMIN)/dx);
prob[ix ]=0.2/dx;
// 初期分布を x=-2 近くでステップ状に与えた
prob[ix+1]=0.2/dx;
prob[ix+2]=0.2/dx;
prob[ix+3]=0.2/dx;
prob[ix+4]=0.2/dx;
sum=0.0;
for (ix=0;ix<NX_CELL;ix++) {
// 平衡分布
x=dx*ix+XMIN;
prob_eq[ix]=exp(-potential(x)/TEMPERATURE);
sum+=prob_eq[ix];
}
factor=1/(dx*sum);
// 規格化しておく
for (ix=0;ix<NX_CELL;ix++) {
熱物理工学講義資料
6
prob_eq[ix]*=factor;
}
}
//------------------------------------------------void output(int nstep)
{
int ix, iv;
char cdum[100];
FILE *fout;
printf("%6d¥n",nstep);
sprintf(cdum,"fp%6.6d.dat",nstep);
fout=fopen(cdum,"w");
for (ix=0;ix<NX_CELL;ix++) {
// 現時点の分布関数と平衡分布をファイルに出力
fprintf(fout,"%6.2f %9.5f %9.5f¥n",ix*dx+XMIN,prob[ix],prob_eq[ix]);
}
fclose(fout);
}
//-------------------------------------------------------int main( )
{
int
ix, nstep;
double dp;
double x,xm,xp,sum;
double f, d;
initial( );
nstep=0;
output(nstep);
for (nstep=1;nstep<=MAX_STEP;nstep++) {
for (ix=1;ix<NX_CELL-1;ix++) {
// 増分を求める
x = dx*ix+XMIN;
xm = x-dx;
xp = x+dx;
dprob[ix] = -(force(xp)*prob[ix+1]-force(xm)*prob[ix-1])/(2*dx)/GAMMA;
dprob[ix]+= (prob[ix+1]-2*prob[ix]+prob[ix-1])/(dx*dx)*(TEMPERATURE/GAMMA);
}
sum=0.0;
for (ix=1;ix<NX_CELL-1;ix++) {
prob[ix] += dprob[ix]*DT;
if (prob[ix] <= 0.0) {
prob[ix] = 0.0;}
else {
sum+=prob[ix];
}
}
// convection term
// diffusion term
// 確率密度の更新と再規格化，ただし両端はゼロのまま
// 確率密度は非負
sum*=dx;
for (ix=1;ix<NX_CELL-1;ix++) {
// 再規格化
prob[ix] /= sum;
}
if (nstep%SAVE_STEP==0) output(nstep); // 分布をファイルに出力
}
return 0;
}
熱物理工学講義資料
7

Download Report