다이내믹 시스템의 시작 - 미적분

보통 다이내믹 시스템을 말하면 미분방정식 (differential equation)을 이야기 하지만 넓게 보았을 때 미분과 적분도 다이내믹 시스템으로 볼 수 있을 것이다. 시간에 따라 변화하는 시스템, 이런 물리 모델을 수학 모델로 바꾸면 미분방정식(differential equation), 미적분(calculous)이 된다. 또한 이런 물리 모델이 역학이든 회로이든 상관없이 동일한 수학 모델로 바뀐다는 사실에 아이디어를 얻어 기계 공학적 문제를 회로 문제로 변환하거나 반대로 회로 모델을 기계 공학 모델로 변환하는 연구도 있다.

문과가 아닌 이과계열을 선택해 대학에 진학했다면 언제 어디서 맞닥뜨릴지 알 수 없는 것이 미적분이다. 벡터에서도 벡터의 미분과 적분을 다루고 미적분과 아무런 관련이 없어 보이는 선형대수(linear algebra)에서도 eigenvalue, eigenvector로 넘어 가면 discrete differential equation(이산 미분방정식)과 만나게 된다. 심지어 복소수의 미적분을 다루는 complex calculus도 있다. 미적분과 관련이 없어 보이는 경제, 환경 이론에도 미적분의 표기는 등장한다. 이런 광범위한 쓰임새와 중요성을 갖고 있지만 전문 수학자들이 아닌 공학도와 일반인들이 이해하기에는 미적분의 시작부터 머리가 아프다.

모든 미적분 책들은 함수 f(x)의 개념과 limit(극한)을 소개하면 시작한다. 문제는 이 limit 부터가 언뜻 보기에 미친 말장난 같이 보인다는 점이다.

(1)

식(1)을 계산하는 것이 미분을 이해하는 첫 걸음인데 이게 아무리 봐도 말이 안 되기 때문이다. 문제는 여기서 그치는 것이 아니다. 이 limit을 설명하는 도중에 어떤 경우는 이건 0이고 이건 0에 가깝게 접근하지만 0은 아니라고 한다. 어떤 것은 0이고 어떤 것은 0이 아니란 말인가? 라는 질문이 계속 떠 다니게 되는 것이다. 세상의 모든 일들이 다 마찬가지지만 숫자도 상대적 가치를 지닌다. 다음 식을 생각해 보자

(2)

식 (2)는 누구나 계산을 할 수 있다. 답은 1이다. 여기서 소수점 자리도 세기 어려운 0.000001은 대단히 중요한 의미를 지닌 숫자이다. 다음 식을 살펴보자.

(3)

이 식의 정확한 값을 구해 낼 수  있는가?  정확한 값의 계산은 계산기를 사용해도 입력하는데 시간이 걸릴 것이다. 하지만 우리는 위의 식에서 0.0000023456789라는 숫자를 무가치한 숫자로 무시해서 0으로 놓고 계산할 수 있다. 왜냐하면 계산에서 사용되는 10000이나 25라는 숫자에 비해 0.... 이란 숫자의 값은 대단히 작아서 결과 값에 거의 영향을 못 미치기 때문이다. 또 이런 이유로 큰 값과 아주 작은 값이 뒤섞인 행렬의 가우스 소거법은 round off error를 발생 시키기도 한다. 여기까지 이해했다고 하더라도 아직도 넘어야 할 산이 남아있다. 수학에서는 실제 수치를 다루지 않기 때문이다. 임의 숫자 n이라던가 다양한 수식으로 표현 되는데 그냥 보아서는 식(2)의 경우인지 식(3)의 경우인지 알 수 없기 때문이다. 식을 판별하기 위해 수학자들은 온갖 방법을 동원해 식을 변형한다. 바로 이 부분에서 미적분에 대한 관심이 사라지기 시작하는 것이다.

시작부터 어려운데 본격적으로 미적분에 들어가면 더욱더 어렵다. 미적분의 계산은 보통 잘 알려진 기본 함수 몇 개의 결과 값을 기본 계산 법칙을 이용해 이리 저리 변형해 계산해 내기 때문이다. 미분은 다음과 같은 기본 함수 미분과 4개의 기본 미분 계산 법칙이 근간을 이룬다.

