📂Julia

ジュリアでの暗示的微分方程式と微分代数方程式の解法

コード

次のコードはヘイスティングス・パウエルシステムを暗示的SINDyで復元して暗示的微分方程式を作り、これをジュリアの微分方程式ソルバーパッケージで解いてみる例である。

using Plots
import DifferentialEquations as DE
import Sundials
import DiffEqBase
function f2(out, du, u, p, t)
    out[1] = 2u[1] + 0.15u[1]^2 - 0.005u[1]^3 - 0.1u[1]*u[2] - 0.1u[1]*du[1] - du[1]
    out[2] = -u[2] - 0.1u[2]^2 + 0.1u[1]*u[2] - 0.15u[2]*u[3] - 0.015u[1]*u[2]*u[3] + 0.01u[1]*u[2]^2 - 0.1u[1]*du[2] - 0.1u[2]*du[2] - 0.01u[1]*u[2]*du[2] - du[2]
    out[3] = -0.7u[3] + 0.065u[2]*u[3] - 0.05u[2]*du[3] - du[3]
end
prob = DE.DAEProblem(f2, zeros(3), rand(3), (0, 200), differential_vars = [true, true, true])
sol = DE.solve(prob, Sundials.IDA(); initializealg = DiffEqBase.BrownFullBasicInit())
plot(sol)
plot(eachrow(stack(sol.u))...)

• Sundialsは暗示的微分方程式を解くためのソルバーであるIDAを提供する¹。
• DiffEqBase.BrownFullBasicInit()は初期値の状態と微分係数を合わせる機能をする。

見てのとおりf2をよく調べると、通常は常微分方程式の左辺にのみあるべき項であるduが右辺にも登場することがわかる。必ずしも暗示的微分方程式でなくとも、outという変数が$0$になるように微分方程式を書けばよい。これはたとえば次のような表現も可能であるという意味になる。

$$ \begin{align*} f(\mathbf{x}, \dot{\mathbf{x}}) =& 0 \\ x_{1} + x_{2} = & 1 \end{align*} $$