DifferentialEquations.jl solution object: save_idxs and interpolation

Specific Domains Modelling & Simulations
1

Hello,
I’m trying to solve a largeish number of largeish systems of linear differential equations. I ultimately only care about one component from each equation, so save_idxs in principle should be exactly what I want: I get all the goodies of a solution object (e.g. interpolation), without having to save all the data.

Problem is, interpolation isn’t working! The following code snippet illustrates the problem I’m trying to solve, and the error I’m getting. The analogue without save_idxs works fine. I’m on Julia 1.4.2 and DifferentialEquations.jl 6.15.

Code

using DifferentialEquations
N = 10

H = randn(N,N)
H -= H'

ψ0 = zeros(N); ψ0[1] = 1

f(ψ :: Array{Float64,1},p,t) = H*ψ

prob = ODEProblem(f, ψ0, (0.0,1.0), save_idxs=[1])
sol = solve(prob, reltol=1e-8, abstol=1e-8)
sol(0.5)

Error:

DimensionMismatch("second dimension of A, 10, does not match length of x, 1")

Stacktrace:
 [1] gemv!(::Array{Float64,1}, ::Char, ::Array{Float64,2}, ::Array{Float64,1}, ::Bool, ::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:456
 [2] mul! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:66 [inlined]
 [3] mul! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:208 [inlined]
 [4] *(::Array{Float64,2}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:47
 [5] f(::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64) at ./In[14]:9
 [6] ODEFunction at /home/christopher/.julia/packages/DiffEqBase/V7P18/src/diffeqfunction.jl:248 [inlined]
 [7] addsteps!(::Array{Array{Float64,1},1}, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::DiffEqBase.NullParameters, ::OrdinaryDiffEq.Vern7ConstantCache{OrdinaryDiffEq.Vern7Tableau{Float64,Float64}}, ::Bool, ::Bool, ::Bool) at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/dense/verner_addsteps.jl:281
 [8] addsteps!(::Array{Array{Float64,1},1}, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::Function, ::DiffEqBase.NullParameters, ::OrdinaryDiffEq.Vern7ConstantCache{OrdinaryDiffEq.Vern7Tableau{Float64,Float64}}) at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/dense/verner_addsteps.jl:266
 [9] ode_interpolation(::Float64, ::Function, ::Nothing, ::Type{T} where T, ::DiffEqBase.NullParameters, ::Symbol) at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/dense/generic_dense.jl:281
 [10] CompositeInterpolationData at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/interp_func.jl:77 [inlined]
 [11] #_#388 at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/composite_solution.jl:16 [inlined]
 [12] ODECompositeSolution at /home/christopher/.julia/packages/OrdinaryDiffEq/VPJBD/src/composite_solution.jl:16 [inlined] (repeats 2 times)
 [13] top-level scope at In[14]:13
2

Meant to add & posted early—I must be missing something v. basic. What’s the best way to get access to those interpolations? Should I just specify the points I want ahead of time with saveat?

3 
using DifferentialEquations
N = 10

H = randn(N,N)
H -= H'

ψ0 = zeros(N); ψ0[1] = 1

f(ψ :: Array{Float64,1},p,t) = H*ψ

prob = ODEProblem(f, ψ0, (0.0,1.0), save_idxs=[1])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
sol(0.5)

is fine. Or using saveat is fine. Vern7(lazy=false) would probably be a good choice here. It’s an obscure issue that comes from the exact combination of features you chose and the algorithm that got automatically picked. It’s kind of amazing you found it :wink:. We’ll try to make this into account in the default algorithm chooser.

4

Here I was thinking I was doing the just about the most straightforward possible thing. Specifying Tsit5 works like a charm—thank you!