IDL에서 적분을 해봅시다 (Part 1)

오늘은 IDL에서 수치적분(Numerical Integration)을 수행하는 방법에 대하여 소개해볼까 합니다. 예를 들어 다음과 같은 2차 곡선이 있고, 일정한 X값 구간에 대한 적분값을 산출하는 상황을 보기로 하겠습니다. 여기서는 함수의 형태가 Y = X^2+4의 형태이며, 적분 구간은 2에서 6까지로 가정하겠습니다. 결국은 그림상에서 붉은 선으로 그어진 영역에 대한 면적을 구하는 셈입니다.

IDL에서는 다양한 기법으로 수치적분을 수행하는 루틴들이 지원됩니다. 위의 예와 같이 함수의 패턴이 정해진 상태에서 원하는 구간에 대한 적분을 수행하고자 하는 경우에는 QROMB, QSIMP와 같은 내장함수들을 사용하면 됩니다. 전자는 Romberg Integration, 후자는 Simpson's Rule이라는 수치적분 기법을 사용합니다. 각각의 기법에 대한 자세한 설명은 여기서는 생략하겠습니다. 다만 수치적분의 속성 자체가 전체 적분구간을 세부적인 구간들로 잘게 쪼개고, 쪼개진 세부구간별로 사다리꼴 또는 2차곡선 등의 근사된 형태 기준의 면적을 구한 다음, 이 면적들을 총합하여 전체 구간에 대한 적분값을 산출하는 방식이라고 간단하게 언급하면 될 것 같습니다. 어쨌든 위의 함수들을 사용하여 적분값을 구하는 예제코드는 아래와 같습니다.

FUNCTION INTFUNC, x

  value = x^2+4

  RETURN, value

END

PRO test_integration

  la = 2.

  lb = 6.

  result_rom = QROMB('intfunc', la, lb)

  result_sim = QSIMP('intfunc', la, lb)

  

  myform = '(A, F12.6)'

  PRINT, FORMAT=myform, 'QROMB Result with Single Precision = ', result_rom

  PRINT, FORMAT=myform, 'QSIMP Result with Single Precision = ', result_sim

  

END

유의할 점은 이 함수들이 필요로 하는 첫번째 인자 함수형 부프로그램의 이름이라는 것입니다. 즉 내가 적분하고자 하는 함수의 형태를 FUNCTION 부프로그램으로 정의를 해놓고 이 이름을 사용해야 한다는 점을 염두에 둬야 합니다. 나머지 두번째 및 세번째 인자값은 구간의 경계값을 넣어주면 됩니다. 사용법 자체는 QROMB나 QSIMP 모두 유사합니다. 출력된 결과값을 보면 둘 다 85.333336으로 나올겁니다. 참고로 이 함수들에는 /DOUBLE이라는 키워드가 지원됩니다. 말 그대로 연산을 2배 정밀도로 수행하여 적분값을 산출해주도록 하는 옵션입니다. 그런데 사용을 해보면 결과값의 차이는 소수점 4~5번째 정도나 가야 겨우 날 정도입니다. 따라서 그냥 일반 정밀도로 계산을 해도 큰 지장은 없을 것 같습니다.

그러면 이 값은 정확한 적분값일까요? 이 결과값에 대한 검증은 적분의 대상이 되는 함수의 원천함수(Integrand)를 구해서 수학적으로 직접 계산할 수 있습니다. 아마 다들 손으로 미적분을 계산해보신지는 꽤 오래 되셨을겁니다. 그래도 더듬더듬 기억을 되짚어서 해보면, 위 함수의 원천함수는 (X^3)/3+4*X로 나올겁니다. 이 원천함수에 lb, la를 대입해서 빼주면 됩니다. 제가 해본바로는 정확하게 일치하는 것으로 확인이 됩니다. 여러분들도 해보시면 아실겁니다.

따라서 위에서 소개된 QROMB, QSIMP와 같은 내장함수를 사용하면, 다양한 형태의 함수들에 대한 적분이 가능합니다. 그런데 적분의 대상이 되는 함수의 형태를 알 수 없는 경우도 있습니다. 예를 들면, 관측 데이터를 일정한 시간간격으로 얻었는데 그 변동패턴을 특정한 analytic 함수의 형태로는 정의하기 어려운 경우들이 꽤 있습니다. 이런 경우에도 방법은 존재합니다. 다음 회에는 이러한 경우의 수치적분 방법에 대하여 살펴보기로 하겠습니다.