Matlab ODE에 추가 파라미터 전달 및 이벤트 설정하기

일반적으로 매트랩(Matlab)에서 수치적분하고자 하는 ODE 함수의 구조는 다음과 같다.

 

dxdt = odefun(t,x)

 

여기서 odefun은 함수의 이름, t 는 시간, x 는 함수의 상태변수다. 함수에서 t 를 계산에 사용하지 않더라도 t 와 x 를 모두 입력값으로 받아야 한다.

그런데 만약 t 와 x 이외에 파라미터가 더 필요하다면 어떻게 해야 할까. 예를 들면, 다음과 같은 함수가 있다.

 

\[ \begin{align} \dot{\mathbf{x}}(t) &= \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{v}_x \\ \dot{v}_y \end{bmatrix} = \mathbf{f}(\mathbf{x}(t)) \tag{1} \\ \\ &= \begin{bmatrix} v_x \\ v_y \\ 2v_y+x- \frac{(1-\mu)(x+\mu)}{r_1^3 }- \frac{\mu (x+\mu-1)}{r_2^3} \\ -2v_x+y-\frac{(1-\mu)y}{r_1^3 }-\frac{\mu y}{ r_2^3 } \end{bmatrix} \end{align} \]

여기서

\[ \begin{align} & r_1= \sqrt{(x+\mu)^2+y^2 } \\ \\ & r_2= \sqrt{(x+\mu-1)^2+y^2 } \end{align} \]

 

이다. 이 함수에는 \(\mathbf{x}\) 이외에 추가적으로 \(\mu\) 가 더 필요하다.

일단 식 (1)의 수식을 매트랩 함수로 작성한다. 이 함수는 추가적인 파라미터 \(\mu\) 도 입력으로 받는다.

 

unction xvec_dot = orbit_ode(t, xvec, mu)
% orbit ODE for additional parameter test
% coded by st.watermelon

xvec_dot = zeros(4,1);

x = xvec(1);
y = xvec(2);
vx = xvec(3);
vy = xvec(4);

r1 = sqrt((x+mu)^2 + y^2);
r2 = sqrt((x+mu-1)^2 + y^2);

xvec_dot(1) = vx;
xvec_dot(2) = vy;
xvec_dot(3) = 2*vy + x - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
xvec_dot(4) = -2*vx + y - (1-mu)*y/r1^3 - mu*y/r2^3;

 

 

 

ode113 을 이용하여 방정식 (1)을 푼다. mu는 미리 정의된 값을 함수 orbit_ode 에 다음과 같은 형식으로 전달하면 된다.

 

% orbit ODE for additional parameter test
% need orbit_ode.m
% coded by st.watermelon

clear all

mu = 0.01215;
tspan = [0 5];
x0 = [0.8 0 0 0.8]';

OPTIONS = odeset('RelTol',1e-13,'AbsTol',1e-14);
[t,x] = ode113(@(t,x) orbit_ode(t,x,mu), tspan, x0, OPTIONS);

plot(x(:,1),x(:,2))
xlabel('x')
ylabel('y')

 

시뮬레이션 결과를 그림으로 그리면 다음과 같다.

 

 

적분구간은 보통 최종시간을 설정하여 결정한다. 예를 들면 앞의 예제의 경우 tspan = [0 5] 으로 적분구간을 0에서 5로 정의하였다. 하지만 간혹 적분구간을 미리 설정된 시간으로 정하기 보다는 특정 이벤트에 의해서 정해야 할 필요가 있다. 예를 들면 앞의 예제에서 궤적이 \(y=0\) 을 통과할 때 수치적분을 중지하고 그 때의 \(t, \ x, \ v_x, \ v_y\) 값을 알고자 하는 경우다. 또는 수치적분은 계속 진행하지만 궤적이 \(y=0\) 을 통과할 때의 \(t, \ x, \ v_x, \ v_y\) 값을 알고자 하는 경우가 있을 수 있다.

이와 같이 수치적분하는 동안 특정 이벤트가 발생하는 시점을 감지하기 위해서는 이벤트 함수를 사용하면 된다. 이벤트 함수는 사용자가 지정한 표현식을 사용하여 해당 표현식이 \(0\) 인 경우에 이벤트를 발생시켜 ode 솔버에 적분을 중단하도록 신호를 보낸다. 이벤트 함수의 형식은 다음과 같다.

 

[value,isterminal,direction] = myEventFcn(t,x)

 

여기서 value 가 0이 되면 이벤트가 발생한다. 이벤트가 발생했을 때 isterminal=1 이면 수치적분이 중단되고 isterminal=0 이면 수치적분은 계속된다. direction=1 이면 0을 찾을 때 이벤트 함수가 중가하는 곳의 0만 찾고 direction=-1 이면 감소하는 것의 0만 찾는다. direction=0 이 디폴트값으로서 이벤트 함수가 증가하거나 감소하는 곳의 0을 모두 찾는다. 이벤트 함수를 작동 시키려면 odeset 의 'Events' 속성에 이벤트 함수이름을 설정해야 한다.

 

options = odeset('RelTol',1e-13,'AbsTol',1e-14, 'Events', @myEventFcn);

 

그러면 ode 함수의 출력은 다음과 같이 3개의 항이 더 추가된다.

 

[t,x,te,xe,ie] = ode113(odefun,tspan,x0,options)

 

출력값에서 te 는 이벤트 발생 시간이고, xe 는 이벤트 발생 시 계산된 해이며, ie 는 이벤트의 인덱스이다.

앞의 예제에서 궤적이 \(y=0\) 일 때를 이벤트로 정의하고 그 때까지만 수치적분해 보자.

 

% orbit ODE for event test
% need orbit_ode.m
% coded by st.watermelon

clear all

mu = 0.01215;
tspan = [0 5];
x0 = [0.8 0 0 0.8]';

OPTIONS = odeset('RelTol',1e-13,'AbsTol',1e-14, 'Events', @events);
[t,x,te,xe,ie] = ode113(@(t,x) orbit_ode(t,x,mu), tspan, x0, OPTIONS);

plot(x(:,1),x(:,2),xe(:,1),xe(:,2),'ro')
xlabel('x')
ylabel('y')

te

function [value,isterminal,direction] = events(t,x)
%
value = x(2); % event y=0
isterminal = 1; % stop integration
direction = -1; % negative direction
end

 

수치적분은 \(y=0\) 에 해당하는 \(t=1.5637\) 에서 중단되었다.

 

 

isterminal = 0 으로 하면 이벤트 발생 여부와 관계없이 수치적분은 설정된 적분구간 tspan = [0 5] 동안 계속된다.

 

 

더 자세한 사항은 다음 사이트를 참고하면 된다.
https://kr.mathworks.com/help/matlab/math/ode-event-location.html