# Maximum of matrix vector product for "large" matrices

I am using MTK and noticed that when I construct an algebraic variable defined as the max of the matrix-vector product of the state vector, the results are accurate for small matrices, but go inaccurate rapidly as the matrix size increases.
Specifically, I get accurate results up 100x10. For 500x10, it’s always inaccurate.
Copy-paste example (you can simply uncomment the definition of C):

``````using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)

function F(;name, α0 = fill(0.0, 10), C = fill(1.0, 10, 10))

sts = @variables T(t)=0.0 α(t)[1:length(α0)]=α0

eqs = [ scalarize(D.(α) .~ 0.0)
T ~ maximum(C*α)
]

ODESystem(eqs, t, [α..., T], []; name=name)
end

α0 = rand(10)

#C = rand(100, 10)
C = rand(500, 10)

@named mySys = F(α0=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
sol = solve(prob)

println(maximum(C*α0))
println(sol[mySys.T][end])
``````

Any clues how to get around that using MTK? If i build the same system directly in DifferentialEquations with an algebraic state defined the same manner, I don’t run into the same problem.

@YingboMa @shashi could this be fixed by using the MuladdMacro expression changing on the generated code?

I see very similar values for the last two lines

``````julia> println(maximum(C*α0))
4.421761822317552

julia> println(sol[mySys.T][end])
4.421761822317553
``````

I started a bigger system (1000x100) running but it hasn’t finished yet…

``````(MTK) pkg> st
Project MTK v0.1.0
Status `~/Desktop/Julia/MTK/Project.toml`
[0c46a032] DifferentialEquations v7.6.0
[961ee093] ModelingToolkit v8.30.0
[91a5bcdd] Plots v1.35.5
[731186ca] RecursiveArrayTools v2.32.1
[0c5d862f] Symbolics v4.13.0
``````

Edit
1000x100 finished!

``````julia> println(maximum(C*α0))
33.386913594597175

julia> println(sol[mySys.T][end])
33.38691359459718
``````

Ah, good find! I checked and made sure the package versions were the same, but it didn’t make a difference in Julia 1.7.2. When I tried in Julia 1.8, I got nearly the same result (not equal but good enough)! Thank you for pointing that out!

Ok, now for part II now that the max issue is resolved by switching to Julia 1.8.2 :

My original code looks a lot more like:

``````using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)

function myComponent(;name, Nm = 10, α0=zeros(Nm))

sts = @variables α(t)[1:Nm]=α0
eqs = scalarize(D.(α) .~ 0.0)

ODESystem(eqs, t, [α...], []; name)
end

function F(;name, α0s = [fill(0.0, 10) for i in 1:20], C = fill(1.0, 10, 200))

components = [myComponent(name = Symbol(string(i)), α0 = α0s[i] ) for i in 1:length(α0s)]

α = reduce(vcat, [component.α for component in components])

@variables T(t)=0.0

eqs = [  T ~ maximum(C*α)
]

compose(ODESystem(eqs, t, [T], [],; name), components...)

end

α0 = fill(rand(10), 20)

#C = rand(100, 10)
C = rand(500, 200)

@named mySys = F(α0s=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
@time solve(prob)
``````

When I run that on Julia 1.7.2. I get:

``````5.666656 seconds (10.51 M allocations: 634.450 MiB, 1.81% gc time, 99.94% compilation time)
retcode: Success
``````

When i try to run that on Julia 1.8.2, it takes > 5min (I’m still waiting for the answer and will edit).

Now I’m sure I’m doing something dumb with this concatenation (if this wasn’t symbolic, i’d at least pre-allocate). Any tips for a newbie at symbolic math/MTK?

Thanks for the help!

It’s been 15min now, so I give up.