몬테칼로 시뮬레이션 (Monte Carlo Simulation)

4분의 1 원이 한 변의 길이가 1인 정사각형에 들어있는 경우를 생각하자. 이 원이 정사각형에서 차지하는 면적의 비율은
 --- (1)
이다.

만약에 정사각형 안에 무작위로 돌을 던진다면 원에 맞을 확률은
 --- (2)
이다.

(1)의 면적 비율과 (2)의 확률은 일치하여야 한다.
 --- (3)

따라서 π는 (3)을 이용하여 계산을 할 수 있다.

돌을 많이 던지면 던질 수록 우리는 정확한 π를 계산할 수 있게 된다.

이 방법은 몬테칼로 시뮬레이션의 기본 개념인 ‘무작위와 확률을 이용하여 계산을 한다‘ 를 정확히 보여주는 대표적인 예이다.

플래시로 보기

Nt는 던진 돌의 수, Nh는 원에 맞은 돌의 수.

앞선 돌을 던져 π를 구하는 방법 보다 더 원조 격인 뷔퐁의 바늘 던지기를 통하여 π를 계산할 수 있다.

바늘의 길이 보다 더 큰 간격으로 떨어진 두 평행선을 두고 바늘을 아무렇게나 던졌을 때, 바늘이 선에 걸칠 확률을 이용하여 π를 구하는 방식이다.

바늘의 길이가 L이고 두 평행선 사이의 거리는 d이다.

바늘의 중심(L/2)이 선에 걸쳐 θ 만큼의 각도를 이루고 있다면 바늘끝이 평행선과 이루는 수직 거리 s는
--- (1)
이다. 이 s안에 들어야 바늘이 선에 걸치게 된다. 따라서 θ에 대해 0부터 π까지 적분을 하면 바로 바늘이 선에 걸치는 범위의 면적이 된다. 선은 아래 위 2개가 있으므로 두배를 해주면 된다.
--- (2)

바늘이 떨어질 수 있는 범위는 s는 0부터 d까지, θ는 0부터 π까지 이다. 따라서 이 범위가 만드는 직사각형의 면적은
πd --- (3)
이다.

(2)와 (3)에서 바늘이 선에 걸칠 확률은
--- (4)
이다. 이 값은 바늘을 N번 던져 N_h번 선에 걸릴 확률과 일치한다.
--- (5)

따라서 우리가 원하는 π는

이다.

플래시로 보기

