적분을 수행하는데 여러가지 방법이 있지만, 어떻게 해야 할지 정말 모르겠다면 몬테 카를로 방법이 있다.


몬테 카를로 방법은 단순히 면적을 계산하는 방법인데 알고리즘은 다음과 같다.


1. 난수 (x, y)를 생성한다.

2. f(x)>y이면 카운트를 +1한다.

3. N번 반복한 후

4. 카운트를 N으로 나누면 적분값이다.


x의 범위는 적분 구간, y의 범위는 f(x)가 해당 적분 구간에서 갖게 되는 최대값과 0 사이의 영역이다.


잘 보면 알겠지만, y가 f(x)와 0으로 둘러싸인 구간에, 임의로 점을 뿌리는 과정이다.


점을 다 뿌린 후, 뿌린 수 중에 몇개나 구간 안에 들어갔는지 갯수를 센다. 그럼, 해당 함수로 둘러싸인 영역 안에 들어갈 확률은 영역의 넓이에 비례하므로, 뿌린 수 대비 들어간 수의 비율을 계산하면 된다.


실제 코드는 다음과 같다.


import random as rd
import numpy as np
import os

fun = lambda x:np.sqrt(1.-x*x)
xA = 0. # x의 구간 시작
xB = 1. # y의 구간 끝
yA = 0. # f(x)의 최소값
yB = 2. # f(x)의 최대값보다 큰 임의의 값. 클수록 정확해지고 작을수록 빨라진다. f(x)의 최대값보다 같지 않고 더 큰 수를 넣어야 한다.
AREA = np.abs((xA-xB)*(yA-yB)) # 구간으로 둘러싸인 영역의 넓이.

rd.seed(os.times())

it = 1000

savefile = open("circle.txt", "w")

counted = 0

for i in range (it):
    if fun(rd.uniform(xA, xB)) > rd.uniform(yA, yB): counted+=1

print str(it)+" iteration, "+str(AREA*float(counted)/(float(it)))+" is integration result.")


빨간색으로 색칠한 부분이 "핵심코드"이다. 너무 짧은거 아니냐고 물어볼수도 있지만, 진짜로 저게 끝이다.


여기서, 함수 안에 들어간 경우의 수를 반복한 횟수로 나눠준 다음 왜 면적을 곱해주는 걸까? 생각해 봅시다.


0부터 1까지 Sqrt(1-x*x)를 적분한 결과. 가로축은 반복 횟수인데, 백만번정도 반복한 것까지 그렸다. 잘 보면 정확히 pi/4에 수렴하는 것을 관찰할 수 있다. 


아무생각없이 짜도 잘 굴러가는 굉장히 아름다운 알고리즘이다.


문제는 f(x)가 음수인 경우에는 아무생각없이 짜면 안되고 생각은 한번 해줘야 한다는 점.

임의의 N차원에서 수행하는 다중적분인 경우에는 난수를 생성하는 구간을 "잘" 잡거나, 또는! 그냥 "충분히 큰" 초직다면체(hyper-rectangular)에 해당하는 난수를 생성한 후 갯수를 셀까 말까 카운트하는 판정 루틴에 f(x)보다 작아야 한다는 조건 뿐만 아니라 "내가 원하는 구간 안에 들어가야 함"까지 조건을 넣으면 된다. 물론 그게 그거겠지만, 여기서는 단순히 "안에 있는가 없는가"만 판단하면 되기 때문에 그냥 다중적분보다는 쉽다.


다시말해서, 내가 어디에서 뭘 적분하는지만 알고 있으면 수치적인 해를 구할 수 있다는 뜻이다. 


...


여담(=개드립)인데.


인생은 내가 어디서 뭘 해야 하는지 알지도 못할 뿐만 아니라, Monte carlo 방법을 쓴답시고 도박하다간 망한다. 수치해석적 인생은 수치스러운 인생으로 귀결된다.

by snowall 2013. 10. 15. 01:57