Detecting and solving small dense linear subsystems symbolically after tearing

[YingboMa]

Consider the RC circuit model in https://github.com/SciML/ModelingToolkit.jl/blob/805811b9ef2da00c3fa5f7460ebf61ecfc1bd9cb/docs/src/tutorials/acausal_components.md#copy-paste-example

We have

julia> sys = structural_simplify(rc_model);

julia> equations(sys)
2-element Vector{Equation}:
 0 ~ capacitor₊v(t) + resistor₊R*capacitor₊p₊i(t) - source₊V
 Differential(t)(capacitor₊v(t)) ~ capacitor₊p₊i(t)*(capacitor₊C^-1)

julia> states(sys)
2-element Vector{Any}:
 capacitor₊v(t)
 capacitor₊p₊i(t)

Note that the algebraic equation is linear w.r.t. the algebraic variable capacitor₊p₊i(t). We could add a post tearing pass that solves linear systems symbolically.

Another example is

using ModelingToolkit, Random

N = 5
@variables xs[1:N]
A = rand(MersenneTwister(123), [-1, 1], N, N)
eqs = xs .~ A * xs
sys′ = NonlinearSystem(eqs, xs, [])
sys = tearing(sys′)

We can see that after tearing, we are left with

julia> equations(sys)
4-element Vector{Equation}:
 0 ~ xs₃ - xs₁ - xs₄ - xs₅ - (2(xs₃ + xs₄ - xs₅))
 0 ~ 2xs₄ - xs₁ - xs₃ - (2xs₅)
 0 ~ -xs₁ - xs₃ - (2xs₄) - xs₅ - (xs₃ + xs₄ - xs₅)
 0 ~ xs₁ + 2xs₃ + 2xs₄ - xs₅

julia> states(sys)
4-element Vector{Any}:
 xs₁
 xs₃
 xs₄
 xs₅

We can solve the resulting system via a linear solve.

[ChrisRackauckas]

JuliaSimCompiler