Is there any faster method to solve differential equations with multidimensional array

Question

I want a solve a complex network system involving higher-order interaction terms acting through a multidimensional array. I have written the corresponding code but it takes too much time to get my results. Is there any possible way to solve the differential equation faster? The code is given below:

  using DelimitedFiles 
     
  using LightGraphs

  using LinearAlgebra

  using PyPlot

   using Random

   using BenchmarkTools


    N=500; #number of nodes

    A=readdlm("sf_simplicial.txt"); # graph adjacency called from a data

   G=Graph(A);

   global deg=degree(G);

   global omega=deg;

    mean(deg);


   A2=zeros(N,N,N); # `adjacency tensor which accounts the connection between three nodes

   for i in 1:N

for j in 1:N

    for k in 1:N

        if (A[i,j]==1 && A[j,k]==1 && A[k,i]==1)

            A2[i,j,k]=1

            A2[i,k,j]=1

            A2[j,k,i]=1

            A2[j,i,k]=1

            A2[k,i,j]=1

            A2[k,j,i]=1

        end;

    end;

end;

end;


   K= sum(p -> A2[:,:,p], 1:N)

   deg_sim= sum(j -> K[:,j], 1:N)/2;


   function Kuramoto(du, u, pp, t)

u1 = @view u[1:N] #### θ
du1 = @view du[1:N]  #### dθ
u2 = @view u[N+1:2*N]  ####### λ
du2 = @view du[N+1:2*N]  ####### dλ

  
α1=0.08
β1=0.04
σ1=1.0
σ2=1.0
λ0=pp

####### local_order
z1 = Array{Complex{Float64},1}(undef, N)
mul!(z1, A, exp.((u1)im))
z1 = z1 ./ deg 

####### generalized_local_order
sum1= sum(p -> A2[:,:,p] * exp(u1[p]im), 1:N)
z2 = Array{Complex{Float64},1}(undef, N)
mul!(z2, sum1, exp.((u1)im))
z2 = z2 ./ deg_sim


####### equ of motion
@. du1  = omega + u2 *( σ1 * deg * imag(z1 * exp((-1im) * u1)) + σ2 * deg_sim * imag(z2 * exp((-1im) * 2*u1)))
@. du2 = α1 *(λ0-u2)- β1 * (abs(z1)+ abs(z2))/2.0
    
  end;

 # setting up time steps and integration intervals

  dt = 0.01 # time step
  dts = 0.1 # save time
   ti = 0.0
   tt = 1000.0 
   tf = 10000.0
   nt = Int(div(tt,dts))
   nf = Int(div(tf,dts))

   tspan = (ti, tf); # time interval

   pp=0.65

   u0=[rand(N)*2*pi;pp*ones(N)];


    using DifferentialEquations

    prob = ODEProblem(kuramoto,u0 , tspan, pp)

    # solving problem using a proper solver
    sol = solve(prob, RK4(), reltol=1e-4, saveat=dts, progress=true);

     r1 = abs.(mean(exp.(sol[1:N,:]*1im),dims=1)[1,:])[nt:nf];

    t=range(0,stop=tf-tt,length=nf-nt+1)

    plot(t,r1)

    ylabel("R(t)")

    xlabel("t")

#= The problem I have found out is the calculation of generalized order parameter in my code it takes too much time. I am new to Julia. Can anyone suggest to me a faster method to calculate that term/ make the code faster somehow so that I can solve the problem in lesser time? =#

Answer 1

I suspect the biggest bottleneck for you right now, at least in the section about computing the generalized_local_order, is that your matrices A2[:, :, p] are dense. You will do better to use a vector of sparse matrices, like A2 = [spzeros(N, N) for i in 1:N], then access them with A2[p].

Answer 2

sol = solve(prob, RK4(), reltol=1e-4, saveat=dts, progress=true) There seems to be a lot going wrong just in this line. RK4 is almost always a bad choice for a solver. It's not an efficient method. Why not choose something that is recommended, or use the automatic algorithm choice? Tsit5 would be better in 99.99% of cases by like half an order of magnitude, so please follow the recommendations if you're looking for performance.

Additionally, saveat with such a dense dts will be much slower than just allowing the adaptive time stepping to choose the save points and relying later on the interpolation. A good chunk of the time here is not spent solving but saving solutions, so using a simpler representation of the save will save you a lot of time.

Additionally, you're allocating a lot of arrays inside of your differential equation definition. This is not recommended if you're looking for performance.

All of these issues are discussed in the Optimizing DiffEq Code tutorial which I would highly recommend you read before continuing.