operator01.jpg

안녕하세요? LGS입니다.
정말 오래간 만에 강좌를 진행하네요 ^^
그동안의 성원과 격려에 진심으로 감사드립니다.

오늘은 ODE(Ordinary Differential Equation 상미분방정식 常微分方程式)에 대해서 강의를 진행하려고 합니다.
공학이나 물리, 수학, 심지여 경제학에서도 이 (상)미분방정식이라는 것이 정말 정말 정말 중요하고 많이 쓰입니니다. 
특히 공대생의 경우 미분방정식을 모르면 대학생활이 꼬이기 시작하죠 -_-;;
오늘 ODE 잘 배우셔서 문과 애들에게 '니네들이 미방을 알어?' 하고 우쭐대 줍시다 ^^
자 그러면 시작하겠습니다~

uU4xBh0h4NBxNSiME0TNi_jZ7uMCvgI4KGQXKvI-qaH0bbLy62UDjDb335n2SkxqojoJ_cjOshH_qdyyZuKodFUcDUe28K2kb7mjFs5bytlCqCV_MWY


그럼 우선 ODE(상미분방정식) 가 무엇인지 알아야겠죠? 
ODE(Ordinary Differential Equation)라는 말을 구성하고 있는 단어들을 하나씩 뜯어 보죠.
Equation(방정식)이란 말은 다들 아실테고(해를 찾는거죠) 우리가 방정식이라는 말에서 주목해야할 부분은 
독립변수와 종속변수 입니다. 함수로 보자면 독립변수는 input 이구요 종속변수는 output 입니다.

Differential Equation 은 방정식은 방정식인데 그 방정식 안에 미분항이 들어 있는 겁니다.
우리가 익히 잘 잘고 있는 Newton의 방정식 F=ma 이 식도 미분 방정식이겠죠? ^^
a(가속도)가 시간에대한 2차 미분이니까요~ 실제로 물리계를 수식으로 표현해보면 
대부분 미분방정식의 형태를 띄게됩니다.

통상적으로 시간에 대한 미분방정식은(즉 시간이 흘러감에 따라 값이 어떻게 변하는지 볼 때는) 
IVP(Initial Value Problem  초기치문제)에 속하구요
공간에 대한 미분방정식은 BVP(Boundary Value Problem 경계치문제)에 속합니다.
또 하나 시간에 대한 미분에는 dot을 쓰구요.  ý= dy/dt   
공간에 대한 미분에는 prime을 쓰는 것도 알아 두세요 y' = dy/dx

미분방정식의 해는 특정한 값(숫자)이 아닙니다.
(ODE의 경우)미분방정식을 만족하는 y = f(x)의 f(x)를 찾아 내는 것이 목적입니다.

Ordinary Differetial Equation 에서 Ordinary 라는 말은 독립변수가 딸랑 1개 라는 말입니다.
즉 ODE는 풀어서 말하자면 독립변수가 1개이고 미분항이 포함된  방정식입니다. 
줄여서 '상미방' 입죠~
반면 독립변수가 2개 이상이면 PDE(Partial Differential Equation) 입니다.
ODE는 주로 Rigid Body Dynamics(강체동력학)에서 사용되구요
PDE는 주로 Flexible Structure Dynamics 등에 사용 됩니다.
(MATLAB에서 PDE Toolbox를 제공하고 있으니 참고하세요)

MATLAB에서 ODE를 푸는 방법은 여러가지가 있습니다.  
Numerical하게(수치적으로 값만 구함) 푸는 방법, Symbolic하게(문자로 정확한 해를 구함) 푸는 방법, 
Simulink로 푸는 방법 ...  사실 Simulink는 Graphical ODE Solver입니다.
오늘은 ode45 함수를 사용한 수치적(Numerical) 풀이법을 설명하겠습니다.
정말 간단하니까 쉽게 따라 오실 수 있을 겁니다. ^^
가장 간단한 1st Order ODE로 시작하죠~

예를 들어 여러분이 공중에 정지해 있는 헬리콥터에서 낙하산을 매고 뛰어 내린다고 합시다.
뛰어 내리기 전에는 속도가 0이고 뛰어내린 후 부터 떨어지는 속도가 증가합니다.
떨어지는 속도가 점점 빨라질 수록 낙하산이 받는 공기저항도 점점 커지게 됩니다.
그럼 땅에 떨어질 때의 속도(Terminal Velocity)는 어떻게 될까요?
이걸 ODE로 풀어 보겠습니다요~ 

우선 위의 상황을 수학식(미분방정식)으로 표현(modeling)해야겠죠?
isJ9jZtor-O5N7OyVTqi7E6BBv-qYpUeOY_-gbrcpwv5i53zQ1ySP9XF7DRz9EXw24Nk0dfi5jKSDFqkRmyvfFpfPSXSGvgRB-YXbW0qj_ns2JGqoq0
<이미지 출처: www.iop.org>

유체역학에 따르면 공기저항(Drag force)를 C*v² 으로 표현 할 수 있습니다. 
(C는 항력계수, v는 속도)
ma = F 에서 F는 -mg + C*v² 이고 ma는 m(dv/dt)로 쓸수 있겠죠?
다시 정리해서 써보면  m(dv/dt) = -mg + C*v² 입니다.
MATLAB에서 미분방정식의 표준 표현식은 ý = f(t,y) 입니다.
양변을 m을 나누면 표준형식으로 바꿀 수 있겠네요 
M-File 하나 만들어 줍시다. paradrop.m으로 저장하시구요

function vdot = paradrop(t,v)
    g = 9.18; %중력가속도 m/s²
    m = 80; %몸무게 kg
    C = 0.4; %항력계수 임의로 0.4를 사용했습니다.
    vdot = C/m*v^2-g; %standard form


그리고 아래 명령어를 command window에 붙여 넣어 보세요

[t v] = ode45(@paradrop,[0 30],0);
plot(t,v);
xlabel('time(s)');
ylabel('velocity(m/s)');

CMyhu8kArXb_S3LWaAKywZnuHi_0lSN-tBizAO-T-xSPPYFHkTHb33Q3WTXX12ufddqlrfu7ik9Vz06ry1Tqmn29E_CNe_pdo5-9PuWK90eFswqpIQ8


[t v] = ode45(@paradrop,[0 30],0);

ode45 함수는 3개의 입력변수를 입력 받습니다.
처음의 @paradrop은 미분방정식을 function handle로 입력하겠다는 말이고
(function handle은 'MATLAB 때려잡기 Data Type 2부'를 참조하세요)
두번째의 [0 30]은 0초 부터 30초까지의 해를 구하겠다는 말입니다.
마지막의 0은 초기값(여기에서는 초기속도)입니다. 왜 IVP 라고 하는지 아시겠죠?
출력변수는 2개입니다. 처음의 t는 시간 vector 이고요 두번째 출력변수 v는 각각의 t에 대한 y값(위의 경우에는 속도) vector 입니다.
위의 그래프를 보니 낙하산의 종말속도는 -43m/s 쯤 되네요

위의 미분방정식은 좀더 정확하게 표현다면 
1st Order Linear Time-Invariant Oridinary Differential Equation (LTI ODE)입니다. 
이제 친구에게 가서 나는 '1차선형시불변상미분방정식'을 MATLAB으로 풀줄 안다고 자랑합시다. ㅋ

오늘을 여기까지 입니다.
다음 시간에는 Order reduction을 통한 2nd Order ODE를 다루겠습니다.
다음 시간에 뵙겠습니다. 좋은 하루 되세요 ^^
profile