Can DifferentialEquations.jl solve a system of matrix DAEs expressed in a Mass-Matrix DAE Format with RecursiveArrayTools

Specific Domains
1

Hello, I have this following code :

using RecursiveArrayTools
using LinearAlgebra
using DifferentialEquations
#
matrix_3 = [-1/0.08 0 -0.2/0.08 ; 1/0.3 -1/0.3 0; 0 50/20 -1/20]
matrix_4 = [-1/0.08 0 -0.2/0.08 ; 1/0.3 -1/0.3 0; 0 50/20 -1/20]
matrix = VectorOfArray([matrix_3, matrix_4])
#
M =VectorOfArray([1.0*I(3), 1.0*I(3)])
#M =1.0*I(6) 
#
p = [0; 0; -0.05*(50/20)]
function powsys!(du, u, p, t)
    du[1] = matrix[1]*u[1] .+ p
    du[2] = matrix[2]*u[2] .+ p
end
u0= VectorOfArray([zeros(3,1), zeros(3,1)])
tspan = (0.0,20.0)
prob = ODEProblem(ODEFunction(powsys!, mass_matrix = M),u0,tspan, p)
sol = solve(prob,Rodas5(), reltol = 1e-8, abstol = 1e-8)

Defining M as a VectorOfArray will cause a MethodError because DiffrentialEquations.jl won’t be able to evaluate if the matrix is singular.

Does anyone know if what I want is possible ? In summary I want to know how the mass-matrix should be defined when the state variable is in a matrix format expressed with VectorOfArray.

PS : I am using VectorOfArray because later on my matrices will have different size.

2

what’s the stacktrace of the error?

3

This is the error I get :

MethodError: no method matching issingular(::VectorOfArray{Float64, 3, Vector{Diagonal{Float64, Vector{Float64}}}})

Closest candidates are:
  issingular(::SymTridiagonal)
   @ ArrayInterface C:\Users\lenovo\.julia\packages\ArrayInterface\ubN3d\src\ArrayInterface.jl:328
  issingular(::Matrix)
   @ ArrayInterface C:\Users\lenovo\.julia\packages\ArrayInterface\ubN3d\src\ArrayInterface.jl:324
  issingular(::Tridiagonal)
   @ ArrayInterface C:\Users\lenovo\.julia\packages\ArrayInterface\ubN3d\src\ArrayInterface.jl:329
  ...
Stacktrace:
 [1] __init(prob::ODEProblem{VectorOfArray{Float64, 3, Vector{Matrix{Float64}}},...
1 Like
4

Open an issue, that’s unintended and fixable.

1 Like