위의 플래시는 서울대학교 이구철 교수님께서 운영하시는 물리로 배우는 플래시 (http://ldap.icpr.or.kr/forum/) 에서 가지고 온 것입니다. 더 자세한 설명과 소스 파일이 "강의19. 뷔퐁의 바늘 - 몬테칼로 기법의 원조" 에 있습니다

플래시로 보기

위의 플래시는 서울대학교 이구철 교수님께서 운영하시는 물리로 배우는 플래시 (http://ldap.icpr.or.kr/forum/) 에서 가지고 온 것입니다. 더 자세한 설명과 소스 파일이 " 강의20. 뷔퐁의 바늘 - 왜 성립하나?
" 에 있습니다

평형상태가 되는 충분조건 (필요조건은 아니다)으로 detailed balance가 있다. 이 조건은 A의 상태에서 B의 상태로 전이하는 확률과 B의 상태에서 A의 상태로 전이하는 확률이 같은 것이다.

W(A→B)는 A에서 B로 가는 전이확률, P(A,t)는 시간 t에 A일 확률이다. 몬테칼로 방법이 유효하기 위해서는 A의 상태에서 B의 상태로 가는 전이확률이 detailed balance를 만족해야 한다.

Boltzmann 분포를 사용하여 detailed balace가 어떻게 적용되는지 살펴보자.

이 계산에서 알아야 할 것은 k, T, E 뿐이다.

전이확률 W를 선택하는 데는 많은 여러 방법이 있지만 Metropolis 알고리듬은 아주 간단한 방법으로 detailed balance를 만족시키고 있다.

따라서 E(B) > E(A)인 경우에

E(B) <= E(A)인 경우는

가 되어 Metropolis 알고리듬은 모든 통계 계에서 detailed balance를 만족한다.

k는 볼츠만 상수이고, T는 절대온도이다. 이 두 값은 쉽게 정의가 가능한 값들이다. 따라서 에너지 E를 계산할 수 있다면 우리는 Metropolis 알고리듬을 이용한 몬테칼로 시뮬레이션을 통하여 쉽게 계의 변화를 관찰할 수 있다.

간단한 예를 위해 에너지가 E(x)=x² 으로 정의되는 계에 놓여진 입자를 몬테칼로 시뮬레이션하자.

입자의 위치 x는 -10과 10 사이에서 시뮬레이션 한다. kT는 0보다 큰 실수이다.
아래는 Metropolis 알고리듬을 이용한 몬테칼로 시뮬레이션 방법이다.

1. 계를 초기화하고, 임의의 초기 위치 x를 선택한다.
2. x에 -1과 1사이의 임의의 난수를 더하여 입자를 임시로 이동시킨다.
• x' = x + 난수[-1,1]
3. 이동 전의 에너지 E(x)와 이동 후의 에너지 E(x')를 계산하여 에너지 변화 ΔE를 계산한다.
• E( x) = x*x
• E(x') = x'*x'
• ΔE = E(x') - E(x)
4. Boltzmann factor를 계산한다.
   • B = Exp[ -ΔE/kT ]
5. Boltzmann factor B와 0과 1사이의 난수를 이용하여 입자 이동을 허락할지 취소할지 결정한다.
   • 난수가 B보다 작으면 입자를 새로운 위치로 이동시킨다 : x=x'  
   • 난수가 B보다 크면 입자는 이동할 수 없다             : x=x
6. 2번으로 이동하여 반복한다.

E=x² 몬테칼로 시뮬레이션

시뮬레이션 C언어 소스 보기 플래시로 보기

그래프는 몬테칼로 시뮬레이션을 10⁶회 실행하여 얻은 입자의 분포곡선이다. X축은 입자의 위치이고 Y축은 그 위치에 입자가 위치한 횟수이다.

kT를 1 (빨간색)부터 9(파란색)까지 1씩 증가시키면서 비교하면, kT=9일 때 분포가 가장 넓게 퍼져있고 kT=1일 때 가장 좁고 뾰족한 분포를 볼 수 있다.

에너지 변화 ΔE가 음수이면 볼츠만 인자는 항상 1보다 크기 때문에 에너지가 낮은 쪽으로의 이동 (x=0 쪽)은 언제나 받아들여 진다. 따라서 분포는 x=0에서 언제나 최대이다.

ΔE가 양수이면, 즉 에너지가 높은 쪽으로의 이동은 볼츠만 인자가 1보다 작기 때문에 이동할 확률은 Exp[ -ΔE/kT ]로 주어진다. 위의 시뮬레이션 과정 5번이 바로 이 확률로 입자를 이동시키는 방법이며 Metroplis 알고리듬을 이용한 것이다.

kT, 즉 온도가 높으면 입자는 더 활발히 움직일 수 있다. 그래프에서 kT가 클 수록 더 넓은 영역으로 퍼지는 것을 확인할 수 있다.

시뮬레이션의 C언어 소스를 제공하니 직접 컴파일해 시뮬레이션 해보기를 바랍니다. 또한, 동일한 몬테칼로 시뮬레이션을 플래시로 구현하였습니다. kT에 따라 분포가 어떻게 변하는지 바로 확인할 수 있습니다.

리눅스에서 컴파일하고 실행하는 방법이다.
gcc -o mc mcx2.c -lm
-lm 은 math 라이브러리를 링크하기 위한 옵션으로 꼭 필요하다 (exp함수를 사용하고 있다). -o mc 는 실행파일 이름을 mc로 만들라는 옵션이다. 따라서 아래처럼 실행하면 histogram으로 시작되는 데이타 파일이 생성된다.
./mc

kT 값을 바꾸기 위해서는 7번째 줄의 값을 바꾸면 되며, 데이타 파일의 이름은 자동으로 바뀐다.
#define KbT 1.0

시뮬레이션 실행 횟수를 바꾸기 위해서는 5번째 줄의 값을 바꾼다. 횟수가 클 수록 깨끗한 곡선을 얻을 수 있다.
#define Trial 1000000

얻은 데이타 파일을 그래프로 그리기 위해 gnuplot을 이용한다. gnuplot을 실행한다.
gnuplot

데이타 파일을 읽어 그래프를 그린다. plot 명령을 사용한다.
gnuplot>plot "histogram_kbt1.00.dat"

추가로 파일을 읽어 그래프를 그린다. 기존 그래프에 replot으로 곡선을 추가한다.
gnuplot>replot "histogram_kbt2.00.dat"
gnuplot>replot "histogram_kbt3.00.dat"
gnuplot>replot "histogram_kbt4.00.dat"
...
gnuplot>replot "histogram_kbt9.00.dat"

* 다른 예 solid-on-solid 결정성장 모델 시뮬레이션 - C언어 소스 보기