計算物理学　第１章序章
２－２．DLAによる樹木状クラスタのフラクタル次元

本稿は、凝縮系物理を中心として広範囲にわたる計算手法に関して、そのアルゴリズムの導出から実装までを丁寧にまとめられていることで知られている計算物理学 J.M.ティッセン (著), (訳) 松田 和典, 道廣 嘉隆, 谷村 吉隆, 高須 昌子, 吉江 友照を元に、 本書で取り扱われているテーマについて、具体的にC++言語にてプログラミングし、計算結果を gnuplot や OpenGL を用いて図や動画にすることを目的とします。

前節でシミュレーションしたDLAによる樹木状クラスタの形状を定量的に取り扱うことを考えます。 吸着率が p=1 の場合と p=0.1 の場合では、粒子数が同じにもかからわす、クラスタの広がり方に違いがあることがわかります。 これは粒子数とクラスタサイズとの関係が得られれば、クラスタの形状の特徴の一端を知ることができることを意味しています。 クラスタサイズを次の式で定義します。

(1)

ただし、

(2)

は、クラスタの重心を表します。 式(1)はクラスタの旋回半径と呼ばれる量で、 クラスタの重心からの位置の平均を表しています。 式(1)は式(2)を組み合わせることにより、

(3)

と変形することができます。 DLAシミュレーションの結果と式(3)から、 旋回半径 R と粒子数 N との関係を次式で表します。

(4)

Duffing振動子の時とは異なるフラクタル次元となります。

計算結果

縦軸：粒子数 N
横軸：旋回半径 R_g
傾きからフラクタル次元 Df を見積もります。 もし Df = 2 ならば、２次元平面を隙間なく埋め尽くしていることになります。

粒子数 N = 10000
吸着率 p = 1

粒子数 N = 10000
吸着率 p = 0.1

p=1.0 の場合のフラクタル次元は、D_f ≒ 1.71 程度になることが知られています。 p=0.1 の場合は、p=1.0 に比べて大きな値を取ることを確認できました。

C++プログラム

前節で計算した樹木状クラスタのデータを読み込んで、旋回半径と粒子数との関係を出力するプログラムです。

 
//////////////////////////////////////////////////////////////////////////
// 旋回半径と粒子数の関係
//////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdlib.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <cstdio>
#include <iomanip>
#include <stdio.h>
#include <direct.h>
using namespace std;
double PI = acos(-1.0);
double e = 2.7182818284590452354;
string folder = "data", ff="";  //保存フォルダ名

const int dimension = 2;//次元
const int N_max = 20000;//クラスタ構成粒子数
const int x_0 = 10000;//中心座標
const int y_0 = 10000;//中心座標
int core_number = 0; // 固まった数

char str[200];
string str1;
int main(){
////////////////////////////////////////////////////////
//ファイルの読み込み ///////////////////////////////////
////////////////////////////////////////////////////////
  for(int n=10; n<=10; n+=1){ //<--------------------------ループの指定(p=0.1～1.0)
    ifstream fin_position;
    sprintf_s(str, "/DLA-%d.data", n); str1 = folder + str; 
    fin_position.open(str1.c_str());

    int r[N_max+1][2]; //粒子の座標
    double r0[2] ={0.0}; //重心位置
    core_number = 1;
    while(!fin_position.eof()){
      int ii, jj;
      fin_position >> ii; 
      fin_position >> jj; 
      if(ii<-x_0 || ii>x_0 || jj<-y_0 || jj>y_0) break;
      r[core_number][0] = ii;//x座標
      r[core_number][1] = jj;//y座標
      r0[0] += r[core_number][0];
      r0[1] += r[core_number][1];
      core_number++;
    }
    fin_position.close();

////////////////////////////////////////////////////////
    //重心位置の決定
    r0[0] = r0[0]/double(core_number); 
    r0[1] = r0[1]/double(core_number);
    ofstream fout_R;
    sprintf_s(str, "/RN-%d.data", n); str1 = folder + str; 
    fout_R.open(str1.c_str());
    double R_sum=0; 
    double R; //旋回半径
    double R_=1.0;
    int i_=1;
    for(int i=1; i<=core_number; i++){
      R_sum += pow(double(r[i][0]) ,2) + pow(double(r[i][1]) ,2);
      if(i%200==0) {
        R = sqrt(R_sum/double(i) - (pow(double(r0[0]) ,2) + pow(double(r0[1]) ,2)));
        fout_R << R << " " << i << " " <<log(double(i)/double(i_))/log(R/R_)  << endl;
        i_ = i;
        R_ = R;
      }
    }
    fout_R.close();
    //cin.get();
  }
}

課題メモ

１．フラクタル次元を直接計算するプログラムを作成する。
２．３次元DLAによる樹木状クラスタのフラクタル次元の計算を行う。

ミスの指摘、質問、コメントなどがございましたら、HTMLフォーム、もしくはinfo:natural-science.or.jpから、お待ちしております。

目次

第１章　序章

• １－１．Duffing振子のカオス的運動のシミュレーション
• １－２．Duffing振子の位相空間上でのカオス的運動のアニメーション
• １－３．Duffing振子の位相空間上のストレンジ・アトラクタに対するフラクタル次元
• ２－１．DLAによる樹木状クラスタのシミュレーション
• ２－２．DLAによる樹木状クラスタのフラクタル次元

第２章　対称ポテンシャルによる量子散乱

ゴール：水素原子のクリプトン原子による散乱に対する散乱断面積の計算

• ・ヌメロフ法による常微分方程式の解法
• ・球ベッセル関数
• ・ルジャンドル多項式
• ・位相シフト→断面積を計算するアルゴリズム

その他

• ・擬似乱数生成器：メルセンヌ・ツイスタの使い方
• ・ランダムウォークによる擬似乱数の比較：C言語rand関数とメルセンヌ・ツイスタ(mt19937ar)