Differential Equations.jl ODEProblem exceedingly slow

I am trying to numerically integrate a set of master equations with DifferentialEquations.jl

I have defined the differential equation like:

function voter_model_ame!(dG,G,t,p)
    N,U  = size(G)
    for n_prime=1:N,u_prime=1:U
        u  = u_prime-1
        n = n_prime-1
        dG[n_prime,u_prime] = 0

        if n == 0

        if u > n 

        (u < n) && (dG[n_prime,u_prime] -= G[n_prime,u_prime]*(n-u)*(u/n))#upflip outflux
        (u > 0) && (dG[n_prime,u_prime] -= G[n_prime,u_prime]*(u)*((n-u)/n))#downflip outflux
        (u < n) && (dG[n_prime,u_prime]+= G[n_prime,u_prime+1]*((u+1)*((n-u-1)/n)))#inwards downflip flux
        (u > 0) && (dG[n_prime,u_prime]+= G[n_prime,u_prime-1]*((n-u+1)*((u-1)/n)))#inwards upflip flux
    return dG 

This is a system of N^2 differential equations.

I can set the parameter N = 500 and I try to numerically integrate like this:

t_max = 10
p = [0]
u₀ = SAME.initialize_u0(N_individuals,p_dist)#creates initial conditions, Matrix{BigFloat} of size NxN
t_span = (0.,t_max)
prob = ODEProblem(voter_model_ame!,u₀,t_span,p)
sol = solve(prob,Tsit5(),saveat=1)

The code take ~3 seconds for N = 50 and upwards of 10min to run with N = 500.
Here are the @time return for N = 50

4.603619 seconds (77.25 M allocations: 3.787 GiB, 8.55% gc time, 46.87% compilation time)

Now this seems to be very slow and takes up an insane amount of memory.

I used the in-place style function instead of allocating a new array every time step as per the optimization instructions in the DifferentialEquations.jl I can’t use static arrays because the size of my arrays is too large.

I have a feeling I am doing something horrendously wrong but I don’t know what. There must be an efficient way to solve large systems of differential equations. What am I missing?

Can you @time or benchmark your rhs! function?

du₀ =zeros(size(u₀))#50x50 Matrix{BigFloat}
@time SAME.voter_model_ame!(du₀,u₀,1,[0])


0.001388 seconds (24.50 k allocations: 1.215 MiB)

For 10 times steps do we expect 4 orders of magnitude slowdown? It seems like it should scale linearly in the number of time steps.

Secondly, it seems like it is taking a huge amount of memory. Does ODEProblem require the use of BigFloat, arbitrary precision floats?

1 Like

No, and in fact it should be much slower than using fixed-precision floats. Is there a specific reason you are generating u₀ with BigFloat? Using e.g. Float64 should get rid of many of the allocations in the ODE function (ideally you want this to have no allocations at all).


Okay wow that makes a crazy difference.

Now when I run the !rhs function I get

0.000009 seconds (1 allocation: 64 bytes)

So the BigFloat is what was slowing it down.

Thanks for the help

1 Like

Nice, I’m glad it helped!

I’m sure that there are some important use cases for arbitrary precision floats, but from my experience, Float64 should be the go-to in virtually all situations (or if memory is scarce perhaps Float32).