# Solving a large system of differential equations

**URL:** <https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870>\
**Category:** Numerics\
**Tags:** question, diffeq\
**Created:** [September 28, 2021, 11:32am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870 "2021-09-28T11:32:23Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 28, 2021, 11:32am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/1 "2021-09-28T11:32:23Z")

</div>

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](https://www.icloud.com/iclouddrive/0v6X9LNJkwihIbebdkbp0IbOQ#L_s)

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 28, 2021, 11:43am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/2 "2021-09-28T11:43:48Z")

</div>

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/#Defining-Linear-Solver-Routines-and-Jacobian-Free-Newton-Krylov)

[https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/#Sundials-Specific-Handling](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.

---

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 28, 2021, 12:03pm UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/3 "2021-09-28T12:03:47Z")

</div>

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](https://diffeq.sciml.ai/stable/solvers/ode_solve/) 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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 28, 2021, 12:13pm UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/4 "2021-09-28T12:13:43Z")

</div>

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](https://diffeq.sciml.ai/stable/solvers/ode_solve/#Low-Storage-Methods)

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

---

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 28, 2021, 12:17pm UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/5 "2021-09-28T12:17:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 28, 2021, 12:18pm UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/6 "2021-09-28T12:18:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 28, 2021, 1:52pm UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/7 "2021-09-28T13:52:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 29, 2021, 9:42am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/8 "2021-09-29T09:42:57Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 29, 2021, 10:08am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/9 "2021-09-29T10:08:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![danielw96](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@danielw96](https://discourse.julialang.org/u/danielw96)\
**Post date:** [September 29, 2021, 10:35am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/10 "2021-09-29T10:35:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 29, 2021, 10:50am UTC](https://discourse.julialang.org/t/solving-a-large-system-of-differential-equations/68870/11 "2021-09-29T10:50:18Z")

</div>

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