📂줄리아

줄리아에서 암시적 미분방정식 및 미분대수방정식 푸는 법

코드

다음의 코드는 헤이스팅-파웰 시스템을 암시적 신디로 복원해서 암시적 미분방정식을 만들고, 이를 줄리아의 미분방정식 풀이 패키지로 풀어보는 예제다.

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*} $$