2

我想要解决一个复杂的网络系统，该系统涉及通过多维数组起作用的高阶交互项。我已经编写了相应的代码，但是获得结果需要花费太多时间。有什么方法可以更快地求解微分方程吗？代码如下：

  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")

#= 我发现的问题是在我的代码中计算广义订单参数需要太多时间。我是朱莉娅的新手。谁能向我建议一种更快的方法来计算该术语/以某种方式使代码更快，以便我可以在更短的时间内解决问题？=#

4

2 回答 2

3

sol = solve(prob, RK4(), reltol=1e-4, saveat=dts, progress=true)仅在这一行中似乎有很多问题。对于求解器来说，RK4 几乎总是一个糟糕的选择。这不是一种有效的方法。为什么不选择推荐的东西，或者使用自动算法选择？Tsit5在 99.99% 的情况下会好半个数量级，因此如果您正在寻找性能，请遵循建议。

此外，saveat这样的密集dts将比仅允许自适应时间步长选择保存点并稍后依赖插值要慢得多。这里的大部分时间不是用来解决而是保存解决方案，因此使用更简单的保存表示可以为您节省大量时间。

此外，您在微分方程定义中分配了很多数组。如果您正在寻找性能，不建议这样做。

所有这些问题都在优化 DiffEq 代码教程中进行了讨论，我强烈建议您在继续之前阅读该教程。

于 2021-12-25T10:17:28.237 回答
2

我怀疑你现在最大的瓶颈，至少在关于计算 的部分generalized_local_order，是你的矩阵A2[:, :, p]是密集的。您最好使用稀疏矩阵的向量，例如A2 = [spzeros(N, N) for i in 1:N]，然后使用A2[p].

于 2021-12-15T00:44:55.440 回答