음이 아닌 ODE 해

여기서는 ODE의 해를 음이 아닌 값으로 제한하는 방법을 보여줍니다. 비음(Non-negativity) 제약 조건을 적용하는 것이 항상 간단하지만은 않지만, 방정식의 물리적 해석이나 해 특성으로 인해 이 제약 조건을 적용해야 하는 경우가 있습니다. 이 제약 조건은 이를 적용하지 않으면 적분에 실패하거나 해를 응용할 수 없게 되는 경우와 같이 필요할 때에만 해에 적용해야 합니다.

해의 특정 성분이 음수가 아니어야 할 경우 odeset을 사용하여 이러한 성분의 인덱스에 NonNegative 옵션을 설정합니다. ode23s와 ode15i, 그리고 질량 행렬을 갖는 문제에 적용되는 음함수 솔버(ode15s, ode23t, ode23tb)에는 이 옵션을 사용할 수 없습니다. 특히, 특이 질량 행렬을 반드시 갖는 DAE 문제에 비음 제약 조건을 적용할 수 없습니다.

예: 절댓값 함수

다음과 같은 초기값 문제를 가정하겠습니다.

$$y^{\prime }=-|y|,$$

초기 조건 $$y(0)=1$$을 사용하여 구간 $$[0,40]$$에 대해 이 방정식을 풉니다. 이 ODE의 해는 0으로 감쇠됩니다. 솔버에서 음의 해 값이 생성되면 솔버는 이 값을 통해 ODE의 해를 추적하기 시작하며, 계산된 해가 $$-\infty$$로 발산하므로 결국 계산이 실패합니다. NonNegative 옵션을 사용하면 이러한 적분이 실패하지 않도록 할 수 있습니다.

ode45를 사용하여 $$y(t)=e^{-t}$$의 해석적 해를 ODE의 해와 비교합니다. 이때 한 번은 추가 옵션을 지정하지 않고 비교하고 다른 한 번은 NonNegative 옵션을 설정한 상태에서 비교합니다.

ode = @(t,y) -abs(y);

% Standard solution with |ode45|
options1 = odeset('Refine',1);
[t0,y0] = ode45(ode,[0 40],1,options1);

% Solution with nonnegative constraint
options2 = odeset(options1,'NonNegative',1);
[t1,y1] = ode45(ode,[0 40],1,options2);

% Analytic solution
t = linspace(0,40,1000);
y = exp(-t);

비교를 위해 이 3개의 해를 플로팅합니다. 해가 $$-\infty$$로 떨어지지 않도록 하기 위해서는 비음(Non-negativity) 제약 조건을 반드시 적용해야 합니다.

plot(t,y,'b-',t0,y0,'ro',t1,y1,'k*');
legend('Exact solution','No constraints','Nonnegativity', ...
       'Location','SouthWest')

예: 무릎 문제(Knee Problem)

음이 아닌 해가 필요한 문제를 보여주는 또 다른 예는 무릎 문제이며, 예제 파일 kneeode에 코딩되어 있습니다. 방정식은 다음과 같습니다.

$$ϵy^{\prime }=(1-x)y-y^2,$$

초기 조건 $$y(0)=1$$을 사용하여 구간 $$[0,2]$$에 대해 이 방정식을 풉니다. 일반적으로 $$0<ϵ\ll 1$$을 충족하는 파라미터 $$ϵ$$이 사용되며, 이 문제에서는 $$ϵ=1\times 10^{-6}$$을 사용합니다. 이 ODE의 해는 $$x<1$$에서 $$y=1-x$$, $$x>1$$에서 $$y=0$$에 접근합니다. 하지만, 디폴트 허용오차를 사용하여 수치 해를 계산하면 해가 전체 적분 구간에서 $$y=1-x$$ 등경사선(Isocline)을 따릅니다. 비음 제약 조건을 적용하면 올바른 해가 결과로 생성됩니다.

비음 제약 조건을 적용한 상태와 그렇지 않은 상태 모두에서 무릎 문제(Knee Problem)를 풉니다.

epsilon = 1e-6;
y0 = 1;
xspan = [0 2];
odefcn = @(x,y,epsilon) ((1-x)*y - y^2)/epsilon;

% Solve without imposing constraints
[x1,y1] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0);

% Impose a nonnegativity constraint
options = odeset('NonNegative',1);
[x2,y2] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0, options);

비교를 위해 해를 플로팅합니다.

plot(x1,y1,'ro',x2,y2,'b*')
axis([0,2,-1,1])
title('The "knee problem"')
legend('No constraints','Non-negativity')
xlabel('x')
ylabel('y')

참고 문헌

[1] Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne, "Non-negative solutions of ODEs," Applied Mathematics and Computation Vol. 170, 2005, pp. 556-569.

참고 항목

odeset

관련 항목

• ODE 솔버 선택하기
• ODE 옵션 요약