確率論

３，確率論

物理的なことをシミュレートするのには、今までのように常微分方程式や偏微分方程式を
使って決定論的に解く方法がありました。しかし、物理現象が決定論的には決まらず
確率的になってしまう場合は、シミュレーションも確率論的に解きます。

モンテカルロ法

確率論的なシミュレーションをするにはモンテカルロ法がもっとも効果的です。
モンテカルロ法では物理学的な量を、乱数を用いて計算します。
計算方法は次のようにします。

●１，初期状態の作成

まず、初期状態を作成します。

●２，試行

乱数を用いてランダムに１つの粒子を選び、これに物理的な変化を加える。
(”変化”とは、位置を移動させる，スピンを変えるなど)

●３，判定

「２」の試行で変化したポテンシャルエネルギーをΔＵとすると
次のような判定をします。

• ΔＵ＜０　　　判定は真
• 乱数Ｒを作成し、次の判定します。(０＜Ｒ＜１)
  • Ｒ＜exp{−βΔＵ}　　判定は真
  • Ｒ＞exp{−βΔＵ}　　判定は偽

※βは(ＫB Ｔ)^-1です。ｋBはボルツマン定数です。

この判定で「真」になったものは「２」の試行を採用します。
(「２」の試行が正しいと判定された)
しかし、「偽」と判定されたものは「２」の試行の結果を採用せず元のままにします。

●４，繰り返し

以上の２と３を繰り返します。
この試行を多く行うことによって実際の物理現象に近づきます。

レポート

物質内には原子があり、スピンを持っています。
このスピンが同じ方向を向くと強磁性体になります。
しかし、常温ではこのスピンの向きはランダムで常磁性体になります。
(これは、温度が絶対０度ではなく、エントロピーが増加しようとするためです。)
これが絶対０度になると、スピンによるポテンシャルエネルギーを小さくするために
スピンの向きをそろえようとします。今回のレポートは
温度Ｔを変化させて、そのときのスピンの状態をシミュレーションします。

●課題

ｍ＝�狽妬j

　　ｓij＝１，−１
　　　　　これは、電子のスピンです。
　　　　　今回は簡単のため−１と１しかないものとします。

このｍ2と温度Ｔとの関係をグラフにせよ

●解答

このシミュレーションプログラムの基本構造は次のようになります。

import java.applet.*;
import java.awt.*;

public class j30c1 extends Applet implements Runnable
{

	boolean move = false;
	double sumis,deltau,ktoj=0.01;
	int [][]is = new int[130][130];
	private Thread	 m_j30c1 = null;

	void setis() 
	{
		int i,j;
		for (i=1; i<=128; i++){
			for (j=1; j<=128; j++){
				if (Math.random()>0.5){
					is[i][j]=1;
				}else{is[i][j]=-1;}
		}}
		for (i=1; i<=128; i++){
			is[0][i]=is[128][i];
			is[129][i]=is[1][i];
			is[i][0]=is[i][128];
			is[i][129]=is[i][1];
		}
		is[0][0]=is[128][128];		is[0][129]=is[128][1];
		is[129][0]=is[1][128];		is[129][129]=is[1][1];
	}


	public void init()
	{
    	resize(500, 400);
	setBackground(Color.white);
	setis();
	}


	public void paint(Graphics g)
	{
	    int i,j,redcount=0;
	    String str;
	    g.setColor(Color.green);g.fillRect(3,3,384,384);
	    g.setColor(Color.red);
	    for (i=1; i<=128; i++){
	        for(j=1; j<=128; j++)
	            if (is[i][j]==1){g.fillRect(i*3,j*3,3,3);redcount++;}
	    }
	    g.setColor(Color.black);
	    g.drawString("ｍｍ："+Math.pow((2*redcount-128*128),2),400,350);
	}

	public void start()
	{
		if (m_j30c1 == null)
		{
			m_j30c1 = new Thread(this);
			m_j30c1.start();
		}
	}
	
	public void stop()
	{
		if (m_j30c1 != null)
		{
			m_j30c1.stop();
			m_j30c1 = null;
		}
	}

	public void run()
	{
	    int indx,indy,i;
	    while (true)
	    {
	        try{Thread.sleep(1000);}
	        catch (InterruptedException e){stop();}
	        if (move){
	            for (i=0; i<128*128; i++){
	                indx=(int)(Math.random()*128)+1;
	                indy=(int)(Math.random()*128)+1;
	                sumis=(double)(  is[indx+1][indy]+is[indx-1][indy]
	                            +is[indx][indy+1]+is[indx][indy-1]);
	                deltau=(2.0/ktoj)*(double)is[indx][indy]*sumis;
	                if (deltau<0){is[indx][indy]*=-1;}
	                else if (Math.random()<=Math.exp(-deltau))
	                {is[indx][indy]*=-1;}

	                if (indy==1){is[129][indy]=is[1][indy];}
	                if (indy==128){is[0][indy]=is[128][indy];}
	                if (indx==1){is[indx][129]=is[indx][1];}
	                if (indx==128){is[indx][0]=is[indx][128];}
	            }
	            repaint();
	        }
	    }
	}

	public boolean mouseDown(Event evt, int x, int y)
	{

		move=!move;
		if (evt.clickCount==2){
			move=true;
			setis();
		}
		if (move){m_j30c1.resume();}else{m_j30c1.suspend();}
		return true;
	}


}

●実行結果

以上のプログラムを元にして制作したのが下のものです。

左側の画面をClickするごとに動作／停止します。
doubleClickで初期状態に戻ります。
さらに「kToJ」の「０．５倍」「１」「２倍」をclickすることで
現在の数値を変更したり、１にプリセットできます。

上へ

Topに戻る