MATLAB ordinary differential equation , ode45

MATLAB 에서 미분 방정식에 대한 해를 구하는 여러 방법 중 수치적인 해를 구할 때 주로 ode23(), ode45() 와 같은 함수를 이용한다.

 

ode 는 ordinary differential equation 의 약자이다.

 

뒤에 붙은 23이나 45 같은 숫자는 Runge-Kutta formulas 의 차수를 말한다.

 

ode23 은 2차 3차 Runge-Kutta 방식을 이용하는 것이고,

 

ode45 은 4차 5차 Runge-Kutta 방식을 이용하는 것이다.

---

  

문제의 유형에 따라 몇 차의 Runge-Kutta 방식을 이용해야 할지 결정 해야 한다.

 

문제의 유형이 stiff, nonstiff 인지 그리고 해의 정확도 가 High 인지 Low 인지에 따라 solver 의 선택을 달리 해야 한다.

 

언제 어떤 solver 를 선택해야 하는지는 다음 페이지를 보면 자세히 나와 있다.

 

http://www.mathworks.com/help/techdoc/ref/ode23.html

 

보통의 경우엔 ode45를 사용하곤 한다.

 

ode45() 함수의 일반적인 사용법은 다음과 같다.

 

[T,Y] = solver(odefun,tspan,y0)

 

odefun : 풀고자 하는 ordinary differential equation 함수의 handle 을 입력한다.

 

tspan : 시간 구간을 입력한다. 일반적으로는 [t0(초기값), tf(끝 값)] 을 입력하며, 때로는 [t0, t1, t2,…, tf] 와 같은 형태로 중간의 값들을 지정해 주기도 한다.

 

y0 : 해의 초기 상태를 의미 한다.

 

T : 시간이 나온다.

 

Y : 해가 출력 된다. 

---

오늘의 포스팅에서는 ode45() 를 이용하여 물통의 수량에 대한 문제를 풀어보려 한다.

 

 

위 문제에서 각 물통의 물의 높이 h1, h2, h3 는 다음과 같은 미분 방정식으로 표현된다.

 

위 식에서 A=3, q_i = 15, C=5 이다.

 

q_d 값은 0 ~ 4 초 동안 0 이다가 4초가 되는 순간부터 9 초까지 q_d = 20 이고 10초 이후에는 다시 q_d = 0 이 된다고 하자.

 

h₁ 의 초기값은5 , h₂ 의 초기값은5, h₁ 의 초기값은0 이다.

 

위와 같은 경우의 정상 상태에서의 각 물통의 물의 높이를 구해 보자.

 

일단 위 식에 따라 미분 방정식 함수는 다음과 같다.

function dhdt = mass_balance(t,h)

 

global C A qd qi % 각 변수를 global 로 설정

 

dhdt = zeros(3,1); % 출력 버퍼 설정

 

dhdt(1) = qi/A - C*sqrt(h(1))/A; % 첫번째 물통

 

dhdt(2) = qd/A + C*sqrt(h(1))/A - C*sqrt(h(2))/A; % 두번째 물통

 

dhdt(3) = C*sqrt(h(2))/A - C*sqrt(h(3))/A; % 세번째 물통

 

위 문제에 대한 풀이 코드는 아래와 같다.

 

 

100 초 까지 시뮬레이션 결과 각 물통의 물의 높이는 아래와 같이 됨을 확인 할 수 있다.

 

steady state h value: 9 , 9 , 9