---

---

---

대부분의 미분 문제는 위의 기본 함수 미분 공식과 미분 계산 법칙을 이용해 풀 수는 있는데 그렇게 단순하다면 미적분 책이 천 페이지 가량 되지는 않을 것이다. 이론적으로는 그렇지만 위 공식을 적용하는 문제가 그렇게 단순하지 않다.  상수 값 α, β 를 판별하는 것부터 단순하지 않다. 그래도 상수 값은 비교적 쉽게 찾을 수는 있지만 우리가 아는 미분 공식을 적용해 풀기 위해 함수 h(x), g(x)를 찾는 것은 상당한 연습을 필요로 한다. 일례로 tan(x)의 미분 값을 위의 공식을 이용해 찾아 보아라. 생각보다 만만치는 않다.

적분은 기본적으로 antiderivative 함수라서 미분 값을 안다면 그 미분 함수의 적분 값은 원 함수 값이란 것을 알 수 있다. 기본 적분 함수는 미분에서 유도된다. 적분도 미분의 계산 법칙과 유사한 선형성이 적용된다. 하나 다른 점이라면 미분과 유사한 직접적인 체인 룰(chain rule)이 존재하지 않기 때문에 복잡한 적분을 풀기 위해서는 미분보다 훨씬 더 복잡한 수학 테크닉을 필요로 한다.

이 글의 시작부터 지금까지 꽤 우울한 이야기를 계속했던 것 같다. 얼마나 고도의 미적분의 테크닉을 익혀야 하는지에 관해 각자의 견해가 있겠지만 python의 sympy pacakge나 maxima 같은 프로그램에서 심볼릭 연산을 할 수 있는 지금 굳이 고도의 미적분 테크닉을 익혀야 할 필요가 있는지는 의문이다. 다음과 같은 식을 생각해 보자

위의 식을 손으로 적분하는 것은 쉬운 일이 아니다. partial fraction이란 테크닉이 이용되는데 partial fracction에 이용되는 함수들의 적분 테이블들을 외우고 있어야 한다. 미분과 적분만을 공부하는 사람이라면 물론 외울 수도 있을 것이다. 하지만 배우고 익혀야 하는 기술들이 넘쳐나는 현대에서 그런 걸 외우고 저런 테크닉을 연습한다는 것이 큰 의미가 있을지는 생각해 볼 일이 아닌가? symbolic integral의 강자 maxima를 이용해 문제를 풀어 보자.

(%i1) f : (x^2-29*x+5)/((x-4)^2*(x^2+3));
                                  2
                                 x  - 29 x + 5
(%o1)                          -----------------
                                      2   2
                               (x - 4)  (x  + 3)
(%i2) integrate(f, x);                                                     
                                       x
                    2        2 atan(-------)
               log(x  + 3)          sqrt(3)                   5
(%o2)        - ----------- + --------------- + log(x - 4) + -----
                    2            sqrt(3)                    x - 4

자신의 컴퓨터에 maxima가 없더라도 웹에서 maxima 프로그램을 이용할 수 있다. Integral Calculator라는 사이트는 maxima를 위한 웹 인터페이스를 제공해 준다. 방문해서 이용해 보기 바란다. 적분하고 싶은 공식을 입력하면 적분 풀이 과정과 결과를 산출해 준다. 단 2줄의 코드로 끝낼 수 있는 것을 머리를 싸잡아 매고 손으로 풀어야 하는 이유가 있는가? 심지어 copy and past 할 수 있는 LaTeX 코드까지 생성해 준다.

(%i3) tex(%o2);
$$-{{\log \left(x^2+3\right)}\over{2}}+{{2\,\arctan \left({{x}\over{
 \sqrt{3}}}\right)}\over{\sqrt{3}}}+\log \left(x-4\right)+{{5}\over{x
 -4}}\leqno{\tt (\%o2)}$$ 

