📂Julia

ジュリアで微分方程式の解法を早期終了する方法

概要

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)