MATLAB에서 4차 Runge Kutta를 이용하여 1차 혹은 2차 미분방정식을 푸는 예제

물론 MATLAB에는 미분방정식을 푸는 멋진 함수도 이미 준비되어 있고, 또 Simulink라고 하는 훌륭한 도구도 있기 때문에 오늘 제가 이야기할 Runge Kutta의 MATLAB 예제 코드는 MATLAB만 놓고 보면 큰 의미가 없습니다. 그러나 예약발생을 할 이 글을 작성하고 있는 현재 시간이 새벽 1시 반인데, 아직도 아가 미바뤼가 잠들지 않은 이유도 있고, 또 C/C++이나 Python과 같은 다른 언어에서 구현하는 것을 고민하는 분들이라면 m-file로 보여드릴 이 예제가 혹시 도움이 되지 않을까 생각합니다. 그래서 오늘은 간편한 1차 미분방정식에서 흔히 RK4라고 하는 4차 Runge Kutta 방법을 소개하고 예제 코드를 보이고, 또 2차 미분방정식에서도 그렇게 하도록 해보겠습니다. 뭐~ 사실 인터넷만 좀 뒤지면 나오는 것입니다만~~~^^

-. 1차 미분 방정식을 Runge Kutta로 풀기

먼저 위와 같은 1차 미분 방정식에서 일차항만 남기고 이항해서 정리하면 

이 되겠죠. 그리고 나서

미소 변위 dx를 구할 수 있게 되는 거지요. 그러면 시간상 다음 값을

이렇게 표현해 볼 수 있게 됩니다.^^. 이제 예제하나 들어 보죠^^

위 일차 방정식을 풀어서 그래프로 나타내고자 합니다. 물론 MATLAB이 제공하는 ODE함수를 사용하는 것이 아니라 위의 전개된 수식을 이용해서 직접 코드를 짜는거지요. 위 예제는 실제 손으로 직접 풀면

위의 식입니다. 그래서 Runge Kutta로 푼것과 위의 실제 값을 비교도 해봐야지요. 아 x와 y가 혼용되었네요. 위의 식도 2y가 아니라 2x이고 그 아래의 y(t)도 x(t)입니다.ㅠㅠ. 먼저 코드를 한 번 보세요...

h = 0.01;
t = 0:h:4;

xReal = 2.5*exp(-2*t)-1.5*exp(-4*t);

% y' + 2y = 3exp(-4t)
func = @(tVal, yVal) 3*exp(-4*tVal)-2*yVal;

x = zeros(size(t));

x(1) = 1; %initial

for n = 1:(length(t)-1)
    k1 = func( t(n), x(n) );
    k2 = func( t(n) + 0.5*h, x(n) + 0.5*h*k1);
    k3 = func( t(n) + 0.5*h, x(n) + 0.5*h*k2);   
    k4 = func( t(n) + h, x(n) + h*k3);   
    
    x(n+1) = x(n) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end

figure
subplot(2,1,1)
plot(t,x)
grid on
title('RK4')
subplot(2,1,2)
plot(t, xReal)
grid on
xlabel('time (sec)')
title('Real')

위에서 for 문안에 있는 k1~k4와 y(k+1)이 바로 Runge-Kutta를 푼 것입니다. 주석 일부분에서 y라는 문자가 보이지만 그냥 x라고 생각해주세요.ㅠㅠ. 그냥 그러려니 해주세요.ㅠㅠ. 아무튼 이코드는 뭐 원체 위 수식과 비교해서 설명할 것이 없네요. 

결과 입니다. 위의 그래프가 m-file로 짠 Runge-Kutta이구요 아랫 그림이 실제 참값입니다. 잘 맞아 떨어지요?^^ Runge-Kutta로 1차 미방을 풀었네요^^.

-. 2차 미분 방정식을 Runge Kutta로 풀기

넵.. 전형적인 2차 미분방정식의 형태이구요.

여기서 상태 x의 미분된 상태를 v로 추가합니다. 

그러면, 위와 같이 표현되구요

그리고 이렇게 되어서 Runge Kutta를 두 번 푸는 것 처럼 볼 수 있지요...^^

이런 느낌은 위의 1차 미방과 동일 하구요^

 

kx1부터 잡아가는데, v쪽 kv1~kv3도 사용된다는 것에 쥐의하구요 

 

 

그리고, kv1부터 또 보면 

 

 

로 표현됩니다. 그리고 이들을 취합하는 것은 

 

 

상태가 x와 그 미분항인 v가 있었으니 두개가 되구요 

 

 

이렇게 1차 처럼 합치면 되지요...

역시나 예제하나 보겠습니다. 그냥 안정한 2차 시스템하나 잡죠. 

 

 

입니다. 

 

 

비교를 위해 simulink에서 위 예제를 꾸몄어요. 그리고, 초기치는 x에만 1을 주었습니다. 이제 MATLAB m-file로 작성한 2차 미방을 푸는 Runge-Kutta 방법의 코드를 보죠 

 

% x" + 2x' + x = 0, x = 1, v = 0

h = 0.1;
t = 0:h:5;

func = @(tVal, xVal, vVal) -2*vVal-xVal;

x = zeros(size(t));
v = zeros(size(t));

x(1) = 1;
v(1) = 0;

for n = 1:(length(t)-1)
    kx1 = v(n);
    kv1 = func( t(n), x(n), v(n) );
    
    kx2 = v(n) + h*kv1/2;
    kv2 = func( t(n) + h/2, x(n) + h*kx1/2, v(n) + h*kv1/2 );
    
    kx3 = v(n) + h*kv2/2;
    kv3 = func( t(n) + h/2, x(n) + h*kx2/2, v(n) + h*kv2/2 );
    
    kx4 = v(n) + h*kv3;
    kv4 = func( t(n) + h, x(n) + h*kx3, v(n) + h*kv3 );
    
    dx = h*(kx1 + 2*kx2 + 2*kx3 + kx4)/6;
    dv = h*(kv1 + 2*kv2 + 2*kv3 + kv4)/6;
    
    x(n+1) = x(n) + dx;
    v(n+1) = v(n) + dv;
end

figure
subplot(2,1,1);
plot(t, x)
grid on
subplot(2,1,2);
plot(t, v)
grid on

 

역시 코드 자체에 다를 것은 없습니다. 너무 심플해서 미분방정식이라는 놈을 이렇게 쉽게 풀어도 되나~ 싶을 정도죠^^ 

 

Simulink와 비교한 겁니다. 뭐 당연하지만 같지요^^ 이제 이 Runge Kutta를 이용해서 C/C++이나 다른 언어에서 좀 재미있는 것을 해볼 수 있지 않을까하는 쓸데 없는 생각을 합니다.^^