물론 그렇다고 analytical method를 배워야 할 필요성을 부정하는 것은 아니다. 분명 기본적인 적분 공식과 기호, 법칙 등은 분명히 공부해야 한다. 그래야 자신에게 필요한 공식을 스스로 찾아낼 수 있고 기본 공학 이론들을 이해할 수 있기 때문이다. 하지만 기본을 넘어선 과도한 문제 풀이와 연습은 개인적으로 시간 낭비이지 않는가 생각한다. 고도의 수학 테크닉보다 기본 개념의 충실한 이해가 공학 문제를 컴퓨터로 풀 때 훨씬 더 중요한 역할을 한다.

다이내믹 시스템을 얘기하면서 상당히 정적인 이야기만을 했다. 아주 다이내믹 하지는 않지만 그래도 위의 이야기 보다 상대적으로 약간은 더 다이내믹한 이야기를 해 보자. 미적분의 창시자 뉴턴은 미분을 개발한 뒤 미분을 이용해 함수의 근(root) 즉, 함수의 그래프가 x축과 만나는 정확한 점을 찾는 방법을 고민했었다. 함수의 root를 찾는 것은 지금도 중요한 문제이고 대부분의 경우 analytical 방법으로는 함수의 근을 찾을 수 없기 때문에 다양한 numerical method가 개발 되었다. 하나의 예로의 정확한 값을 계산해 보자 이 값은 의 해이다. 또는 y=x 직선과 y=cos(x)의 그래프가 만나는 지점은 어디일까? 그것은 x-cos(x)=0 방정식의 근이다. 하지만 이런 실수 값을 계산할 analytical method는 없다. 여기에 numerical method가 등장한다. 뉴턴 라프슨(newton-raphson) 알고리듬은 간단하고 대단히 빠르게 수렴하는 장점이 있어서 단변수(single variable) 함수의 근을 찾는 상당히 효과적인 방법이다.

초기 값 x0의 함수 값 f(x0)에서 접선을 그렸을 때 그 접선과 x축이 만나는 점을 x1이라고 하자. 함수 f(x)의 미분 f'(x)는 바로 저 접선의 기울기이다. 접선의 기울기 m은 위의 기하학적 모델에서 유도하면

따라서 m = f'(x0) 이다.

바로 이 공식이 뉴턴 라프슨(newton-raphson) 알고리듬이다. 현재의 x 값에서 함수 f(x) 값이 0에 가까워지는 다음 x 좌표 값을 반복적으로 계속 구하는 방식이다. 이 과정을 그래프로 기록하면 다음과 같다.

자동으로 0에 수렴하는 x 값을 찾아가는 과정이 상당히 흥미롭다. 하지만 newton method도 만능은 아니라서 특정 함수에서 잘못된 초기값 선택은 알고리듬이 수렴하지 못하고 계속 같은 자리를 멤돌 수 있다. 또한 미분을 사용하기 때문에 연속인 함수에서만 이용할 수 있다는 점에 주의를 기울여야 한다. 위의 그래프를 만드는데 maxima 프로그램을 이용했다. 아래는 그 코드이다.

f: x^3-6*x^2+12*x-6$
/* starting x value */
s_val: 0.5$
/* margin of x coordinate in graph */
pad: 0.2$
/* store the history of caculation by newton raphson method*/
array(y, 10)$ /* store y coordinate */
array(x, 10)$ /* store x coordinate */ 
array(slope, 10)$ /* store the derivative of f(x) */

