eulode 사용법 : 매트랩 수치적분 Euler ODE solver

수치해석을 배우게 되었어요. 너무 어려워요. 하쥐만 저는 이과의 자랑(전혀 아님), 한번 해보겠어요 

 

Euler법

 상미분방정식의 초기값 문제를 푸는 방법 중 하나로 현재 값 하나로 다음 값을 예측하는 것이다. 즉, 현재의 위치에서의 접선(기울기)을 따라 아주 조금 이동하며 다음 위치를 예측하는 방법이다. 이 Euler법의 원리는 곡선을 아주 작은 직선들의 연결로 가정하는 것이다. 

 

 언제나 그렇듯 세상의 많은 미분방정식은 수학적으로 딱 떨어지는 해가 없다. 로켓의 궤도 또한 공식으로 풀기 너무 복잡하다. 그러나 공식이 없어도 현재의 변화율(기울기)만 알면, 오일러 방법으로 값을 하나씩 계산해 나가며 해를 찾아낼 수 있다. 물론 오일러 방법은 굉장히 간단하기에 실제로 많이 사용하지는 않는다. 그런데 왜 배우는가에 대해서 고민해보았다. 

 

  오일러 방법은 가장 단순하지만 실무에서 많이 쓰이는 정확하고 복잡한 방법 ( ex. Runge-Kutte 룬게-쿠타)의 기본 원리를 담고 있다. 또한 구현이 매우 쉽고 빨라 아주 정밀한 값이 필요 없는 초기 시뮬레이션에서는 계산 속도가 빨라 유용하다. 즉 어떤 복잡한 미분방정식이라도 일단 답을 구해낼 수 있는 가장 쉬운 해결책이기에 사용한다. 

 

출처 : 응용수치해석

 

새로운 값 = 이전 값 + 기울기 x 간격크기

- y_n+1 : 예측할 다음 값
- y_n : 현재 값 

- h : 단계 크기 (Step size, 옆으로 이동할 간격)

- f(x_n, y_n) : 현재 위치에서의 미분값(기울기) 

 

 

그럼 이 Euler법을 matlab에서 어떻게 쓸 수 있는지에 대해 알아보자 ! 

 

MATLAB eulode정의 

function [t,y] = eulode(dydt, tspan, y0, h)
% [t,y] = eulode(dydt, tspan, y0, h, p1,p2,...)
%       uses EULER's method to Integrate at ODE
%       (uses the slope at the beginning of the stepsize to grpah the function)
% Input:
%       dydt = name of hte M-file that evaluates the ODE
%       tspan = [ti, tf] where ti and tf = initial and final values of independent values
%       y0 = initial value of dependent variable
%       h = step size
%       p1,p2 = additional parameter used by dydt
% Output:
%       t = vector of independent variable
%       y = vector of solution for dependent variable

if nargin<4, error('at least 4 input arguments required'), end
ti = tspan(1); tf = tspan(2);
if ~ (tf>ti), error('upper limit must be greater than lower limit'), end

t = (ti:h:tf)';
n = length(t);

if t(n)<tf
        t(n+1) = tf;
    n = n+1;
    t(n) = tf;
end

y = y0*ones(n,1);

for i = 1:n-1
        y(i+1) = y(i) + dydt(t(i),y(i))*(t(i+1)-t(i));
end

plot(t,y);

 

 

MATLAB eulode 사용하기 

예제

dydt = @(t,y) 4*exp(0.8*t)-0.5*y;

[t,y] = eulode(dydt,[0,4], 0, 1)

 

dy/dt를 정의하고 tspan, y0, h의 값을 각각 정의해주면 값이 잘 나오는 것을 확인할 수 있답니다. 

시간에 따른 y의 변화를 Euler법을 사용하여 계산한 결과 

t =

     0
     1
     2
     3
     4


y =

    2.0000
    5.0000
   11.4021
   25.5132
   56.8493

잘 나오는 것을 확인할 수 있다! 

해석해와 Euler로 예측한 다음 해

이렇게 수학적 정의만 코드로 잘 옮겨오고 eulode만 잘 사용하면 된다.