Simpson's Rule 은 적분이 들어가는 모든 문제에 캡숑 유용하게 쓰이는 방법이다. 목적함수를 2차함수로 근사하기 때문에, 실제 목적함수의 차수가 2차 이하의 다항식인 경우에는 귀찮게 손으로 적분할 필요 없이 정확한 답을 구할 수 있다. 당연하지만 이 이외의 경우 -목적함수가 다항식이 아니거나 한 경우-에는 오차가 매우 커진다.
알고스팟 모의고사 준비하면서 sqrt(r^2 - x^2) 를 적분할 일이 있었다. 물론 닥치고 울프람 형아 적분헤 주떼영 하면 되지만 나도 남자의 자존심이 있어서 그건 싫고 -_- 해서 numerical method 를 쓸 양으로 Simpson's rule 을 써 보았다.
lang:py
def f(x):
return sqrt(1 - x*x)
def simpson(a, b):
c = (a+b)/2.0
h3 = abs(b-a) / 6.0
return h3 * (f(a) + 4.0*f(c) + f(b))
하고, -1 부터 1 까지 적분하면 반원의 넓이니까 pi/2 가 나와야겠지.
In [2]: print simpson(-1,1)
------> print(simpson(-1,1))
1.33333333333
깨갱 -_- 0.2나 차이나네; 이렇게 되면 질수없다. 인터벌을 1000개로 나눠서 각각 구간에서 simpson's rule 을 돌리기로 한다.
In [9]: xs = arange(-1,1,0.002)
In [10]: sum([simpson(xs[i],xs[i+1]) for i in range(len(xs)-1)])
Out[10]: 1.5707083802290707
오 꽤 정확해짐. 하지만 아직도 오차 범위가 1e-5 라서 많이 부족하다. (대개 프로그래밍 대회에서는 1e-7 이하의 정확도를 요구한다.) 인터벌을 몇개로 늘리면 될까. 1만개로 테스트.
In [13]: xs = arange(-1,1,0.0002)
In [14]: pi/2-sum([simpson(xs[i],xs[i+1]) for i in range(len(xs)-1)])
Out[14]: 2.7818300534221407e-06
아니 시발 이게 뭐야 -_-; 결국 테스트를 해보니 5만개 정도로는 나눠야 1e-7 정확도를 얻을 수 있다. 물론 일루옹이 낸 문제답게 n 이 작지도 않고.. 이대로는 절대 시간 안에 안나온다. 그래서 써먹은 것이 Adaptive simpson's method.
기본 규칙은 매우 간단하다. 이 함수가 생긴 모양을 봐서 2차 함수로 적절히 근사가 된다 싶으면 그냥 simpson's rule 을 써서 근사하고, 아니면 절반으로 나눠서 divide & conquer. 위키피디아 엔트리 링크
lang:py
def adaptive(a, b, sum, eps):
c = (a+b) / 2.0
left = simpson(a, c)
right = simpson(c, b)
if abs(left + right - sum) <= 15 * eps:
return left + right - (left + right - sum) / 15
return adaptive(a, c, left, eps/2) + adaptive(c, b, right, eps/2)
15 로 나누는 부분이 매우 맘에 걸리지만 지금은 수치해석시간이 아니니 자세한 분석은 넘어가기로 하고.. (나 수치해석 학점 B 다..) 말할 것 없이 돌려보면
In [31]: pi/2 - adaptive(-1, 1, simpson(-1, 1), 1e-8)
Out[31]: 7.0505905558349014e-09
실제 simpson() 함수가 계산되는 회수를 세어보면 800회정도밖에 안된다. 우훗! 물론 이거야 이 함수가 매우 스무스하게 생긴거라 그러니 오해하면 안됩니다. 여튼 적분하기 싫은 여러분, 이렇게 이겨냅시다 [....]
