IDL로 연립방정식을 풀어봅시다

IDL의 기본 라이브러리에는 매우 많은 종류의 수치계산(Numerical Calculation) 루틴들이 존재합니다. 너무 많아서 사실 저같은 사람도 그 많은 루틴들 중 직접 사용을 해본 것들은 그리 많지는 않습니다. 어떤 종류의 루틴들이 있는가는 IDL 도움말에서 Contents - Routines(by Topic) - Mathematic 섹션을 보면 잘 나와 있습니다. 특정한 수치계산 알고리즘을 필요로 하는 작업을 하는데 있어서 IDL의 어떤 루틴을 사용하면 되는가에 대해서는 이 섹션의 내용을 참조하면 됩니다. 오늘은 연립방정식을 푸는 루틴들 몇 가지만 간단하게 소개해보겠습니다.

먼저 예제로서 다음과 같은 3원 1차 연립방정식을 푼다고 가정해봅시다. 이와 같은 연립방정식은 결국은 미지수의 계수들만으로 이루어진 A라는 행렬, 미지수들로 이루어진 X라는 행렬, 그리고 방정식의 우변의 값들에 해당되는 B라는 행렬로 단순화하여 표현됩니다. 결국 X를 구하는 것이 최종 목적이 됩니다.

위의 연립방정식을 실제 IDL 코드상에서는 다음과 같이 A와 B라는 배열로 먼저 정의해주면 됩니다.

A = [[2.0,  1.0, 1.0], $

        [4.0, -6.0, 0.0], $

        [-2.0, 7.0, 2.0]]

B = [3.0, 10.0, -5.0]

이와 같은 연립방정식을 푸는 가장 일반적인 방법은 Gauss 또는 Gauss-Jordan 소거법이라 할 수 있습니다. 하지만 IDL에 내장된 관련 루틴들은 수치계산적인 알고리즘들을 사용을 합니다. 주로 Numerical Recipes라는 수치계산 분야의 바이블이라 할 수 있는 이 책의 내용을 기반으로 한 것들이 많습니다. 어차피 세부적인 기법에 관한 설명은 여기서는 하지는 않겠습니다(어차피 하기도 힘듭니다). 그냥 어떤 루틴을 쓰면 내부적으로 어떤 알고리즘을 사용하여 해를 구하는가에 대한 소개만 드리려고 합니다.

먼저 Cramer's Rule이라는 기법을 사용하는 방법인데요. 이 경우에는 CRAMER라는 내장함수를 사용합니다. 다음과 같은 방식으로 사용하면 됩니다.

result = CRAMER(A, B)

PRINT, result

사실 복잡한 작업들은 내부적으로 다 이루어지기 때문에 함수 사용법 자체는 매우 간단합니다. 또 다른 기법으로는 LU Decomposition이라는 것도 있는데, 이 기법을 사용하는 LUDC 프로시저 및  LUSOL 함수도 있습니다. 사용법은 다음과 같습니다.

LUDC, A, tmp

result = LUSOL(A, tmp, B)

PRINT, result

두 경우 모두 출력된 결과값들을 보면 1.0, -1.0, 2.0이라는 값들로 산출되었음을 확인할 수 있습니다. 즉, 이 값들은 미지수였던 x1, x2, x3에 해당되며, 바로 위의 연립방정식의 해(Solution)입니다. 물론 이와 같이 유일한 해가 존재하는 평이한 연립방정식은 몇몇 내장함수들을 사용하면 아주 간단하게 풀립니다.

사실 연립방정식의 해를 구하는 가장 직접적인 방법은 A의 역행렬을 B에 곱해주는(행렬곱) 것인데요. 물론 이러한 방법도 사용할 수 있습니다. 이 경우에는 역행렬을 구해주는 INVERT라는 내장함수를 사용하여 다음과 같이 해를 구하면 됩니다.

invA = INVERT(A)

result = invA ## TRANSPOSE(B)

PRINT, result

물론 이와 같이 역행렬을 사용하는 방법에서는, 역행렬의 존재여부를 판가름해주는 Determinant라는 값이 유효한 값으로 존재해야 한다는 단서는 붙습니다. 만약 역행렬을 구하기에 앞서서 이 determinant의 존재여부를 먼저 확인하려면 DETERM이라는 내장함수를 A에 적용하여 의미있는 값이 나오는가를 먼저 확인하면 됩니다. 그리고 산출된 역행렬을 B에 행렬곱을 취할 때, 그냥 B를 바로 사용하지 않고 TRANSPOSE(B)를 사용하였는데요. 이유는 처음에 B라는 배열을 정의할 때 가로방향 배열로 정의가 된 상태이기 때문입니다. 실제로 해를 계산하기 위하여 A 역행렬을 B에 곱할 때에는 B가 세로방향 배열이어야 합니다. 이러한 변환, 즉 행렬의 원소값들의 방향을 대각선 기준 대칭으로 바꿔주는 역할을 하는 함수가 TRANSPOSE라는 사실도 참고해두시기 바랍니다.

연립방정식의 세계는 사실 여기까지가 다는 아닙니다. 미지수가 훨씬 더 많은 연립방정식의 경우에 대해서는 또 다른 기법들이 사용되기도 하는데요. Singular Value Decomposition, Cholesky Decomposition, Gauss-Seidel Iteration 등의 기법들을 지원하는 내장함수들도 마련되어 있습니다. 앞서 언급된 IDL 도움말에서의 해당 섹션에 이러한 기법들에 해당되는 내장함수들에 대한 설명 및 사용법이 나와 있다는 점 참조하시면 좋을 것 같습니다.