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

 

odeordinary 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, qi = 15, C=5 이다.

 

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

 

h1 의 초기값은5 , h2 의 초기값은5, h1 의 초기값은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

 


  1. 장누리 2016.11.12 00:08

    안녕하세요. 좋은 글 잘봤습니다. 저는 생물공정에 대해서 공부하는 대학원생 입니다. 현재 제가 운영하는 생물 반응 시스템을 수치적으로 시뮬레이션하기위해서 전산 프로그램을 공부중인데요..문제는 이러한 전산 프로그램에 대한 기초가 전혀 없습니다.. 학교에서 제공하는 매스매티카를 설치해서 조금씩 공부해본 결과 간단한 ( 제가 말하는 간단한이란 의미는 dxdt 가 딱 한 종류만 있는) 미분방정식을 풀고 그래프화하는 방법까지는 알아냈습니다. 하지만 제가 시뮬레이션해야 할 시스템은 dxdt,,dydt,dzdt... 이런식으로 파라미터들이 거의 6,7 개 정도가 되며 각각의 변화량이 다른 파라미터들의 함수가 되는 형태입니다. 또 한 종류의 파라미터는 시간에 대한 함수가 아니라 다른 파라미터 값에 대한 함수인데... 어디서부터 손을 대야 할 지 잘 모르겠습니다.. 사실 지금 위의 시스템을 엑셀 파일로 구현을 하는데에는 성공했습니다만.. 좀 더 편하게 매트랩이나, 매스매티카 같은 툴을 이용해서 쉽게 풀고 GUI를 만들고 싶습니다.. 혹시 도움을 얻을 수 있을까요? ...

    • 남성 2016.11.12 00:33 신고

      구현하면서 MATLAB 사용과 관련하여 궁금한 부분이 있다면 답변은드릴 수 있지만 구현을 대신 해 드릴수는 없습니다.

  2. 이태빈 2018.05.13 18:06

    mass_balance 함수에서 입력 인자로 t,h를 쓰는데 함수 구현부에서 t가 왜 안쓰이는 지 궁금합니다.

    • 남성 2018.05.13 20:49 신고

      원래 mass_balance 함수에서는 위 수식을 보시면 변수 t 값을 사용하지 않았습니다.그런데 ode45 를 사용하기 위해서는 첫번재인자가 함수 핸들인경우 tspan 을 넣어줘야 하기 때문에 mass_balance 함수에 첫번째 인자로 t 값을 준 것입니다.

  3. 김진희 2018.11.08 11:20

    궁금한 것이 있는데요, 보니까 ode45의 미분방정식을 푸는데 t를 1~100의 벡터로 놓고,

    1과 2사이를 0.2씩 5구간으로 나눠서 for문으로 100번 반복하여 풀던데요,


    아예 처음부터 ode45함수에 500구간을 모두 넣지 않고, 5구간 씩 for문으로 100번 반복하게 하신 이유는 무엇인가요?


    저는 다음과 같이 풀었습니다.


    global C A qi

    % 출처: http://iamaman.tistory.com/71

    C = 5;
    A = 3;
    qi = 15;

    t = 0 : 0.2 : 100;

    h0 = [5 5 0];


    [t,h] = ode45 (@mass_balance, t, h0);

    figure, plot(t, h), grid on, legend('h1','h2','h3')

    fprintf('steady state h value: %g , %g , %g \n',h(end,:))


    여기서 qd는 qd.m파일을 만들어서 mass_balance.m파일에서 불러다 쓰는 함수로 만들어 놓았습니다.


    여기에 답변해주신다면 고맙겠습니다.

    doskid2@naver.com

  4. 김진희 2018.11.08 11:50

    문제를 다시 읽어보니까 for문을 쓰신 이유는 h값이 시간에 따라 변하기 때문이신것 같아요.

    이렇게 접근해보니까 for문을 100번 돌리는 것보단 아예 t를 0.2씩 증가시켜가며 500번 돌리는 것이 더 정확할 것 같아요.

  5. 김진희 2018.11.08 17:03

    변형시켜서 해보니까 님의 해법이 모범예시였군요.
    감사합니다.

+ Recent posts