operator01.jpg

안녕하세요? LGS입니다.
오늘은 2nd order ODE를 풀어봅시다. 
앞의 강의를 잘 따라 오셨다면 이번 강의도 매우 쉬울 겁니다.

nA7DXDMOZL0TAcvskBSTrnFwPCtyt9gY5bBEfhz2JjF4Oc7imLn_RJiQfwcdUaHqqbnmq-CemHizfRKgfjAALNI6FTTCh17VUpG39DdzcHDJn7g_Kig

위의 표를 보시고 우리가 어디에 와있는지를 잘 파악해보세요~
보통 물리, 공학에서는 2차 ODE 까지만 쓰입니다. 
거리에 대한 시간 1차미분이 속도, 2차미분이 가속도죠~
물리적으로 3차 미분은 jerk(가가속도)라는 개념이 있지만 잘 쓰이지 않습니다.
오늘 우리가 가지고 놀 녀석은 mck System 이란 놈입니다. 
(Mass Spring Damper로 구성된 System)
mck system은 2nd Order ODE로 나타낼 수 있습니다.

sCV__m9JTHP6u4ffrY7VgUTvSx3z4cz5lQ1j5Wb5Kn-P5Ok7vbUkkAc4oA9iMEvfr7RZQrdWjCM6kEqWUJXp46PN8Iv2R_AJr2oEjgBwEWTN6GxXkj0


2차 이상의 ODE를 MATLAB에서 ode45 함수를 사용해서 Numerical하게 푸실려면 아래 3단계를 거쳐야 합니다.
(미분방정식을 세우는 방법은 Modern Control Theory with MATLAB의 Mathematical Modeling 강좌 시리즈를 참조하시기 바랍니다.)

1. n차 ODE를 n개의 1차 ODEs로 변환한다 (Order Reduction)
2. M-file ODE Function을 만들어 준다. (State Space)
3. ode45 함수에 적용한다. (ode45(@ODE함수,[시작시간, 끝시간],[초기값들])



지난 시간에 3번은 이미 했습니다.
그럼 첫번째 step인 Order Reduction 에 대해서 이야기하겠습니다.

차원축소(order reduction)은 n차의 미분방정식을 n개의 1차 미분방정식으로 만든는 것을 의미합니다. 

왜 order reduction을 해 주는 걸까요? 그 이유는 컴퓨터는 깡통이기 때문입니다. 

컴퓨터(MATLAB)은 정형화된 ODE에 대해서만 아주 잘 풀 수 있습니다. 

만약 우리가 정형화되지 않은(즉 2차 3차... 고차)의 ODE를  직접 사용한다면 

MATLAB에도 2차, 3차... 고차 ODE에 상응하는 solver 함수를 만들어 두어야 합니다. 

사용자가 100차 ODE를 풀고 싶은 것을 가정해서 MATLAB에 100차 ODE에 해당하는 solver 함수를 만들어 둘껍니까? -_-;; 

왜 Order Reduction 하는지 아시겠죠?

자 그럼 2nd order ODE를 2개의 1st order ODEs로 만들어 봅시다.
약간의 Trick을 쓰면 쉽게 order reduction을 할 수 있습니다.

b3uEbgUYbrdqKOYBW_HqRjv_oUA2M31SMzR4FdNMjP1NjTYeEmo0La_ZAAbNa2IaXMn5AHBF-3doiwZFUbikncmMBAnSvbhK1T7eANLet-Vf2I7ltqY

어떻습니까?(x1과 x2를 어떻게 잡았는지 잘 보세요)
Order Reduction을 통해 2차 ODE가 1차 ODE 두개로 바뀌었죠?
이제 우리가 풀고자 했던 ODE가 MATLAB이 알아먹을 수 있는 형태로 바뀌었습니다.
그럼 ODE Function M-file을 한번 만들어 봅시다.

fucntion dx = ode_fun(t,x)
   m=1;c=1;k=1; %상수는 임으로 정했습니다.
   u=1 %Unit Step 으로 힘을 가하는 경우를 보겠습니다.
   dx = zeros(2,1) %dx를 column vector로 고정시키기위해
   dx(1) = x(2);
   dx(2) = -k/m*x(1) - c/m*x(2) + 1/m*u;

간단하죠? 위의 그림에 있는 식을 그대로 썼습니다.
다른 것은 이해하시는데 문제가 없으실테고 저 u란 놈이 좀 아리송하겠네요~
상자가 있는 그림을 잘 보시면 바퀴가 달린 상자를 일정한 크기의 힘(1N)으로 계속 당기는 상황입니다.
마지막으로 ode45 함수에 적용해봅시다.

[t x] = ode45(@ode_fun,[0 15],[0 3]);
x_position = x(:,1);
x_velocity = x(:,2);
plot(t,x_position,t,x_velocity);
legend('Position','Velocity');
xlabel('Time(s)');
ylabel('Output');

alC3V8B-fws-91bW6Ah5UBu-RBdPYqLIjPnqvL8Y_1fnkXh9K7v7vpm3-bPXGnHeEKwL6LjdeH7ZCYR0b8f9P11DnC-UA8noeZmkIPy1Knv2WfTcMTI

[t x] = ode45(@ode_fun,[0 15],[0 3]) 에서 맨 마지막의 [0 3]은 무슨 의미일까요? 
처음 0은 상자의 초기 위치, 그다음 3은 상자의 초기 속도가 되겠습니다.

좀더 뽀대나게 ode_fun을 만들고 싶으시면 State Space Form을 이용해보세요
결과는 똑같습니다.

EA0CrcOtXkxbT1P_fwLWPsfSnp_cCPyUWPzYrW_BGvn339x9JMb0626NVcPixMI1WvMnJmUcKqJsH96ZkPq-aAuegolDsO2l6Fs8LqelXD_NOeMs24c


fucntion dx = ode_fun(t,x)
   m=1;c=1;k=1;u=1;
   A = [0  1 ; -k/m    -c/m];
   B = [0 ; 1/m];
   dx = A*x+B*u

고차 ODE에대해서도 위의 과정을 반복하시면 됩니다.
BVP(Boundary Value Problem)과 관련해서는 MATLAB에서  bvp4c 라는 함수를 제공하고 있으니 참고하세요


P.S. State Space에 대한 강좌 Update 되었습니다.

11강 - MCT 때려잡기 - State Space 1부 (Creating SS Model in MATLAB)


P.S. Boundary Value Problem (BVP)에 관한 강좌 Update 되었습니다.

http://www.matlabinuse.com/index.php?mid=Mastering_MATLAB&document_srl=9776


이상입니다. 다음시간에 뵙겠습니다

profile