本文へスキップ
モンテカルロシミュレーション
11. DLAシミュレーション
 DLA(Diffusion-limited aggregation)とは日本語では拡散律速凝集といい、ブラウン運動する粒子が、核となる粒子に結合することで粒子の集合体(クラスタ)が成長する過程のことです。その集合体を凝集体といいます。雪の結晶や雷の成長などに見られ、美しいフラクタル構造となっていることが知られています。ブラウン運動は6章で詳しく勉強したので、その応用としてDLAを描いてみましょう。

なお、この章のプログラムは、16ビットのTurboCでは、配列のメモリが不足するため、RasoberryPIのUNIXで実行します。

 はじめの例は元になる核は地表とし、空中の1点を発した雪が降りそそぎ、地表に到達すると、下から結合され、徐々に積み上がり、樹氷のようになる姿を描きます。雪は直線的に下に落ちるのではなく、右に左にゆらゆらと落下します。その様子はランダムウォークそのものです。

プログラムのポイント
シミュレーションプログラムをs1110.cに示します。
行15のように、画面に相当する状態配列を用意します。値は1あるいは0で、その位置にすでに雪があるかないかを示すこととします。初期状態は行26のように、地表のみ、1としておきます。
雪は、行31のように、画面中央トップから降りそそぐとしました。行39から40が雪はゆらゆら降るさまを示したもので、ランダムウォークです。つぎのタイミングで、粒子の位置を、x , yとし、右に行く割合が40%、左に行く割合が40%、残り20%が真下に落ちるとしています。雪が下に位置を変えるたびに、行41でその位置の状態がすでに雪で占められているかどうか調べます。雪がすでにあれば、行46のように一つ手前 x0, y0 にもどり、その位置の状態を、1に変えます。これが凝集体に捕捉された処理です。
適当な高さに成長したら、処理を終えます。

アウトプット
図11.1.1に示します。

  s1110.c   樹氷
 /* S1110.C
 *        DLA 樹氷
 *         (C) H.Ishikawa  2018/11/10
 */

#include "compati.h"
#include "window.h"

#define X_MAX    640
#define Y_MAX    480

int main(void){
	int x0,y0;    //雪の現在地     
	int x,y;     //次の位置              
	int s[X_MAX + 10][Y_MAX + 10] = {};  //状態配列 1:雪がある 0:雪がない
	float random;//乱数
	/*グラフィクの準備*/
	opengraph();
	SetWindow(0,  0, 0, X_MAX, Y_MAX, 0,HIGHT, WIDTH,0);  
	srand(17);
	x = 0.0;
	y = 0.0;
	x0 = 0.0;
	y0 = 0.0;

	for (x = 0; x <= X_MAX ; x ++) {
		s[x][0] = 1;      //初期状態 地上を状態1に
	}
	while(y0 < Y_MAX - 20){//高さがY_MAX - 20まで成長したらとめる
		//雪の起点
		x = X_MAX / 2;
		y = Y_MAX;
		while (y > 0){
			random = RAND();
			if(random > 0.6){//雪は右へ
				if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
			}else if (random < 0.4){雪は左へ
				if (x > 0) { x= x - 1;}else{ x = 0;}
			}
			y = y - 1;//雪は下に
			if (s[x][y] == 0){//この位置に雪がないとき
				PutPixel(0, (double)x, (double)y ,1);//雪の軌跡を青で描く
				x0 = x;
				y0 = y;
			}else{//この位置にすでに雪があれば一つ前の位置を状態1に
				s[x0][y0] = 1;//捕捉
				PutPixel(0, (double)x0, (double)y0 ,15);//白をプロット
				break;//次の雪を降らせる
			}
		}
  	}	
	getch();
	closegraph();
	return 0;
}


図11.1.1  s1110.c   樹氷


 次の例はDLAとして有名な凝集体です。金属イオンが溶けた溶液に電極を入れると、金属が析出し、徐々に植物のような形に育っていきます。これを金属樹、金属葉とよばれます。これもフラクタル構造になることが知られています。

プログラムのポイント
 下の図11.2.1に従って凝集体のシミュレーションを説明します。まず核を一つ画面中央におきます。粒子は中央から等距離の円周上から、ランダムに発生させます。そこから2次元のランダムウォークをさせます。その位置x ,yの状態を調べ空いていれば(0であれば)ランダムウォークを続けます。もしすでに粒子で占められていれば(1であれば)、ひとつ戻って、隣接する位置x0 ,y0の状態を1に遷移させます。この粒子はここで捕捉されましたので、次の粒子を発生させ処理を繰り返します。
 プログラムS1120.Cは、凝集体のシミュレーションです。x0 ,y0の状態を1に変化させるときに、白をプロットします。中央の核に粒子が捕捉され、凝集がはじまります。
図11.2.1 凝集体のシミュレーションの方法 文献(16)

  s1120.c  凝集体のシミュレーション   
/* s1120.c
 *        凝集体のシミュレーション
 *         (C) H.Ishikawa  2018/11/20
 */

#include <math.h> //  sin()  cos()  M_PI
#include "compati.h"
#include "window.h"

