# OrdinaryDiffEq with a mass\_matrix and jac\_prototype

**URL:** https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462
**Category:** Numerics
**Tags:** ordinarydiffeq
**Created:** [December 10, 2025, 7:56am UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462 "2025-12-10T07:56:50Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)
#### Post date: [December 10, 2025, 7:56am UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462/1 "2025-12-10T07:56:51Z")

</div>

Hi everyone,

Short story : see the MWE below.

Using `OrdinaryDiffEq`, I’m trying to solve an ODE involving a `mass_matrix` that (may) evolve depending on the solution, the time, the parameters. So I set up a `SciMLOperators.MatrixOperator` and provided it to the `ODEFunction`. So far, so good.

The problem I’m facing is that if I provide a `jac_prototype` (of my rhs function) to the `ODEFunction`, then the `solve` errors here [OrdinaryDiffEq.jl/lib/OrdinaryDiffEqDifferentiation/src/alg\_utils.jl at b80bd3b002927c63ba1d0eec9fdcbcee3d560873 · SciML/OrdinaryDiffEq.jl · GitHub](https://github.com/SciML/OrdinaryDiffEq.jl/blob/b80bd3b002927c63ba1d0eec9fdcbcee3d560873/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl#L120). Without `jac_prototype`, on a MWE, everything is fine. Am I doing something wrong?

Moreover, I don’t understand the quoted lines in `OrdinaryDiffEqDifferentiation` : the `jac_prototype` is related to my `rhs` function, why are `mass_matrix` and `jac_prototype` connected together here?

Here is a MWE:

```julia
using OrdinaryDiffEq
using SparseArrays
using SparseDiffTools

const M_mwe = sparse([1.0 2.0; 0.0 1.0])

function update_func(A, u, p, t)
    return M_mwe
end

function run()
    mass_matrix = MatrixOperator(M_mwe; update_func = update_func)
    f!(du, u, p, t) = du .= u
    alg = ImplicitEuler()
    tspan = (0.0, 1.0)
    u0 = ones(2)

    # Without jac_prototype
    odeFunction = ODEFunction(f!; mass_matrix)
    prob = ODEProblem(odeFunction, u0, tspan)
    solve(prob, alg) #-------> OK

    # With jac_prototype
    jac_prototype = sparse([1.0 0.0; 0.0 1.0])
    colorvec = [1, 1]
    odeFunction = ODEFunction(f!; mass_matrix, jac_prototype, colorvec)
    prob = ODEProblem(odeFunction, u0, tspan)
    solve(prob, alg) #-------> KO
end

run()

```

PS : passing the `jac_prototype` of the rhs is not an option for me because in my real case the automatic tools (such as `SparseConnectivityTracer`) fail to find it.

---

<div class="post-metadata">

### Author: ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [December 10, 2025, 8:26am UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462/2 "2025-12-10T08:26:51Z")

</div>

Hi!  
Can you give the error stack trace?  
Also, what do you get when you try automatic sparsity detection? Just for curiosity

---

<div class="post-metadata">

### Author: ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)
#### Post date: [December 10, 2025, 10:17am UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462/3 "2025-12-10T10:17:08Z")

</div>

Hi Guillaume, here is the stack trace:

```julia
MethodError: no method matching keys(::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
The function `keys` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  keys(::Cmd)
   @ Base process.jl:716
  keys(::Base.TermInfo)
   @ Base terminfo.jl:232
  keys(::LibGit2.GitTree)
   @ Revise .julia/packages/Revise/p9Zlq/src/git.jl:52
  ...

Stacktrace:
  [1] pairs(collection::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
    @ Base ./abstractdict.jl:178
  [2] findall(testf::ComposedFunction{typeof(!), typeof(iszero)}, A::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
    @ Base ./array.jl:2700
  [3] prepare_user_sparsity(ad_alg::AutoForwardDiff{nothing, ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem})
    @ OrdinaryDiffEqDifferentiation .julia/packages/OrdinaryDiffEqDifferentiation/DE4z5/src/alg_utils.jl:120
  [4] prepare_alg(alg::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}, u0::Vector{Float64}, p::SciMLBase.NullParameters, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem})
    @ OrdinaryDiffEqDifferentiation .julia/packages/OrdinaryDiffEqDifferentiation/DE4z5/src/alg_utils.jl:53
  [5] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}; kwargs::@Kwargs{})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1200
  [6] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1183
  [7] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::@Kwargs{})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1096
  [8] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1086
  [9] run()
    @ Main ./REPL[11]:18
 [10] top-level scope
    @ REPL[12]:1

```

Regarding the automatic sparsity detection, the problem is this one : [Bcube and DifferentiationInterface · Issue #189 · bcube-project/Bcube.jl · GitHub](https://github.com/bcube-project/Bcube.jl/issues/189#issuecomment-3273798705). I could use your trick, but as a first step I’d prefer not to.

---

<div class="post-metadata">

### Author: ![jClugstor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jclugstor/32/32689_2.png) [@jClugstor](https://discourse.julialang.org/u/jClugstor)
#### Post date: [December 10, 2025, 2:47pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462/4 "2025-12-10T14:47:44Z")

</div>

If I recall correctly, this probably has something to do with the fact that the `jac_prototype` and sparsity pattern need to have certain entries that correspond to the non-zeroes of the mass matrix set to be potentially non-zero. It has something to do with the way the WOperator is calculated.

Probably a related issue: [Strict in-place sparse Jacobians now cause errors when they don't anticipate the mass matrix · Issue #2653 · SciML/OrdinaryDiffEq.jl · GitHub](https://github.com/SciML/OrdinaryDiffEq.jl/issues/2653)

Can you file an issue on the github repo?

---

<div class="post-metadata">

### Author: ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)
#### Post date: [December 10, 2025, 5:13pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-with-a-mass-matrix-and-jac-prototype/134462/5 "2025-12-10T17:13:08Z")

</div>

Yes, I read again the fonction and indeed the jacobian of the whole system `M dq/dt=f(x)` is something like `M - dt J` xhere `J=df/dq`.

Ok I’ll open an issue on the repo with this MWE and in the meantime I will write the specialization for `findall(!iszero, M::MatrixOperator` in my script.