ModelingToolkit - JuliaCon 2024

Trying out code from a talk by Chris Rackauckas

Acausal modeling tool, what does acausal mean?

Replicate some code:

using ModelingToolkit, StochasticDelayDiffEq, Plots, Latexify

@variables t, x(..) # What does the x(..) do?
@parameters a=-4.0 b=-2.0 c=10.0 α=-1.3 β=-1.2 γ=1.1
D = Differential(t)
@brownian  η
τ = 1.0
eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η]
@mtkbuild sys = System(eqs, t)
prob_mtk = SDDEProblem(sys, [x(t) => 1.0 + t], (0.0, 1.0); constant_lags = (τ,));
display(latexify(sys))
sol_mtk = solve(prob_mtk, RKMil())
plot(sol_mtk)

\[\begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= c + a x\left( t \right) + b x\left( -1 + t \right) \end{align}\]

(a) White noise

(b)

Figure 1

using Catalyst, DifferentialEquations, Plots, Latexify

rn = @reaction_network begin
    @parameters c₁ c₂ c₃ c₄ Ω
    (c₁/Ω^2), 2X + Y --> 3X
    (c₂), X --> Y
    (c₃*Ω, c₄), 0 <--> X
end
display(latexify(rn))
p = [:c₁ => 0.9, :c₂ => 2, :c₃ => 1, :c₄ => 1, :Ω => 100]
u₀ = [:X => 1, :Y => 1]
tspan = (0., 100.)

dprob = DiscreteProblem(rn, u₀, tspan, p)
jprob = JumpProblem(rn, dprob, Direct())
sol = solve(jprob, SSAStepper(), saveat=10.)
plot(sol)

\[\begin{align*} 2 \mathrm{X} + \mathrm{Y} &\xrightarrow{\frac{c_1}{\Omega^{2}}} 3 \mathrm{X} \\ \mathrm{X} &\xrightarrow{c_2} \mathrm{Y} \\ \varnothing &\xrightleftharpoons[c_4]{c_3 \Omega} \mathrm{X} \end{align*}\]

(a) catalyst

(b)

Figure 2

prob = ODEProblem(rn, u₀, tspan, p)
sol = solve(prob)
plot(sol)

Figure 3: catalyst-ode

using MomentClosure
central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false)
closed_raw_eqs = moment_closure(central_eqs, "normal")
u₀map = deterministic_IC([1, 1], closed_raw_eqs)
p = [0.9, 2, 1, 1, 100]
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob)
plot(sol)

Figure 4: moment-closure

using MethodOfLines
import ModelingToolkit: Interval, infimum, supremum

@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = Dxx(u(x,y)) + Dyy(u(x,y)) ~ sin(pi*x)*sin(pi*y)

bcs = [
    u(0,y) ~ 0.f0,
    u(1,y) ~ -sin(pi*1)*sin(pi*y),
    u(x,0) ~ 0.f0,
    u(x,1) ~ -sin(pi*x)*sin(pi*1)
]

domains = [
    x ∈ Interval(0.0, 1.0),
    y ∈ Interval(0.0, 1.0)
]
@named pde_system = PDESystem(eq, bcs, domains, [x,y], [u])

# dx = 0.05
# dy = 0.05
# discretization = MOLFiniteDifference([x=>dx, y=>dy], nothing, centered_order=2)

# prob = discretize(pde_system, discretization)
# sol = solve(prob)
# xs,ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
# u_sol = reshape(sol.u, (length(xs), length(ys)))
# plot(xs, ys, u_sol, linetype=:contourf, title="solutions")

PDESystem
Equations: Equation[Differential(y)(Differential(y)(u(x, y))) + Differential(x)(Differential(x)(u(x, y))) ~ sin(πx)*sin(πy)]
Boundary Conditions: Equation[u(0, y) ~ 0.0, u(1, y) ~ -1.2246467991473532e-16sin(πy), u(x, 0) ~ 0.0, u(x, 1) ~ -1.2246467991473532e-16sin(πx)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0.0 .. 1.0), Symbolics.VarDomainPairing(y, 0.0 .. 1.0)]
Dependent Variables: Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}[u⋆]
Independent Variables: Num[x, y]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

Figure 5: finite-diff

using NonlinearSolve
@variables u1 u2 u3 u4 u5

eqs = [
    0 ~ u1 - sin(u5),
    0 ~ u2 - cos(u1),
    0 ~ u3 - hypot(u1, u2),
    0 ~ u4 - hypot(u2, u3),
    0 ~ u5 - hypot(u4, u1)
]
vars = [u1, u2, u3, u4, u5]
@named sys = NonlinearSystem(eqs, vars, [])
display(latexify(sys))
sys = structural_simplify(sys)
display(latexify(sys))
display(observed(sys))
u0 = [u5 .=> 1.0]
prob = NonlinearProblem(sys, u0)
sol = solve(prob, NewtonRaphson())
[sol[u] for u in vars]

\[\begin{align} 0 &= u1 - \sin\left( u5 \right) \\ 0 &= u2 - \cos\left( u1 \right) \\ 0 &= u3 - \mathrm{hypot}\left( u1, u2 \right) \\ 0 &= u4 - \mathrm{hypot}\left( u2, u3 \right) \\ 0 &= u5 - \mathrm{hypot}\left( u4, u1 \right) \end{align}\]

\[\begin{align} 0 &= u5 - \mathrm{hypot}\left( u4, u1 \right) \end{align}\]

4-element Vector{Equation}:
 u1 ~ sin(u5)
 u2 ~ cos(u1)
 u3 ~ hypot(u1, u2)
 u4 ~ hypot(u2, u3)

5-element Vector{Float64}:
 0.9993449829954304
 0.540853367725015
 1.1363154317431527
 1.2584653852200776
 1.6069926947050053