#define X_MAX    480
#define Y_MAX    480


#include <time.h>

int main(void){
	int x0,y0;    //粒子現在位置   
	int x,y;    //粒子次の位置             
	int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
	long t = 0;
	int color;
	double th;
	double r;

	/*グラフィクの準備*/
	opengraph();
	SetWindow(0,  0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2, HIGHT, (WIDTH + HIGHT) / 2,0);  

	/*メイン*/
	srand(107);
	x = 0;
	y = 0;
	x0 = 0;
	y0 = 0;
	
	s[X_MAX / 2][Y_MAX / 2] = 1;//中央に核をおく
	
	while(t <= 15000){
		/*粒子発生*/
		th = 2.0 * M_PI *RAND(); 
		x = (int)(X_MAX /2.0 * cos(th) + X_MAX /2.0) ;		
		y = (int)(Y_MAX /2.0 * sin(th) + Y_MAX /2.0) ;
		while (1){
			/*ランダムウォーク*/
			r = RAND();
			if(r>0.67){//右へ
				if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
			}else if (r < 0.33){//左へ
				if (x > 0) { x= x - 1;}else{ x = 0;}
			}//このほかはとどまる			
			r = RAND();
			if(r>0.67){//上へ
				if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
			}else if (r < 0.33){//下へ			
				if (y > 0) { y= y - 1;}else{ y = 0;}
			}//このほかはとどまる
			
			/*判定*/
			if (s[x][y] != 1){
				PutPixel(0, (double)x, (double)y,1);
				x0 = x;
				y0 = y;
			}else{/*捕捉*/			
				s[x0][y0] = 1;
				t = t + 1;
				PutPixel(0, (double)x0, (double)y0 ,15);
				break;/*次の粒子へ*/
			}
			
		}
  	}
	getch();
	closegraph();
	return 0;
}



アウトプット
s1120.c
実行します。結果を図11.2.2に示します。処理時間約10分
図11.2.2  s1120.c   凝集体のシミュレーション



 忍石(しのぶいし)とは、石灰岩・頁岩(けつがん)などの割れ目に二酸化マンガンなどの樹枝状結晶が付着して、シダ類の化石のように見えるものです。江戸時代の石のコレクター、木内石亭の著書「雲根志」には、忍草石として登場しており珍重されました。

 この形で成長する結晶は多く、冬の窓に付く霜の雪片もこの一種です。

図11.3.1  しのぶ石の例 (Wikipedia)


プログラムのポイント

しのぶ石のシミュレーションプログラムはs1130.c のとおりです。

初期化は行26のように地表のみ、1としておきます。粒子は、行30のように画面全体の任意の位置からランダムに発生するとします。粒子のランダムウォークは2次元で、1/3の割合で右、左、とどまるとします。上下方向も同様です。凝集体に捕捉される処理はs1120.cと同様です。

  s1130.c   しのぶ石
/* S1130.C
 *        しのぶ石
 *         (C) H.Ishikawa  2018/11/10
 */

#include "compati.h"
#include "window.h"

#define X_MAX    640
#define Y_MAX    480

int main(void){
	long x0,y0;    //現在地     
	long x,y;     //次の位置              
	long s[X_MAX + 1][Y_MAX + 1] = {};  //状態配列 1:粒子がある 0:粒子がない
	float r;//乱数
	/*グラフィクの準備*/
	opengraph();
	SetWindow(0,  0, 0, X_MAX, Y_MAX, 10,HIGHT-10, WIDTH-10,10);  
	srand(1003);
	x = 0.0;
	y = 0.0;
	x0 = 0.0;
	y0 = 0.0;
	//初期化
	for (x = 0; x <= X_MAX ; x ++) {
		s[x][0] = 1;      //地上のところどころを状態1に
	}
	while(y0 < Y_MAX - 20){//高さがY_MAX - 20まで
		//粒子の起点は画面全体
		x = (int)(RAND() * X_MAX);
		y = (int)(RAND() * Y_MAX);
		while (y > 0){	
			//2次元のランダムウォーク
			r = RAND();
			if(r>0.67){//右へ
				if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
			}else if (r < 0.33){//左へ
				if (x > 0) { x= x - 1;}else{ x = 0;}
			}//このほかはとどまる			
			r = RAND();
			if(r>0.67){//上へ
				if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
			}else if (r < 0.33){//下へ			
				if (y > 0) { y= y - 1;}else{ y = 0;}
			}//このほかはとどまる
			
			if (s[x][y] == 0){//この位置に粒子がないとき
				PutPixel(0, (double)x, (double)y ,1);//粒子の軌跡
				x0 = x;
				y0 = y;
			}else{//この位置にすでに粒子があれば一つ前の位置を状態1に
				s[x0][y0] = 1;
				PutPixel(0, (double)x0, (double)y0 ,15);//白をプロット
				break;//次の粒子を発生させる
			}
		}
  	}
	getch();
	closegraph();
	return 0;
}



アウトプット
s1130.c
を実行します。結果を図11.3.2に示します。これを描くのに、約30分かかります。