/*
newton raphson method
======================
f      : equation
s_val  : starting x coordinate
return : iteration count
*/
newton_root(f, s_val) :=
block ([cnt: 0, df: diff(f,'x,1), val: s_val],
    display(f),
    display(df),
    printf(true,"| ~2a | ~10a | ~10a | ~%","id", "x", "y"),
    for i: 0 while (i < 10 ) 
        do(
        slope[i]: float(at(df, x=val)),
        y[i]: float(at(f, x=val)),
        x[i]: val,
        val: val - y[i]/slope[i],
        printf(true, "| ~2d | ~10f | ~10f | ~%",i, x[i], y[i]),
        cnt: i,
        if (abs(y[i]) < 0.00001) then return (cnt)
        ), 
    cnt 
)$

cnt: newton_root(f, s_val)$

/* find minimum x coordinate value for graph */
xmin: lmin(
    listarray(
        rearray(x,cnt)
    )
)$

/* find maximum x coordinate value for graph */
xmax: lmax(
    listarray(
        rearray(x, cnt)
    )
)$

/* find minimum y coordinate value for graph */
ymin: lmin(
        [lmin(
            listarray(
                rearray(y,cnt)
            )
        ),
        at(f, x=xmin-pad),
        at(f, x=xmax+pad)]
    )$

/* find maximum y coordinate value for graph */
ymax: lmax(
      [lmax(
            listarray(
                rearray(y,cnt)
            )
        ),
        at(f, x=xmin-pad),
        at(f, x=xmax+pad)]
    )$

/* tangent equation */
line[m, x0, y0] := m*(x - x0) + y0$

legend1: printf(false,"~a",f)$
for i: 0 thru cnt do (
    /* 
        create file name for graph 
        newton#.png
        ex) newton0.png
    */
    fname: concat("/home/main/newton", i, ".png"),
    legend2: printf( false, "~6f * ( x - ~6f ) + ~6f",slope[i], x[i], y[i] ),
    plot2d(
        [f, line[slope[i], x[i], y[i]]],
        [x, xmin-pad, xmax+pad], [y, ymin, ymax],
        [legend, legend1, legend2],
        [plot_format, gnuplot], 
        [gnuplot_term, png],
        [gnuplot_out_file, fname]
    )
)$
kill(all)$

코드가 조금은 길지만 사실 뉴턴 알고리듬을 구현한 코드는 굵게 표시한 부분이 전부이다. 나머지는 과정을 text로 출력하고 그래프를 animated gif로 만들기 위해 다수의 그래프 파일을 생성하는 코드인데 배보다 배꼽이 더 커졌다.

| id | x          | y          | 
|  0 |        3.5 |      5.375 | 
|  1 |  2.7037037 |  2.3484733 | 
|  2 | 1.12287541 | 1.32518635 | 
|  3 |  0.5487152 | -1.0567363 | 
|  4 |  0.7159551 |  -0.117096 | 
|  5 |  0.7396285 |  -0.002145 | 
|  6 |  0.7400788 | -0.0000007 |  

6번의 연산만으로 빠르게 f(x)값이 0가 되는 지점을 찾아가는 것을 볼 수 있다. 실험해 본 바로는 초기 값이 적당하다면 대부분 4~6번 연산에서 root 값을 찾았다.

위의 코드에 조금은 문제가 있는데 항상 사람이 적당한 초기 값을 지정해 주어야 하고 함수의 근이 다수인 경우 자동으로 다 찾아 주지 않는다. 이런 문제를 해결하기 위해 인간이 초기 값을 찾는 방식과 유사한 방법을 적용해 봤다. 인간은 초기 값을 함수의 그래프를 관찰한 뒤 지정한다. 마찬가지로 컴퓨터도 뉴턴 알고리듬을 실행하기 전 먼저 그래프에서 함수 값 f(x)부호가 바뀌는 근사점을  찾아 초기 값으로 자동 지정하게 프로그램 해주면 된다. 이것을 octave 프로그램을 이용해 구현했다. 특히 아래 octave 프로그램은 함수 미분을 analytical method가 아닌 numerical method로 계산한다. 미분을 numerical method로 구현하는 건 미분의 analytical method가 무지 무지 복잡한 것에 비해 완전 중학교 수준의 수학으로 간단하다.

#script

% clear a memeory
clear -a

% adjust display precision of floating point
format long

% call back function for close event of plot window  
function close_function()
    disp('closing...')
    closereq()
    exit()
endfunction 

function [result, iter] = newtroot (fname, x)

% usage: newtroot (fname, x)
%
%   fname : a string naming a function f(x).
%   x     : initial guess

    %tol  :tolerance value  
    %deta : dx
    delta = tol = sqrt (eps);
    % max iteration
    maxit = 20;
    fx = feval (fname, x);
    for i = 1:maxit
        if (abs (fx) < delta)
            result = x;
            iter = i;
            return;
        endif

        % f(x)
        fx_new = feval (fname, (x + delta)) ;

        %derivative f(x)
        deriv = (fx_new - fx) / delta;

        %newton raphson algorithm
        x = x - fx / deriv;
        fx = f(x);
    endfor

    iter = i;
    result = x;
endfunction

function ret = f(x)
% Test equations
%    ret = x.**7 + 3*x.**6 + 17*x.**5 - 57*x.**3 -10*x.**2 + 29*x +20;
%    ret = x.**3 - 3*x.**2 +2;
    ret = x.**3-6*x.**2+12*x-6;
%    ret = x.**3 - cos(x);
%    ret = x.**2 - 3;
%    ret = exp(-1*x).*(3*cos(2*x)+4*sin(2*x));
%    ret = x - 2*sin(x);
% ret = x.**2 - 1;
endfunction

% more memory efficient than linespace(-2:2:40) 
x= 0:0.1:3.5;
y=f(x);

% number of row and col 
[nr, nc] = size(y);
idx=1;

% find the closest guess and caculate newton roots
printf("root = \n");
printf("%10s | %15s | %20s\n", "iteration", "eitamate", "newton root");
i = 2;
while (i <= nc)
    % exact root
    if ( y(i) == 0 )
        root(idx) = x(i);
        printf("%10d | %15g | %20.15g\n", 0, x(i), x(i) );
        i++;
        idx++;
    elseif( sign(y(i-1)) != sign(y( i )) )
        % caculate a newton root
        [root(idx), iter] = newtroot("f", x(i));
        printf("%10d | %15g | %20.15g\n", iter, x(i), root(idx) );
        idx++;
    endif
    i++;
endwhile


% plot the equeation
plot(x,y, "linewidth", 1.6)
hold on

% auxillary graph for x.**3 - cos(x) 
% plot(x, x.**3, "color", "c", "linewidth", 2)
% plot(x, cos(x), "color", "g", "linewidth", 2)

% if idx <= 1 then no closest guess of root, 
% and therefore no caculation of newton roots
if( idx > 1) 
    % mark roots
    plot(root, f(root), 
        "linestyle", "none", 
        "marker", "o",
        "markeredgecolor", "r",
        "markerfacecolor", "r",
        "markersize", 10)
endif
hold off
grid on

% saving a plot in a png format
print -dpng figure1.png

% register call back function 
set(gcf, 'closerequestfcn', @close_function)

% infinit loop for vi shell command execute
% wait for closing the plot window 
% exit at window close event
ginput ()

위의 코드도 다른 잡다한 기능들 때문에 길지만 실제적인 알고리듬을 구현한 코드는 사실 몇 줄 안 된다. 위에 굵게 표시한 부분이 numerical 미분을 계산하는 곳이다. delta 값으로 사용한 eps 값은 머신의 최소 숫자 단위이다. 컴퓨터 환경에 따라 다르지만 내 컴퓨터의 경우는 2.2204e-16이다. 이 값을 square root 했기 때문에 delta 값은 1.4901e-08이다.이니 dx가 0에 근접한다는 미분의 정의에 적당하지 않는가? 이 코드의 출력 값은 아래와 같다.

root = 
 iteration |        eitamate |          newton root
         4 |             0.8 |    0.740078950062632

그래프로부터 자동으로 0.8 이란 root 근사값을 찾았고 newton method로 0.740078950062632이란 root를 찾았다.  근이 하나가 아닌 여러 개인 함수  의 근을 찾아보자.

root = 
 iteration |        eitamate |          newton root
         5 |            -1.7 |    -1.75403724912773
         5 |             1.1 |     1.06292835927139
         6 |             1.4 |     1.30913820115207

총 3개의 근사값과 3개의 newton root를 찾았다. 총 16번의 연산으로 3개의 루트를 찾았으니 상당히 효율적이지 않은가? 이걸 사람 손으로 푼다는 건 비효율의 극치이다. 지금 로그 테이블을 손으로 구하고 그 로그 테이블로부터 로그 값을 구하는 사람이 없듯이 더이상 미분과 적분의 해도 사람 손으로 구하는 시대는 지났다. 미분방정식으로 넘어가면 아예 analytical method로 해를 구할 수 없는 경우가 허다하다. 개념에 충실하고 방정식을 구하는데 집중하자. 나머지 계산은 이제 컴퓨터가 알아서 해 준다.