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? =#

Can you fix the formatting? This is really hard to read currently — Oscar Smith, Dec 14 '21 at 20:57
There's a lot of global variables here too. Code that uses non-`const` global variables can't be optimized well by the compiler, that's why this is explained at the very top of the Performance Tips page of the Julia docs. It's possible to "patch" the code up with `const` and type annotations to allow some optimization, but refactoring as local variables in methods is better and far easier. That said I didn't profile your code so I can't say how much it would help, your bottleneck could very well be something else. — BatWannaBe, Dec 14 '21 at 22:59

score 3 · Answer 1 · answered Dec 25 '21 at 10:17

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.

Cardano · Answer 2 · 2021-12-15T15:05:08.503

2

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].

edited Dec 15 '21 at 15:05

answered Dec 15 '21 at 00:44

Cardano

931
1
8
14

Also, `N` appears to be a non-constant global variable – Oscar Smith Dec 15 '21 at 05:37
Thanks for your suggestions. But is it possible to define a multidimensional sparse array in julia?? As far as I know Sparsearray package is only for matrices not for higher dimensional arrays( tensors). – MSA Dec 15 '21 at 14:04
1

@MSA you're right, I didn't realize that. You can still use a vector of sparse matrices for your specific problem. – Cardano Dec 15 '21 at 15:05

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

2 Answers2