Solving a large system of differential equations

Specific Domains Numerics
,
1

Hey :slight_smile: ,
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

Would be really happy about any help you can give me.

2

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:

https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Defining-Linear-Solver-Routines-and-Jacobian-Free-Newton-Krylov

https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Sundials-Specific-Handling

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

3

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?

4

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?

5

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.

6

sol = solve(prob, Tsit5(), adaptive=true)
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]
[6] similar at .\broadcast.jl:197 [inlined]
[7] similar at .\broadcast.jl:196 [inlined]
[8] similar at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\structuredbroadcast.jl:130 [inlined]
[9] copy at .\broadcast.jl:862 [inlined]
[10] materialize at .\broadcast.jl:837 [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

7

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.

8

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)

9

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.

10

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.

11

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