matlab의 ode45 : 미분방정식의 해를 구하는 함수

   미분된 식이 존재할 때, 미분되지 않은 원래의 함수를 계산하는 기능이 matlab 의 ode45 로 제공되고 있다. ode45의 사용법을 살펴 본다.

ode45는 다음과 같은 형식으로 사용한다.

ode45(@미분된 함수 이름, 시간 구간, 초기값)

시간 구간은 row 벡터이고, 초기값은 column 벡터이다. 즉,

시간은 [0 10] 과 같은 형식이고,
초기값은 [1; 2] 와 같은 형식

이다. 미분된 함수의 이름^[각주:1]이란 사용자가 직접 작성한 함수로, matlab 내부에서는 함수 주소를 입력받아 그 함수를 호출하는 형식으로 진행된다. 일반적으로 우리는 다음과 같은 형식으로 문제를 접하게 된다. 즉, 

각 변수의 미분된 함수식만을 알고 있는 상황이며, 이 상황에서 미분되지 않은 함수를 찾아 내야 하는 것. 사용자가 작성해 넣어야 할 함수는 벡터를 반환하는 함수로, 위에서 f_i 를 각각의 요소로 반환을 하면 된다. 실제 예를 하나 보자. 모델은 다음과 같다.

위 모델은 S-system 으로 나타내기는 조금 까다로운데, 다음과 같은 간단한 상황이다.

Y1 : Y1의 생성은 S의 속도로 일정하게 계속 이루어지고 있다. 이 상황에서 Y2로의 변환은 alpha의 속도로 진행되고 있다. 만약 Y2가 충분히 만들어지면 beta의 정도로 Y1에서 Y2로의 변환이 Y2 자체에 의해 억제된다. 즉, Y2는 자신에 대한 negative feedback 을 이루고 있다.
Y2 : Y1에서 alpha의 속도로 생성되고 있다. 또한 생성된 Y2는 d 의 속도로 분해가 일어난다.

위와 같은 상황은 다음과 같이 모델링이 가능하다.

즉, 각각의 요소의 '변화량' 중에서 '생성'되는 것과 '소멸'되는 것을 따로 생각한 다음 그 두 식을 더하면 위와 같이 표현이 가능한 것이다^[각주:2]. 이 상황에서 우리는 y1'과 y2' 를 구현해서 ode45로 넘겨 주어야 한다. 다음과 같이 구현 가능할 것이다.

% negative_feedback.m 파일
function dydt = negative_feedback(time, y)
global alpha beta S deg;

dydt = zeros(2,1);
dydt(1) = S - alpha*y(2) + beta*y(1);
dydt(2) = alpha*y(1) - beta*y(2) - deg*y(2);

아주 간단한다. 이 때, alpha/beta/S/deg 변수는 외부에서 변화시킬 수 있도록 global 변수로 선언을 해 놓자. 즉, 우리는 각각의 변수를 일정한 값으로 시뮬레이션 한 후, 그 중 어느 값을 변화시켜서 다시 시뮬레이션 하는 식으로 여러 값들이 시스템에 어떤 영향을 주는지 보아야 하므로 모델에 사용된 변수들은 외부에서 접근 가능해야 할 것이다. 위처럼 negative_feedback.m 파일을 만들자. 이 때, 이 파일이 저장된 경로를 matlab의 working directory 로 변환시키거나 matlab 환경 변수 중 path 에 이 경로를 추가해야 한다. 그렇게 하기 위해서는 matlab 메인 메뉴의 File -> Set Path... 를 선택해서 작업한다. 이제 직접 사용하는 부분을 보자.

% process.m
clear;
clf;
global alpha beta S deg
alpha = 0.1;
beta = 0.03;
S = 0.5;
deg = 0.01;
hold on;
for i = 0:8
    alpha = i/8*0.25;
    [t y] = ode45(@negative_feedback,[0 10], [10; 0]);
    subplot(3,3,i+1);
    plot(t,y(:,1),'r',t,y(:,2),'g');
    title(sprintf('alpha=%.2f',alpha));
end;
hold off;

가장 중요한 부분은 붉은 색 글자 부분이다. 주어진 각각의 변수 값에서, 위처럼 alpha 값을 변화시켜 가면서 y1 과 y2 가 어떻게 변하는지를 그림으로 그리고 있다. 시간은 [0 10] 에서 알 수 있듯이 10 까지 하게 해 놓았고, 초기값은 y1은 10, y2는 0 으로 되어 있다. 나머지 부분의 자세한 설명이나 기본적인 matlab 사용법은 여기서 설명하지 않는다. 위처럼 해서 process.m 을 실행시키면 다음과 같은 그림이 생성된다.

그래프에서 붉은 색이 y1 이고, 초록색이 y2 이다. y1에서 y2 로의 변환 속도가 0, 즉 변환하지 않으면 y1만 계속 증가하는 것을 알 수 있다, 왜냐 하면 y1 은 S 의 속도로 계속 만들어지고 있다고 모델링 해 놓았기 때문이다. 보면 알 수 있듯이, 변환 속도에 따라 y1과 y2의 시간에 따른 양의 전체적인 양상이 다르다는 것을 알 수 있다. 위 코드에서 다른 변수들의 지정된 값(내가 임의로 넣었다)에서는, alpha가 0.13만 되어도 y1의 계속적인 증가는 일어나지 않는다는 것을 알 수 있다.

   y_i' 가 각 문제마다 다를 것이므로, 위의 경우에 negative_feedback.m 과 같은 함수^[각주:3]를 각 문제에 맞게 구현한 이후 위처럼 ode45를 이용하면 될 것이다.


ps. 유입 검색어 때문에 작성함. 미분방정식에 대한 수치해석학적 해에 대한 원리와 실제 구현 및 사용에 관한 글을 작성했기 때문인지 계속 matlab ode45 로 유입이 되고 있는데, 한두번이면 모르겠지만 꾸준히 들어와서 결국 작성한다. 난 proof-of-concept 의 확인을 위해 빠르게 해보아야 할 때만 다른 프로그램(maple이나 matlab, spss 등)을 이용하고, 일단 된다 싶으면 대부분 구현을 하기 때문에 각 프로그램을 그렇게 능숙하게 사용하지는 못한다.

1. 원래는 함수 '주소'이다. 그래서 함수 이름 앞에 @ 이 붙는다. [본문으로]
2. 물론 Y2가 자신의 생성을 억제할 때 beta 의 비율로 줄어드는 것으로 되어 있는데, 그렇지 않아도 된다. 그렇다면 (2) 번 식은 조금 달라지겠지. 문제 상황에 잘 맞게 수식을 세우면 된다. [본문으로]
3. matlab 에선 함수를 작성할 시, 함수 이름.m 파일로 저장해야 그 함수를 외부에서 호출해 사용할 수 있다. [본문으로]