NEWTON 함수를 사용하여 비선형 방정식의 해 구하기

오늘은 수치계산과 관련된 주제를 다뤄보기로 하겠습니다. IDL에서 비선형(Non-linear) 방정식의 해를 구하는 예제를 소개해볼까 하는데요. 여기서는 NEWTON이라는 내장함수를 사용하고자 합니다. NEWTON 함수는 그 이름에서 이미 짐작할 수 있듯이 수치계산 분야에서 얘기하는 Newton's Method 또는 Newton-Raphson Method라는 기법을 사용하여 방정식의 해를 구하는 역할을 합니다. 특히 계산이 좀 난감한 비선형 방정식의 해를 구하는데 있어서 유용하게 사용됩니다. 이 기법 자체에 관해서 제가 여기서 자세히 설명을 하지는 않겠습니다(구글링을 좀 해보면 많이 나옵니다). 다만 밑에서 설명하는 과정에서 필요한 내용이 있다면 간단히 언급하기로 하겠습니다.

설명을 위하여 먼저 다음과 같이 예제 방정식을 하나 가정하고 그림을 그려 봅시다. 방정식의 형태는 3차식을 가정하였고, 여기서는 다음과 같이 독립된 FUNCTION 프로그램의 형태로 만들었습니다. 이렇게 한 이유는 NEWTON 함수를 사용하는데 있어서 이러한 방식이 꼭 필요하기 때문입니다. 아래 예제코드에서는 정의된 함수를 사용하여 XY축 공간상에 플롯을 그리고, 이해를 돕기 위하여 Y=0인 수평선도 추가로 그리도록 하였습니다.

FUNCTION func_nonln, x

RETURN, (x-1.4)^3-5

END

PRO test_nonlinear

xp = [-10.:10.:0.1]

yp = func_nonln(xp)

win = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl = PLOT(xp, yp, THICK=2, YRANGE=[-10, 10], /CURRENT)

line = POLYLINE(pl.XRANGE, [0, 0], /DATA)

x = -5.

result = NEWTON(x, 'func_nonln')

PRINT, result

point = SYMBOL(x, 0, 'circle', /SYM_FILLED, SYM_COLOR='red', /DATA)

line = POLYLINE([result, result], pl.YRANGE, COLOR='red', /DATA)

END

위 내용을 test_nonlinear.pro라는 프로그램 파일로 저장한 후 컴파일 및 실행하면 다음과 같은 그림이 표출됩니다.

일단 이 그림을 보면 우리가 찾고자 하는 해가 대략 어떤 값인지는 짐작이 가능합니다. 곡선과 X축이 만나는 지점의 X좌표값을 찾으면 되는데 대략 3보다 약간 큰 값일 것으로 보입니다. 물론 정확한 계산을 위해서는 이제 NEWTON 함수를 사용하면 됩니다. 이를 위하여 다음과 같은 내용을 위 프로그램의 END 직전에 추가하면 됩니다.

x = 5.

result = NEWTON(x, 'func_nonln')

PRINT, result

여기서는 NEWTON 함수를 사용하였는데, 두 개의 인자들이 사용되었습니다. 첫번째 인자는 변수 x인데 그 값은 5.0으로 정의되어 있습니다. 이 값은 방정식의 해를 찾기 위하여 처음에 가정하는 일종의 출발점에 해당되는 값입니다. 사실 Newton-Raphson 기법에서는 임의의 포인트에서부터 방정식의 패턴에 따라 해에 가까운 값으로 서서히 근접하는 방법을 사용합니다. 이를 위해서는 처음 출발점은 있어야하기 때문에 이러한 값을 주게 됩니다. 어떤 값을 사용해도 상관이 없는데, 일단 여기서는 5.0으로 가정하였습니다. 그리고 두번째 인자는 따옴표로 문자값을 명시하였는데 바로 FUNCTION으로 정의된 함수의 이름입니다. 해를 구하고자 하는 방정식이 정의된 함수 부프로그램의 이름을 이렇게 두번째 인자로 줘야 합니다. 이 정도면 NEWTON 함수를 사용할 수 있으며, 그 결과를 result라는 변수에 담아 출력할 수 있습니다.

이대로 실행하여 출력된 값을 보면 아마 3.11 정도의 값이 나올겁니다. 우리가 앞서 그림상으로 대략 짐작했던 값과 비슷합니다. 좀 더 도식적으로 쉽게 확인하기 위하여 END 직전에 다음과 같은 내용도 추가해서 프로그램을 다시 실행해 봅시다. 그러면 다음과 같은 그림을 얻을 수 있습니다.

point = SYMBOL(x, 0, 'circle', /SYM_FILLED, SYM_COLOR='red', /DATA)

line = POLYLINE([result, result], pl.YRANGE, COLOR='red', /DATA)

여기서는 출발점 포인트의 위치를 붉은 점으로 표시하였습니다. 그리고 NEWTON 함수로 계산된 result 값을 사용하여 수직 방향의 선을 표시함으로써, 제대로 해가 구해졌는지를 눈으로 쉽게 확인해볼 수 있습니다. 이 정도면 NEWTON 함수를 사용하여 방정식의 해를 구하고 그 결과를 그림으로 표출하는 작업을 수행하는 프로그램의 골격은 완성된 셈입니다. 그리고 한가지 더 짚고 넘어갈 부분이 있는데, 바로 출발점의 위치입니다. 앞선 예제에서는 5.0을 사용했는데, 사실 다른 값을 사용해도 마찬가지입니다. -5.0을 사용했을 때 결과는 다음 그림과 같습니다.

물론 그렇다고 해서 어떤 경우든 출발점의 위치가 전혀 상관없는 것은 아닙니다. 방정식의 패턴에 따라서는 출발점에 따라 결과가 달라질 수 있습니다. 위의 프로그램에서 FUNCTION 부프로그램의 내용만 다음과 같이 변경해서 2차 곡선에 해당되는 패턴의 방정식으로 테스트해봅시다. 그리고 이번에는 출발점 x의 값을 8.0으로 가정하겠습니다.

FUNCTION func_nonln, x

RETURN, (x-3.7)^2-5

END

그 결과는 다음 그림과 같습니다. 해는 약 5.94 정도의 값으로 계산됩니다.

사실 이 방정식은 그림을 보면 알 수 있듯이 해가 둘입니다. 그리고 출발점의 위치에 따라 어느 해가 계산될지가 결정됩니다. 따라서 만약 출발점을 -5.0로 가정할 경우에는 그 결과는 다음 그림과 같습니다. 해는 약 1.46 정도의 값으로 계산됩니다.

그리고 출발점을 5.0으로 가정할 경우에는 다음 그림과 같은 결과를 얻게 됩니다. 출발점의 위치에 따라 결과가 달라지는 것을 쉽게 확인할 수 있습니다.

따라서 좀 더 복잡한 형태의 함수일수록 출발점의 위치에 따라 얻게 되는 결과값이 달라질 수 있다는 점을 유의할 필요가 있습니다. 다음과 같이 X축과 세번 교차하는 형태의 방정식으로도 한번 테스트해보시기 바랍니다.

FUNCTION func_nonln, x

RETURN, 4*x^3 - 5*x^2 - 8*x + 1

END

이 방정식을 사용한 결과 예제 그림을 첨부하면서 마치도록 하겠습니다.