# How many differential equations can I solve simultaneously? Julia is crashing for more than a few thousand odes

**URL:** <https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [February 1, 2022, 5:20pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588 "2022-02-01T17:20:49Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 1, 2022, 5:20pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/1 "2022-02-01T17:20:49Z")

</div>

I was trying to solve a discretized partial differential equation when I realized that Julia is crashing when I try to solve more than a few thousand differential equations. When I try to run the following code, I get a message that says Julia “Julia has exited. Press Enter to start a new session.” in Juno and if I run the same code in the terminal then I simply get an output saying “Killed”. If I change the grid to 20 x 20 the code runs, so there is no issue with the code. Is there a trick to solving more differential equations in Julia?

```julia

using DifferentialEquations
using Sundials

function eoms!(du,u,p,t)
    println(t)
    for i in 1:200
        for j in 1:200
            du[i,j]=-1.0*u[i,j]
        end
    end
end

u0=ones(Float64,(200,200))
tspan=(0.0,1.0e-7)
p=[]
prob=ODEProblem(eoms!,u0,tspan,p)

tol=1e-2
reltolp=tol
abstolp=tol

@time sol=solve(prob,KenCarp4(),abstol=abstolp,reltol=reltolp,saveat=1e-7)

using Plots
plot(sol,vars=(0,1))

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 1, 2022, 5:44pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/2 "2022-02-01T17:44:18Z")

</div>

Probably depends on the RAM of your computer. With 64 GB on mine the MWE is barely running.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 1, 2022, 5:49pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/3 "2022-02-01T17:49:38Z")

</div>

But the same code works when I implement it in c++; even if I increase the dimensions to 2000 \times 2000. So if it is a RAM issue there must be a way around it because the problem itself does not require that much RAM.

Edit: for 100 \times 100, the allocation is 1.90 M. It cant be running out of RAM at 200 \times 200, if it is easily running for 100 \times 100 grid.

---

<div class="post-metadata">

**Author:** ![Keno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/keno/32/285_2.png) [@Keno](https://discourse.julialang.org/u/Keno)\
**Post date:** [February 1, 2022, 6:00pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/4 "2022-02-01T18:00:01Z")

</div>

I don’t believe the low level interfaces are automatically computing the jacobian sparsity. You’ll either need to declare it manually (easy in your case, since the equations are not coupled) or use one of the automatic sparsity detection packages. Or just use one of the higher level packages like MTK that does this automatically.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 1, 2022, 6:07pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/5 "2022-02-01T18:07:07Z")

</div>

I have posted the minimal example that reproduces the problem. The actual problem I am trying to solve is coupled and the jacobian is not very sparse. It is sparse for the first index but not the second.

Edit: You were right though. If I change the linear solver to GMRES. Then the code runs without any problem for 2000 \times 2000 system.

---

<div class="post-metadata">

**Author:** ![Keno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/keno/32/285_2.png) [@Keno](https://discourse.julialang.org/u/Keno)\
**Post date:** [February 1, 2022, 6:15pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/6 "2022-02-01T18:15:34Z")

</div>

> [@Shashank](#):
>
> The actual problem I am trying to solve is coupled and the jacobian is not very sparse.

Ok, but a dense jacobian for a 2000x2000 system is 128TB, which though in reach for an HPC system is not something you’d run on a single node, so if that’s the size you care about I’d take another look at the sparsity of your jacobian.

---

<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:** [February 1, 2022, 6:44pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/7 "2022-02-01T18:44:34Z")

</div>

> [@Shashank](#):
>
> The actual problem I am trying to solve is coupled and the jacobian is not very sparse.

A partial differential equation that doesn’t give a sparse Jacobian? I’ve never heard of that. The quintessential example of sparse matrices is PDEs.

Anyways, here’s a sparse PDE code optimization example.

[https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/)

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 1, 2022, 6:58pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/8 "2022-02-01T18:58:16Z")

</div>

Thanks for sharing the example. Actually, the problem I am trying to solve is an integro-differential equation. So the pde part of it is sparse but I have integration (or sum in the discretized setup) over the second index.

---

<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:** [February 2, 2022, 4:05am UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/9 "2022-02-02T04:05:25Z")

</div>

Then you might need to develop a matrix-free preconditioner

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 2, 2022, 12:03pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/10 "2022-02-02T12:03:59Z")

</div>

I don’t know how to do that. Can you point me towards an example on the internet where this is done for simple differential equations? If I can understand the syntax I can modify it for my case.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 2, 2022, 12:35pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/11 "2022-02-02T12:35:37Z")

</div>

Can’t tell if this is satire or not ^^

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 2, 2022, 12:51pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/12 "2022-02-02T12:51:02Z")

</div>

> [@josuagrw](#):
>
> > [@ChrisRackauckas](#):
> >
> > Then you might need to develop a matrix-free preconditioner
> 
> Can’t tell if this is satire or not ^^

It’s not satire.

A “matrix-free” method is one in which you act a linear operator on a vector without _explicitly_ constructing a matrix. For example, a you can do a convolution in O(n \log n) time and O(n) memory via FFTs, whereas an explicit matrix would be O(n^2). For integral operators — relevant to the present case with an integro-differential equation — there are things like the fast multipole method and related tricks (e.g. hierarchical low-rank approximations) to multiply quickly and with much less memory than an explicit matrix.

If you are doing iterative solves on each step for a stiff system, then you probably additionally want a preconditioner, and again you may want to use a matrix-free method if your explicit matrices are not sparse.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 2, 2022, 1:13pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/13 "2022-02-02T13:13:04Z")

</div>

Nope, and here is an example implementation:

[https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/#Adding-a-Preconditioner](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Adding-a-Preconditioner)

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 2, 2022, 1:20pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/14 "2022-02-02T13:20:45Z")

</div>

Thank you — I actually meant to respond to a different comment 🙈

---

<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:** [February 2, 2022, 1:48pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/15 "2022-02-02T13:48:49Z")

</div>

So to start, there are two other methods I thought of.

One is to build the preconditioner based on an approximation. If you do sparsity detection without the integral, you’ll get a sparse matrix that approximates your Jacobian. You can then use that inside of DifferentialEquations.jl as the `jac_prototype` and that will generate an approximate Jacobian through color differentiation to that sparse matrix. That approximate Jacobian can be used inside of an ILU/AMG preconditioner. This is fine because preconditioners are only supposed to be an approximation anyways: whether this works is really whether it’s a “good enough” approximation, which depends on your equation and the parameter values. This only takes like 5 lines of code to try (copy-pasted out of the tutorial) so I would try this first.

The second thing you can do is use an IMEX integrator:

[https://diffeq.sciml.ai/stable/types/split\_ode\_types/](https://diffeq.sciml.ai/stable/types/split_ode_types/)  
[https://diffeq.sciml.ai/stable/solvers/split\_ode\_solve/#OrdinaryDiffEq.jl](https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl)

I’d recommend KenCarp47, where you split the integral to be part of the explicit part. If the integral has well-behaved eigenvalues, this would let you treat it explicitly, and then your true Jacobian is all of the non-integral stuff, and then once again sparse matrix preconditioners etc. etc. that tutorial applies.

One of those two “should” work. But let’s address the other thing just for full clarity:

> [@Shashank](#):
>
> Nope, and here is an example implementation:
> 
> [https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/#Adding-a-Preconditioner](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Adding-a-Preconditioner)

Note that those preconditioners are not matrix-free. ILU needs to know the sparse matrix and a `tau` cutoff value in order to construct the preconditioner. Algebraic Multigrids also need to know the structure of the equations through the sparse matrix. That’s really the issue: preconditioners need to match details of the equation since they need to be an approximate inverse, so the easiest way to construct an algorithm that knows “some structure of the equations” is to have as input a sparse matrix because that naturally will have much of the structure for many PDEs.

That said, as @stevengj discussed, there are preconditioners that do not require you have the sparse matrix. All of those will naturally require that you write a function based solely on your PDE, for example fast multiple methods and geometric multigrids, so it’s hard to ever give someone a software that does this automatically like you can with AMG and ILU (well, there is a way to do it, but that’ll be a research project over the next year 😉).

I think the easiest to write might be a geometric multigrid, since it just requires knowing how to write your PDE with a different dx, along with expanding and contracting the grid. Gil has some good notes on it:

> **[Readings | Mathematical Methods for Engineers II | Mathematics | MIT...](https://ocw.mit.edu/courses/18-086-mathematical-methods-for-engineers-ii-spring-2006/pages/readings/)**
>
> The readings section provides information on textbooks, and supplementary readings for the course.

But other than that, for now it’s kind of “may the force be with you” in terms of automation, which I hope can be solved better in the future.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [February 2, 2022, 2:15pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/16 "2022-02-02T14:15:15Z")

</div>

Awesome! Thanks for the help.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 2, 2022, 2:42pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/17 "2022-02-02T14:42:18Z")

</div>

For the record, boundary element methods (BEM) produce PDEs with dense Jacobian

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 2, 2022, 3:08pm UTC](https://discourse.julialang.org/t/how-many-differential-equations-can-i-solve-simultaneously-julia-is-crashing-for-more-than-a-few-thousand-odes/75588/18 "2022-02-02T15:08:36Z")

</div>

> [@josuagrw](#):
>
> For the record, boundary element methods (BEM) produce PDEs with dense Jacobian

I would say that BEM solves a surface-integral equation (SIE), not a PDE.

(The SIE is generally _equivalent_ to a PDE, so in that sense you are also solving a PDE. But it is the integral and not the differential equation that is the starting point for the numerics, and nontrivial computation is required to get the PDE solution at an arbitrary point from the SIE solution. Conversely, traditional FEM solves a PDE, which of course also gives you the corresponding SIE solution, but most people would not describe FEM as an “integral-equation solver”.)

Put another way, Galerkin BEM and FEM both assemble a matrix A\_{ij} = \langle b\_i, \hat{L} b\_j \rangle where b\_i and b\_j are localized polynomial basis functions on a mesh (a surface mesh in BEM and a volume mesh in FEM) and \langle \cdot, \cdot \rangle is an inner product (integral), but the linear operator \hat{L} is a differential (PDE) operator for FEM that yields a sparse A, and \hat{L} is an integral (SIE) operator for BEM that yields a dense A (which you represent in a matrix-free manner for large-scale problems).