図11.3.2  s1130.c   しのぶ石シミュレーション


.核をどのように初期条件としてあたえるかによって、さまざまなDLAを描くことができます。核を円周上に置き、粒子のスタートは、中央としたとき、どのようなDLAが描けるか、トライしなさい。
  (解答) 
  a1110.c   円周上に核
/* a1110.c
 *        DLA 円周上に核
 *         (C) H.Ishikawa  2018/11/20
 */

#include <math.h> // sin()  cos() 
#include "compati.h"
#include "window.h"

#define X_MAX    480
#define Y_MAX    480

int main(void)
{
	
	int x0,y0;    //粒子現在位置   
	int x,y;    //粒子次の位置             
	int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
	long t = 0;
	int color;
	double th;
	double r;

	/*グラフィクの準備*/
	opengraph();
	SetWindow(0,  0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2, HIGHT, (WIDTH + HIGHT) / 2,0);  

	/*メイン*/
	srand(13);
	x = 0;
	y = 0;
	x0 = 0;
	y0 = 0;
	//初期化 円周上に核
	for (th = 0 ; th < 2.0 * M_PI; th = th + 0.0001){
		x = (int)(X_MAX /2.0 * cos(th) + X_MAX /2.0) ;
		y = (int)(Y_MAX /2.0 * sin(th) + Y_MAX /2.0) ;
		s[x][y] = 1;
	}
	
	while(t < 25000){
		/*粒子発生は画面中央*/
		x = (int)(X_MAX /2.0) ;
		y = (int)(Y_MAX /2.0) ;
		while (1){
			/*ランダムウォーク*/
			r = RAND();
			if(r>0.67){//右へ
				if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
			}else if (r < 0.33){//左へ
				if (x > 0) { x= x - 1;}else{ x = 0;}
			}//このほかはとどまる			
			r = RAND();
			if(r>0.67){//上へ
				if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
			}else if (r < 0.33){//下へ			
				if (y > 0) { y= y - 1;}else{ y = 0;}
			}//このほかはとどまる
			
			/*判定*/
			if (s[x][y] != 1){
				PutPixel(0, (double)x, (double)y,1);
				x0 = x;
				y0 = y;
			}else{/*捕捉*/			
				s[x0][y0] = 1;
				t = t + 1;
				PutPixel(0, (double)x0, (double)y0 ,15);
				break;/*次の粒子へ*/
			}
		}
  	}
	getch();
	closegraph();
	return 0;
}

図11.a.1  a1110.c の結果  


.核を正方形の4辺 に置き、粒子のスタートは、画面全体としたとき、どのようなDLAが描けるか、トライしなさい。
  (解答) 
  a1120.c   正方形の4辺に核を置く
/* a1120.c
 *        DLA 4辺に核
 *         (C) H.Ishikawa  2018/11/27
 */

#include <math.h> // sin()  cos() 
#include "compati.h"
#include "window.h"

#define X_MAX    480
#define Y_MAX    480

int main(void)
{
	
	int x0,y0;    //粒子現在位置   
	int x,y;    //粒子次の位置             
	int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
	long t = 0;
	int color;
	double th;
	double r;

	/*グラフィクの準備*/
	opengraph();
	SetWindow(0,  0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2 +10, HIGHT-10, (WIDTH + HIGHT) / 2 -10 ,10);  

	/*メイン*/
	srand(103);
	x = 0;
	y = 0;
	x0 = 0;
	y0 = 0;
	//初期化 正方形の4辺に核を置く	
	for (x = 0; x <= X_MAX ; x ++) {
		s[x][0] = 1;      
	}
	for (x = 0; x <= X_MAX ; x ++) {
		s[x][Y_MAX] = 1;      
	}
	for (y = 0; y <= Y_MAX ; y ++) {
		s[0][y] = 1;      
	}
	for (y = 0; y <= Y_MAX ; y ++) {
		s[X_MAX][y] = 1;      
	}
	while(t <100000){
		//粒子の起点は画面全体
		x = (int)(RAND() * X_MAX);
		y = (int)(RAND() * Y_MAX);
		while (1){
			/*ランダムウォーク*/
			r = RAND();
			if(r>0.67){//右へ
				if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
			}else if (r < 0.33){//左へ
				if (x > 0) { x= x - 1;}else{ x = 0;}
			}//このほかはとどまる			
			r = RAND();
			if(r>0.67){//上へ
				if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
			}else if (r < 0.33){//下へ			
				if (y > 0) { y= y - 1;}else{ y = 0;}
			}//このほかはとどまる
			
			/*判定*/
			if (s[x][y] != 1){
				PutPixel(0, (double)x, (double)y,0);
				x0 = x;
				y0 = y;
			}else{/*捕捉*/			
				s[x0][y0] = 1;
				t = t + 1;
				//経過時間にしたがって色を変える
				PutPixel(0, (double)x0, (double)y0 ,(int)(t/6000 + 1));
				break;/*次の粒子へ*/
			}
		}
  	}
	getch();
	closegraph();
	return 0;
}


図11.a.2  a1120.c の結果