# Solving a large system of differential equations

Hey ,
I am currently trying to solve a system of differential equations. I create the coefficient matrix as a sparse matrix, but when i try to use the ode solver, it gives me an out of memory error. The code is:

``````ff(u,p,t) = (L_s + I(63504).*p*Î±*cos(Ď‰ * t)) * u
prob = ODEProblem(ff, collect(Complex{Float64}, u_0), time, p)
sol = solve(prob, AutoTsit5(ABDF2()), adaptive=true, saveat = 0.1)
``````

Where p is a vector with complex integer values (not sure how important the actual vector is but hereâ€™s the code)

``````driving = [false true true true false]
for (i, (state, state_t)) in enumerate(basis_states_s)
particles = sum([parse(Int, i) for i in [x for x in state][collect(Iterators.flatten([[x x] for x in driving]))]])
particles_t = sum([parse(Int, i) for i in [x for x in state_t][collect(Iterators.flatten([[x x] for x in driving]))]])
p[i] = -1im*(particles - particles_t)
end
``````

Link to the â€śL_sâ€ť matrix â†’ iCloud

Is your system stiff? `AutoTsit5(ABDF2())` is almost certainly a bad idea (even for stiff), but if itâ€™s non-stiff you can cut out the generated Jacobian by just using the non-stiff ODE solver. If it is stiff, well your memory issue is that the Jacobian requires n^2 memory so youâ€™ll want to use a Jacobian-free method. One easy thing is ROCK2, but you can also follow the tutorial on how to choose Jacobian-free methods:

Additionally, defining sparse Jacobians helps both memory and efficiency a ton, so you might want to specify a sparsity pattern.

I dont know if it is stiff our not i got to admit, but when just using Tsit5(), or any other non-stiff solver i tried so far mentioned here ODE Solvers Â· DifferentialEquations.jl it ran out of memory as well. So if I use a non-stiff solver and it still runs out of memory then what is the issue there?

Tsit5 will use 7*u0 memory, so are you close to the memory limit before solving? You can use a low-memory RK method if you need to.

https://diffeq.sciml.ai/stable/solvers/ode_solve/#Low-Storage-Methods

Is `L_s` being read into a large dense matrix?

Well actually it gives me the out of memory error after a few seconds, or if i run the command with the same solver again, instantly. It does even show a spike in my memory usage when i check it in the task manager.
â†’ I also posted the error message.
So L_s is certainly a sparse matrix, because otherwise it would run out of memory even before i use the ode solver.

ERROR: OutOfMemoryError()
Stacktrace:
[1] Array at .\boot.jl:424 [inlined]
[2] Array at .\boot.jl:432 [inlined]
[3] Array at .\boot.jl:439 [inlined]
[4] similar at .\abstractarray.jl:675 [inlined]
[5] similar at .\abstractarray.jl:674 [inlined]
[11] ff(::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}, ::Float64) at .\none:1
[12] ODEFunction at C:\Users\Daniel.julia\packages\SciMLBase\UIp7W\src\scimlfunctions.jl:334 [inlined]
[13] initialize!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,false,Array{Complex{Float64},1},Nothing,Float64,Array{Complex{Float64},1},Float64,Float64,Float64,Float64,Array{Array{Complex{Float64},1},1},ODESolution{Complex{Float64},2,Array{Array{Complex{Float64},1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Complex{Float64},1},1},1},ODEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},false,Array{Complex{Float64},1},ODEFunction{false,typeof(ff),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{false,typeof(ff),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Array{Array{Complex{Float64},1},1},Array{Float64,1},Array{Array{Array{Complex{Float64},1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},DiffEqBase.DEStats},ODEFunction{false,typeof(ff),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,PIController{Rational{Int64}},typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),Nothing,CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.FasterForward},DataStructures.BinaryHeap{Float64,DataStructures.FasterForward},Nothing,Nothing,Int64,Tuple{},Tuple{},Tuple{}},Array{Complex{Float64},1},Complex{Float64},Nothing,OrdinaryDiffEq.DefaultInit}, ::OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}) at C:\Users\Daniel.julia\packages\OrdinaryDiffEq\PIjOZ\src\perform_step\low_order_rk_perform_step.jl:565
[14] __init(::ODEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},false,Array{Complex{Float64},1},ODEFunction{false,typeof(ff),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem}, ::Tsit5, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::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::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{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), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\Daniel.julia\packages\OrdinaryDiffEq\PIjOZ\src\solve.jl:456
[15] #__solve#465 at C:\Users\Daniel.julia\packages\OrdinaryDiffEq\PIjOZ\src\solve.jl:4 [inlined]
[16] #solve_call#56 at C:\Users\Daniel.julia\packages\DiffEqBase\QiFNl\src\solve.jl:61 [inlined]
[17] solve_up(::ODEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},false,Array{Complex{Float64},1},ODEFunction{false,typeof(ff),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem}, ::Nothing, ::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}, ::Tsit5; kwargs::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:adaptive,),Tuple{Bool}}}) at C:\Users\Daniel.julia\packages\DiffEqBase\QiFNl\src\solve.jl:82
[18] #solve#57 at C:\Users\Daniel.julia\packages\DiffEqBase\QiFNl\src\solve.jl:70 [inlined]
[19] top-level scope at none:1

Do you need to save every step of the solution in a way that is dense and continuous? If not, use `saveat`, `save_everystep`, etc. to control the saving behavior.

I tried saveat, but it does not work as well. So I used really huge time step over small overall time (e.g. time = (0, 10), and saveat = 0.5)

Itâ€™s not a time step, itâ€™s a save point. The main difference being, using a large `saveat` wonâ€™t induce any error so itâ€™s safe. If youâ€™re running out of memory and thatâ€™s as much as your computer can save then thatâ€™s probably a good option. The other option would be to kick out parts of the solution to disk while solving.

Sorry I misspoke â†’ I have figured out already that there is a difference between timestep and saveat, and i did actually use saveat, I just wrote it wrong before. But it just seems odd to me that i get an error very fast, without seeing a continuus rise in memory usage. It just tells me quite fast it ran out of memory.

Yeah thatâ€™s odd because you should see it allocating one after another until it fills up?