Error when using LawsonEuler method on SplitODE

Specific Domains Modelling & Simulations
1

Hi all. I have a Split/Semiinear ODE I want to solve efficiently. The docs suggest the LawsonEuler algorithm. I tried to implement but am running into an error I don’t understand. Here is a MWE

using DifferentialEquations
using LinearAlgebra

function lineardynamics!(du, u, p, t)
    mul!(du, p[1], u)
end

function nonlineardynamics!(du, u, p, t)
    broadcast!(*, du, u, u)
end

n = 10
u0 = zeros(n)
hvals = Diagonal(rand(n));
tspan = (0, 1);

p = (hvals)
prob = SplitODEProblem(lineardynamics!, nonlineardynamics!, u0, tspan, p);
sol = solve(prob, alg=LawsonEuler(), dt=1e-2)

This gives the following MethodError

MethodError: no method matching size(::typeof(lineardynamics!))

Closest candidates are:
  size(::Union{QR, LinearAlgebra.QRCompactWY, QRPivoted})
   @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/qr.jl:582
  size(::Union{QR, LinearAlgebra.QRCompactWY, QRPivoted}, ::Integer)
   @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/qr.jl:581
  size(::Union{LinearAlgebra.QRCompactWYQ, LinearAlgebra.QRPackedQ})
   @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/qr.jl:585
  ...


Stacktrace:
  [1] alg_cache(alg::LawsonEuler{10, true, Val{:forward}, true, nothing}, u::Vector{Float64}, rate_prototype::Vector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::Vector{Float64}, uprev2::Vector{Float64}, f::SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(nonlineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, t::Float64, dt::Float64, reltol::Float64, p::Diagonal{Float64, Vector{Float64}}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/NBaQM/src/caches/linear_nonlinear_caches.jl:194
  [2] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Diagonal{Float64, Vector{Float64}}, SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(nonlineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, alg::LawsonEuler{10, true, Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Int64, abstol::Nothing, reltol::Nothing, qmin::Int64, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Int64, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, LawsonEuler{0, true, Val{:forward}, true, nothing}, Tuple{Symbol}, NamedTuple{(:alg,), Tuple{LawsonEuler{0, true, Val{:forward}, true, nothing}}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/NBaQM/src/solve.jl:341
  [3] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/NBaQM/src/solve.jl:10 [inlined]
  [4] #__solve#746
    @ ~/.julia/packages/OrdinaryDiffEq/NBaQM/src/solve.jl:5 [inlined]
  [5] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/NBaQM/src/solve.jl:1 [inlined]
  [6] solve_call(_prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Diagonal{Float64, Vector{Float64}}, SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(nonlineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, args::LawsonEuler{10, true, Val{:forward}, true, nothing}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:alg, :dt), Tuple{LawsonEuler{0, true, Val{:forward}, true, nothing}, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:608
  [7] solve_up(::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Diagonal{Float64, Vector{Float64}}, SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(nonlineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, ::Nothing, ::Vector{Float64}, ::Diagonal{Float64, Vector{Float64}}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:alg, :dt), Tuple{LawsonEuler{0, true, Val{:forward}, true, nothing}, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:1057
  [8] solve_up
    @ ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:1043 [inlined]
  [9] solve(::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Diagonal{Float64, Vector{Float64}}, SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(nonlineardynamics!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:alg, :dt), Tuple{LawsonEuler{0, true, Val{:forward}, true, nothing}, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:980
 [10] top-level scope
    @ In[3]:1
Selection deleted

I’m using Julia 1.9.3 and DifferentialEquation 7.12.0. I’m not sure where I’ve gone wrong or why size is being called on lineardynamics!. Any help is much appreciated!

2

Solved. Instead of passing in lineardynamics! I pass in a MatrixOperator as the first function in the ODEProblem. The MatrixOperator defines the linear dynamics.

using DifferentialEquations
using LinearAlgebra
using SciMLOperators

function nonlineardynamics!(du, u, p, t)
    broadcast!(*, du, u, u)
end

n = 10
u0 = rand(n)
hvals = Diagonal(rand(n));
H = MatrixOperator(hvals);
tspan = (0, 1);

prob = SplitODEProblem(M, nonlineardynamics!, u0, tspan);
sol = solve(prob, alg=LawsonEuler(), dt=1e-2)

Helpful example in the docs here