Interpolation of solution fails when mass matrix is the identity matrix

[natlampen]

Interpolation using the solution struct errors when the mass matrix is the identity matrix.

A bit of context. For my use case I was passing some mass information to a function which would be passed to the mass matrix.

Minimal Reproducible Example 👇

using OrdinaryDiffEq 

function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = 1
end

fun = ODEFunction(f!; mass_matrix = Float64[1 0; 0 1])
prob = ODEProblem(fun, [0, 0], (0, 10))
sol = solve(prob)
sol(10)

Error & Stacktrace ⚠️

ERROR: LoadError: Derivative order too high for interpolation order. An interpolation derivative is
only accurate to a certain derivative. For example, a second order interpolation
is a quadratic polynomial, and thus third derivatives cannot be computed (will be
incorrectly zero). Thus use a solver with a higher order interpolation or compute
the higher order derivative through other means.

You can find the list of available ODE/DAE solvers with their documented interpolations at:

* https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/
* https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Stacktrace:
  [1] _ode_interpolant!(out::Vector{…}, Θ::Float64, dt::Float64, y₀::Vector{…}, y₁::Vector{…}, k::Vector{…}, cache::OrdinaryDiffEq.Tsit5Cache{…}, idxs::Nothing, T::Type{…}, differential_vars::BitVector)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/dense/interpolants.jl:26
  [2] ode_interpolant(Θ::Float64, dt::Float64, y₀::Vector{…}, y₁::Vector{…}, k::Vector{…}, cache::OrdinaryDiffEq.Tsit5Cache{…}, idxs::Nothing, T::Type{…}, differential_vars::BitVector)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/dense/generic_dense.jl:945
  [3] ode_interpolation(tval::Int64, id::OrdinaryDiffEq.InterpolationData{…}, idxs::Nothing, deriv::Type{…}, p::SciMLBase.NullParameters, continuity::Symbol)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/dense/generic_dense.jl:784
  [4] InterpolationData
    @ ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/interp_func.jl:169 [inlined]
  [5] AbstractODESolution
    @ ~/.julia/packages/SciMLBase/rR75x/src/solutions/ode_solutions.jl:186 [inlined]
  [6] #_#473
    @ ~/.julia/packages/SciMLBase/rR75x/src/solutions/ode_solutions.jl:177 [inlined]
  [7] AbstractODESolution
    @ ~/.julia/packages/SciMLBase/rR75x/src/solutions/ode_solutions.jl:175 [inlined]
  [8] (::ODESolution{…})(t::Int64)
    @ SciMLBase ~/.julia/packages/SciMLBase/rR75x/src/solutions/ode_solutions.jl:175
  [9] top-level scope
    @ ~/Development/DiffEqBug/test.jl:11
 [10] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [11] top-level scope
    @ REPL[3]:1
in expression starting at /home/natlampen/Development/DiffEqBug/test.jl:11
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [1dea7af3] OrdinaryDiffEq v6.85.0

• Output of versioninfo()

Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 24 × AMD Ryzen 9 5900X 12-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 24 virtual cores)

[ChrisRackauckas]

I see... you should be writing this as fun = ODEFunction(f!; mass_matrix = LinearAlgebra.I) because allocating that mass matrix is not useful for performance, and if you use I then that's the same as the default. Right now we're checking for compatibility with === I, but we probably should be checking if the matrix is equal to the identity matrix in general.

However, is there a reason to not use I here?

[natlampen]

Yes. That is a good point. This was the most basic example I could make trigger the error. This was an edge case of a spring-mass-damper system where mass is 1. So essentially the mass matrix would be LinearAlgebra.Diagonal([1,m]).

[oscardssmith]

we should be able to support arbitrary nonzero diagonal matrices here.