📂줄리아

줄리아에서 미분방정식 풀이 조기종료하는 법

개요

SIR 모델: $$ \begin{align*} {{d S} \over {d t}} =& - {{ \beta } \over { N }} I S \\ {{d I} \over {d t}} =& {{ \beta } \over { N }} S I - \mu I \\ {{d R} \over {d t}} =& \mu I \end{align*} $$

가령 위와 같은 감염병 확산 모델 등을 시뮬레이션한다고 할 때, 고정점은 $I = 0$ 에 있을 것이고 감염자가 충분히 $0$ 에 가까워졌다고 생각한다면 그 시점에서 시뮬레이션을 종료하고 싶을 수 있다.

코드

줄리아의 미분방정식 풀이 패키지에서는 이벤트 핸들러¹의 기능 중 terminate!를 사용해서 미분방정식의 풀이를 조기에 종료시킬 수 있다².

기존

prob = ODEProblem(SIR, rand(3), (0, 200))
sol = solve(prob)

임의의 초기조건에서 $t \in [0, 200]$ 구간까지 미분방정식을 풀면 그 해는 위와 같을 수 있다. 보다시피 $t > 50$ 이후로는 별달리 의미 있는 변화가 없고, 이러한 시뮬레이션이 수없이 반복되는 경우라면 굳이 $t = 200$ 까지 풀 필요가 없다.

조기종료

condition(u,t,integrator) = u[2] - 0.01
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
sol = solve(prob, callback = cb);

한편 $I < 0.01$ 이라는 조건 하에서 종료하라는 콜백이 있는 경우, condition 함수의 부호가 바뀌는 시점에서 affect!가 호출되고, 그 안에서 호출된 terminate!에 의해 미분방정식 풀이가 종료된다.

전체코드

using DifferentialEquations, Plots

function SIR(du, u, p, t)
    S, I, R = u
    du[1] = -0.5 * S * I
    du[2] = 0.5 * S * I - 0.1 * I
    du[3] = 0.1 * I
end
prob = ODEProblem(SIR, rand(3), (0, 200))

sol = solve(prob)
plot(sol)

condition(u,t,integrator) = u[2] - 0.01
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
sol = solve(prob, callback = cb);
plot(sol)