I wrote an ad hoc script for a method of lines problem where I discretized the spatial domain and used and ODE solver to integrate in the time domain. At the moment, I’m not interested in a discussion of whether this is the best way to solve a given PDE system, more so a discussion of the issue that I ran into.
I have several variables that are functions of both time and position, and my initial instinct was to structure the problem as one where u is a matrix, and each column of u represents the value of each variable at different points in the spatial discretization at the current time step of the ODE solver.
In order to satisfy various boundary conditions, I needed to set this up as a DAE system where algebraic equations are used to impose a specific constraint at various nodes. If u is a vector, this works fine. If u is a matrix, I get an error that I suspect is arising from using a mass matrix with a system where du is also a matrix, even though, in theory, this should give the correct result as long as each column of du has differential and algebraic equations at the same indices. See example below for an illustration of the type of error I am seeing.
##Example of an ODEsystem with u in matrix form. This will run as expected.
function odefun!(du, u, p, t)
du[:,:] = -1*u
end
u0 = [1 2; 3 4]
span = (0.0, 10.0)
f = ODEFunction(odefun!)
prob = ODEProblem(f, u0, span)
sol = solve(prob)
plt1 = plot(sol)
display(plt1)
#Analogous example of a DAE system with u in matrix form. This will cause an exception
function daefun!(du, u, p, t)
du[:,:] = -1*u
du[end,:] = 10 .- u[end,:]
end
u0 = [1 2; 3 4] #will work for u0 = [1, 2]
span = (0.0, 10.0)
M = [1 0; 0 0] #In my mind this says, first variable is differential, second is algebraic
f2 = ODEFunction(daefun!, mass_matrix = M)
prob = ODEProblem(f2, u0, span)
sol = solve(prob, Rodas4())
plt2 = plot(sol)
display(plt2)
Stacktrace is a bit long, but it kicks off with this:
BoundsError: attempt to access 2×2 Matrix{Float64} at index [Bool[0, 1]]
Stacktrace:
[1] throw_boundserror(A::Matrix{Float64}, I::Tuple{Base.LogicalIndex{Int64, Vector{Bool}}})
@ Base .\abstractarray.jl:651
[2] checkbounds
@ .\abstractarray.jl:616 [inlined]
[3] view
@ .\subarray.jl:177 [inlined]
[4] _initialize_dae!(integrator:
I could alternatively construct this problem with u0 = [1, 2, 3, 4], and M = diagm([1,0,1,0]) which would resolve the issue, but I’m curious why this doesn’t work. If you multiply the mass matrix above by the du in the DAE system, it seems to generate the result that I want, but that doesn’t necessarily mean that it is behaving correctly with the integrator.