Speeding up solution to the large system of non-linear, complex-valued ODEs

General Usage
, , ,
1

First of all, thank you to the people who developed Julia and DifferentialEquations.jl these both are really great tools!

Now, cutting to the chase:

I need to solve a large system of non-linear complex-valued ODEs (N~1000). For now, I am testing my code on a small system, and currently am battling with the memory allocation issue, which I can not improve despite following (hopefully, correctly) performance improvement tutorials. MWE:

module Dif
using DifferentialEquations, JLD2, StaticArrays

const tspan=(0.0,100.0)
const t_name=tspan[2]
function Nodes(u,p,t) #equations to be solved. For large N, they get very long and ugly
A1=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[1]-p[2]/2*(p[1])^2*conj(u[7]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[4])+conj(p[1])*u[1]*u[4]+2*p[1]*u[2]*conj(u[5])+conj(p[1])*u[2]*u[3]+2*p[1]*u[3]*conj(u[6])+conj(p[1])*u[3]*u[2]+2*p[1]*u[4]*conj(u[7])+conj(p[1])*u[4]*u[1] + u[1]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[2]+u[2]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[3]+u[2]*conj(u[3])*u[2]+u[3]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[4]+u[2]*conj(u[4])*u[3]+u[3]*conj(u[4])*u[2]+u[4]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[5]+u[2]*conj(u[5])*u[4]+u[3]*conj(u[5])*u[3]+u[4]*conj(u[5])*u[2]+u[5]*conj(u[5])*u[1]+u[1]*conj(u[6])*u[6]+u[2]*conj(u[6])*u[5]+u[3]*conj(u[6])*u[4]+u[4]*conj(u[6])*u[3]+u[5]*conj(u[6])*u[2]+u[6]*conj(u[6])*u[1]+u[1]*conj(u[7])*u[7]+u[2]*conj(u[7])*u[6]+u[3]*conj(u[7])*u[5]+u[4]*conj(u[7])*u[4]+u[5]*conj(u[7])*u[3]+u[6]*conj(u[7])*u[2]+u[7]*conj(u[7])*u[1]))/im
A2=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[2]-p[2]/2*(p[1])^2*conj(u[6]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[3])+conj(p[1])*u[1]*u[5]+2*p[1]*u[2]*conj(u[4])+conj(p[1])*u[2]*u[4]+2*p[1]*u[3]*conj(u[5])+conj(p[1])*u[3]*u[3]+2*p[1]*u[4]*conj(u[6])+conj(p[1])*u[4]*u[2]+2*p[1]*u[5]*conj(u[7])+conj(p[1])*u[5]*u[1] + u[1]*conj(u[1])*u[2]+u[2]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[3]+u[2]*conj(u[2])*u[2]+u[3]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[4]+u[2]*conj(u[3])*u[3]+u[3]*conj(u[3])*u[2]+u[4]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[5]+u[2]*conj(u[4])*u[4]+u[3]*conj(u[4])*u[3]+u[4]*conj(u[4])*u[2]+u[5]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[6]+u[2]*conj(u[5])*u[5]+u[3]*conj(u[5])*u[4]+u[4]*conj(u[5])*u[3]+u[5]*conj(u[5])*u[2]+u[6]*conj(u[5])*u[1]+u[1]*conj(u[6])*u[7]+u[2]*conj(u[6])*u[6]+u[3]*conj(u[6])*u[5]+u[4]*conj(u[6])*u[4]+u[5]*conj(u[6])*u[3]+u[6]*conj(u[6])*u[2]+u[7]*conj(u[6])*u[1]+u[2]*conj(u[7])*u[7]+u[3]*conj(u[7])*u[6]+u[4]*conj(u[7])*u[5]+u[5]*conj(u[7])*u[4]+u[6]*conj(u[7])*u[3]+u[7]*conj(u[7])*u[2]))/im
A3=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[3]-p[2]/2*(p[1])^2*conj(u[5]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[2])+conj(p[1])*u[1]*u[6]+2*p[1]*u[2]*conj(u[3])+conj(p[1])*u[2]*u[5]+2*p[1]*u[3]*conj(u[4])+conj(p[1])*u[3]*u[4]+2*p[1]*u[4]*conj(u[5])+conj(p[1])*u[4]*u[3]+2*p[1]*u[5]*conj(u[6])+conj(p[1])*u[5]*u[2]+2*p[1]*u[6]*conj(u[7])+conj(p[1])*u[6]*u[1] + u[1]*conj(u[1])*u[3]+u[2]*conj(u[1])*u[2]+u[3]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[4]+u[2]*conj(u[2])*u[3]+u[3]*conj(u[2])*u[2]+u[4]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[5]+u[2]*conj(u[3])*u[4]+u[3]*conj(u[3])*u[3]+u[4]*conj(u[3])*u[2]+u[5]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[6]+u[2]*conj(u[4])*u[5]+u[3]*conj(u[4])*u[4]+u[4]*conj(u[4])*u[3]+u[5]*conj(u[4])*u[2]+u[6]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[7]+u[2]*conj(u[5])*u[6]+u[3]*conj(u[5])*u[5]+u[4]*conj(u[5])*u[4]+u[5]*conj(u[5])*u[3]+u[6]*conj(u[5])*u[2]+u[7]*conj(u[5])*u[1]+u[2]*conj(u[6])*u[7]+u[3]*conj(u[6])*u[6]+u[4]*conj(u[6])*u[5]+u[5]*conj(u[6])*u[4]+u[6]*conj(u[6])*u[3]+u[7]*conj(u[6])*u[2]+u[3]*conj(u[7])*u[7]+u[4]*conj(u[7])*u[6]+u[5]*conj(u[7])*u[5]+u[6]*conj(u[7])*u[4]+u[7]*conj(u[7])*u[3]))/im
A4=((-p[2]/2*abs2(p[1]))*u[4]-p[2]/2*(p[1])^2*conj(u[4]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[1])+conj(p[1])*u[1]*u[7]+2*p[1]*u[2]*conj(u[2])+conj(p[1])*u[2]*u[6]+2*p[1]*u[3]*conj(u[3])+conj(p[1])*u[3]*u[5]+2*p[1]*u[4]*conj(u[4])+conj(p[1])*u[4]*u[4]+2*p[1]*u[5]*conj(u[5])+conj(p[1])*u[5]*u[3]+2*p[1]*u[6]*conj(u[6])+conj(p[1])*u[6]*u[2]+2*p[1]*u[7]*conj(u[7])+conj(p[1])*u[7]*u[1] + u[1]*conj(u[1])*u[4]+u[2]*conj(u[1])*u[3]+u[3]*conj(u[1])*u[2]+u[4]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[5]+u[2]*conj(u[2])*u[4]+u[3]*conj(u[2])*u[3]+u[4]*conj(u[2])*u[2]+u[5]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[6]+u[2]*conj(u[3])*u[5]+u[3]*conj(u[3])*u[4]+u[4]*conj(u[3])*u[3]+u[5]*conj(u[3])*u[2]+u[6]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[7]+u[2]*conj(u[4])*u[6]+u[3]*conj(u[4])*u[5]+u[4]*conj(u[4])*u[4]+u[5]*conj(u[4])*u[3]+u[6]*conj(u[4])*u[2]+u[7]*conj(u[4])*u[1]+u[2]*conj(u[5])*u[7]+u[3]*conj(u[5])*u[6]+u[4]*conj(u[5])*u[5]+u[5]*conj(u[5])*u[4]+u[6]*conj(u[5])*u[3]+u[7]*conj(u[5])*u[2]+u[3]*conj(u[6])*u[7]+u[4]*conj(u[6])*u[6]+u[5]*conj(u[6])*u[5]+u[6]*conj(u[6])*u[4]+u[7]*conj(u[6])*u[3]+u[4]*conj(u[7])*u[7]+u[5]*conj(u[7])*u[6]+u[6]*conj(u[7])*u[5]+u[7]*conj(u[7])*u[4]))/im
A5=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[5]-p[2]/2*(p[1])^2*conj(u[3]) - (p[2]/2)*(2*p[1]*u[2]*conj(u[1])+conj(p[1])*u[2]*u[7]+2*p[1]*u[3]*conj(u[2])+conj(p[1])*u[3]*u[6]+2*p[1]*u[4]*conj(u[3])+conj(p[1])*u[4]*u[5]+2*p[1]*u[5]*conj(u[4])+conj(p[1])*u[5]*u[4]+2*p[1]*u[6]*conj(u[5])+conj(p[1])*u[6]*u[3]+2*p[1]*u[7]*conj(u[6])+conj(p[1])*u[7]*u[2] + u[1]*conj(u[1])*u[5]+u[2]*conj(u[1])*u[4]+u[3]*conj(u[1])*u[3]+u[4]*conj(u[1])*u[2]+u[5]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[6]+u[2]*conj(u[2])*u[5]+u[3]*conj(u[2])*u[4]+u[4]*conj(u[2])*u[3]+u[5]*conj(u[2])*u[2]+u[6]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[7]+u[2]*conj(u[3])*u[6]+u[3]*conj(u[3])*u[5]+u[4]*conj(u[3])*u[4]+u[5]*conj(u[3])*u[3]+u[6]*conj(u[3])*u[2]+u[7]*conj(u[3])*u[1]+u[2]*conj(u[4])*u[7]+u[3]*conj(u[4])*u[6]+u[4]*conj(u[4])*u[5]+u[5]*conj(u[4])*u[4]+u[6]*conj(u[4])*u[3]+u[7]*conj(u[4])*u[2]+u[3]*conj(u[5])*u[7]+u[4]*conj(u[5])*u[6]+u[5]*conj(u[5])*u[5]+u[6]*conj(u[5])*u[4]+u[7]*conj(u[5])*u[3]+u[4]*conj(u[6])*u[7]+u[5]*conj(u[6])*u[6]+u[6]*conj(u[6])*u[5]+u[7]*conj(u[6])*u[4]+u[5]*conj(u[7])*u[7]+u[6]*conj(u[7])*u[6]+u[7]*conj(u[7])*u[5]))/im
A6=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[6]-p[2]/2*(p[1])^2*conj(u[2]) - (p[2]/2)*(2*p[1]*u[3]*conj(u[1])+conj(p[1])*u[3]*u[7]+2*p[1]*u[4]*conj(u[2])+conj(p[1])*u[4]*u[6]+2*p[1]*u[5]*conj(u[3])+conj(p[1])*u[5]*u[5]+2*p[1]*u[6]*conj(u[4])+conj(p[1])*u[6]*u[4]+2*p[1]*u[7]*conj(u[5])+conj(p[1])*u[7]*u[3] + u[1]*conj(u[1])*u[6]+u[2]*conj(u[1])*u[5]+u[3]*conj(u[1])*u[4]+u[4]*conj(u[1])*u[3]+u[5]*conj(u[1])*u[2]+u[6]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[7]+u[2]*conj(u[2])*u[6]+u[3]*conj(u[2])*u[5]+u[4]*conj(u[2])*u[4]+u[5]*conj(u[2])*u[3]+u[6]*conj(u[2])*u[2]+u[7]*conj(u[2])*u[1]+u[2]*conj(u[3])*u[7]+u[3]*conj(u[3])*u[6]+u[4]*conj(u[3])*u[5]+u[5]*conj(u[3])*u[4]+u[6]*conj(u[3])*u[3]+u[7]*conj(u[3])*u[2]+u[3]*conj(u[4])*u[7]+u[4]*conj(u[4])*u[6]+u[5]*conj(u[4])*u[5]+u[6]*conj(u[4])*u[4]+u[7]*conj(u[4])*u[3]+u[4]*conj(u[5])*u[7]+u[5]*conj(u[5])*u[6]+u[6]*conj(u[5])*u[5]+u[7]*conj(u[5])*u[4]+u[5]*conj(u[6])*u[7]+u[6]*conj(u[6])*u[6]+u[7]*conj(u[6])*u[5]+u[6]*conj(u[7])*u[7]+u[7]*conj(u[7])*u[6]))/im
A7=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[7]-p[2]/2*(p[1])^2*conj(u[1]) - (p[2]/2)*(2*p[1]*u[4]*conj(u[1])+conj(p[1])*u[4]*u[7]+2*p[1]*u[5]*conj(u[2])+conj(p[1])*u[5]*u[6]+2*p[1]*u[6]*conj(u[3])+conj(p[1])*u[6]*u[5]+2*p[1]*u[7]*conj(u[4])+conj(p[1])*u[7]*u[4] + u[1]*conj(u[1])*u[7]+u[2]*conj(u[1])*u[6]+u[3]*conj(u[1])*u[5]+u[4]*conj(u[1])*u[4]+u[5]*conj(u[1])*u[3]+u[6]*conj(u[1])*u[2]+u[7]*conj(u[1])*u[1]+u[2]*conj(u[2])*u[7]+u[3]*conj(u[2])*u[6]+u[4]*conj(u[2])*u[5]+u[5]*conj(u[2])*u[4]+u[6]*conj(u[2])*u[3]+u[7]*conj(u[2])*u[2]+u[3]*conj(u[3])*u[7]+u[4]*conj(u[3])*u[6]+u[5]*conj(u[3])*u[5]+u[6]*conj(u[3])*u[4]+u[7]*conj(u[3])*u[3]+u[4]*conj(u[4])*u[7]+u[5]*conj(u[4])*u[6]+u[6]*conj(u[4])*u[5]+u[7]*conj(u[4])*u[4]+u[5]*conj(u[5])*u[7]+u[6]*conj(u[5])*u[6]+u[7]*conj(u[5])*u[5]+u[6]*conj(u[6])*u[7]+u[7]*conj(u[6])*u[6]+u[7]*conj(u[7])*u[7]))/im

@SVector [A1, A2, A3, A4, A5, A6, A7]
end

function Diff(N::Int64,psi::Float64,g::Float64)
println(tspan)

Analytic_N=Int64(floor((N-1)/2))

k=sqrt(g*psi) #defining parameters
p = [psi,g,k] 

u1=[0.0+0.0*im for i in 1:N]  # simple initial conditions
u1[Int64((N-1)/2)]=1.0+0.0*im    
u1[Int64((N-1)/2)+1]=0.0+0.0*im      
u1[Int64((N-1)/2)+2]=1.0+0.0*im  
u0=SVector{N}(u1)

path=pwd()
path1="/path/JLD2 for g=$g, psi=$psi, N=$Analytic_N,t=$t_name"
rm(path1,force=true,recursive=true)
mkdir(path1)
cd(path1)

prob = ODEProblem(Nodes,u0,tspan,p)
sol = solve(prob,RK4(), adaptive=true,dense=false,calck=false)
@save "sol_tf=$(tspan[2]).jld2" sol
cd(path)
return "success"
end
end

Benchmarking this gives:

@benchmark Dif.Diff($7,$1.0,$16.0)                                                         
BenchmarkTools.Trial:                                   
  memory estimate:  2.53 GiB                            
  allocs estimate:  58532496                            
  --------------                                        
  minimum time:     1.857 s (11.35% GC)                 
  median time:      1.867 s (11.93% GC)                 
  mean time:        1.876 s (11.73% GC)                 
  maximum time:     1.903 s (11.90% GC)                 
  --------------                                        
  samples:          3                                   
  evals/sample:     1

Is there anything else I can do to improve the memory allocation? Perhaps I could use a better-suited solver? Write down equations in matrix representation (I would expect this would slow the program down)?
When I run the above for the system of 13 equations, memory allocation jump to 26.06GB. I have also tried using saveat=0.01 but found that it increases memory allocation for adaptive=true, not by much, but somewhat. If I run this code for system of 41 equations, I get out of memory error.
My supervisor solves very similar equations using RK4 on C++, and is able to solve system of 1000 equations without any issues in around 30 minutes, which makes me conclude that I am doing something wrong.

Any advice is greatly appreciated!

Thank you.

2

Put parenthesis around each 4th term of so. There’s 36 terms nested deep in * which causes an issue. Parse 1 + 2 + 3 as +(1, +(2, 3)) · Issue #34999 · JuliaLang/julia · GitHub is the proposed fix.

For now, an automated way to work around this is to use ModelingToolkit to regenerate your function:

using ModelingToolkit, StaticArrays
@variables t u[1:7](t) p[1:3]
@derivatives D'~t
function Nodes(u,p,t) #equations to be solved. For large N, they get very long and ugly
  A1=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[1]-p[2]/2*(p[1])^2*conj(u[7]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[4])+conj(p[1])*u[1]*u[4]+2*p[1]*u[2]*conj(u[5])+conj(p[1])*u[2]*u[3]+2*p[1]*u[3]*conj(u[6])+conj(p[1])*u[3]*u[2]+2*p[1]*u[4]*conj(u[7])+conj(p[1])*u[4]*u[1] + u[1]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[2]+u[2]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[3]+u[2]*conj(u[3])*u[2]+u[3]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[4]+u[2]*conj(u[4])*u[3]+u[3]*conj(u[4])*u[2]+u[4]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[5]+u[2]*conj(u[5])*u[4]+u[3]*conj(u[5])*u[3]+u[4]*conj(u[5])*u[2]+u[5]*conj(u[5])*u[1]+u[1]*conj(u[6])*u[6]+u[2]*conj(u[6])*u[5]+u[3]*conj(u[6])*u[4]+u[4]*conj(u[6])*u[3]+u[5]*conj(u[6])*u[2]+u[6]*conj(u[6])*u[1]+u[1]*conj(u[7])*u[7]+u[2]*conj(u[7])*u[6]+u[3]*conj(u[7])*u[5]+u[4]*conj(u[7])*u[4]+u[5]*conj(u[7])*u[3]+u[6]*conj(u[7])*u[2]+u[7]*conj(u[7])*u[1]))/im
  A2=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[2]-p[2]/2*(p[1])^2*conj(u[6]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[3])+conj(p[1])*u[1]*u[5]+2*p[1]*u[2]*conj(u[4])+conj(p[1])*u[2]*u[4]+2*p[1]*u[3]*conj(u[5])+conj(p[1])*u[3]*u[3]+2*p[1]*u[4]*conj(u[6])+conj(p[1])*u[4]*u[2]+2*p[1]*u[5]*conj(u[7])+conj(p[1])*u[5]*u[1] + u[1]*conj(u[1])*u[2]+u[2]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[3]+u[2]*conj(u[2])*u[2]+u[3]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[4]+u[2]*conj(u[3])*u[3]+u[3]*conj(u[3])*u[2]+u[4]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[5]+u[2]*conj(u[4])*u[4]+u[3]*conj(u[4])*u[3]+u[4]*conj(u[4])*u[2]+u[5]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[6]+u[2]*conj(u[5])*u[5]+u[3]*conj(u[5])*u[4]+u[4]*conj(u[5])*u[3]+u[5]*conj(u[5])*u[2]+u[6]*conj(u[5])*u[1]+u[1]*conj(u[6])*u[7]+u[2]*conj(u[6])*u[6]+u[3]*conj(u[6])*u[5]+u[4]*conj(u[6])*u[4]+u[5]*conj(u[6])*u[3]+u[6]*conj(u[6])*u[2]+u[7]*conj(u[6])*u[1]+u[2]*conj(u[7])*u[7]+u[3]*conj(u[7])*u[6]+u[4]*conj(u[7])*u[5]+u[5]*conj(u[7])*u[4]+u[6]*conj(u[7])*u[3]+u[7]*conj(u[7])*u[2]))/im
  A3=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[3]-p[2]/2*(p[1])^2*conj(u[5]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[2])+conj(p[1])*u[1]*u[6]+2*p[1]*u[2]*conj(u[3])+conj(p[1])*u[2]*u[5]+2*p[1]*u[3]*conj(u[4])+conj(p[1])*u[3]*u[4]+2*p[1]*u[4]*conj(u[5])+conj(p[1])*u[4]*u[3]+2*p[1]*u[5]*conj(u[6])+conj(p[1])*u[5]*u[2]+2*p[1]*u[6]*conj(u[7])+conj(p[1])*u[6]*u[1] + u[1]*conj(u[1])*u[3]+u[2]*conj(u[1])*u[2]+u[3]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[4]+u[2]*conj(u[2])*u[3]+u[3]*conj(u[2])*u[2]+u[4]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[5]+u[2]*conj(u[3])*u[4]+u[3]*conj(u[3])*u[3]+u[4]*conj(u[3])*u[2]+u[5]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[6]+u[2]*conj(u[4])*u[5]+u[3]*conj(u[4])*u[4]+u[4]*conj(u[4])*u[3]+u[5]*conj(u[4])*u[2]+u[6]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[7]+u[2]*conj(u[5])*u[6]+u[3]*conj(u[5])*u[5]+u[4]*conj(u[5])*u[4]+u[5]*conj(u[5])*u[3]+u[6]*conj(u[5])*u[2]+u[7]*conj(u[5])*u[1]+u[2]*conj(u[6])*u[7]+u[3]*conj(u[6])*u[6]+u[4]*conj(u[6])*u[5]+u[5]*conj(u[6])*u[4]+u[6]*conj(u[6])*u[3]+u[7]*conj(u[6])*u[2]+u[3]*conj(u[7])*u[7]+u[4]*conj(u[7])*u[6]+u[5]*conj(u[7])*u[5]+u[6]*conj(u[7])*u[4]+u[7]*conj(u[7])*u[3]))/im
  A4=((-p[2]/2*abs2(p[1]))*u[4]-p[2]/2*(p[1])^2*conj(u[4]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[1])+conj(p[1])*u[1]*u[7]+2*p[1]*u[2]*conj(u[2])+conj(p[1])*u[2]*u[6]+2*p[1]*u[3]*conj(u[3])+conj(p[1])*u[3]*u[5]+2*p[1]*u[4]*conj(u[4])+conj(p[1])*u[4]*u[4]+2*p[1]*u[5]*conj(u[5])+conj(p[1])*u[5]*u[3]+2*p[1]*u[6]*conj(u[6])+conj(p[1])*u[6]*u[2]+2*p[1]*u[7]*conj(u[7])+conj(p[1])*u[7]*u[1] + u[1]*conj(u[1])*u[4]+u[2]*conj(u[1])*u[3]+u[3]*conj(u[1])*u[2]+u[4]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[5]+u[2]*conj(u[2])*u[4]+u[3]*conj(u[2])*u[3]+u[4]*conj(u[2])*u[2]+u[5]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[6]+u[2]*conj(u[3])*u[5]+u[3]*conj(u[3])*u[4]+u[4]*conj(u[3])*u[3]+u[5]*conj(u[3])*u[2]+u[6]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[7]+u[2]*conj(u[4])*u[6]+u[3]*conj(u[4])*u[5]+u[4]*conj(u[4])*u[4]+u[5]*conj(u[4])*u[3]+u[6]*conj(u[4])*u[2]+u[7]*conj(u[4])*u[1]+u[2]*conj(u[5])*u[7]+u[3]*conj(u[5])*u[6]+u[4]*conj(u[5])*u[5]+u[5]*conj(u[5])*u[4]+u[6]*conj(u[5])*u[3]+u[7]*conj(u[5])*u[2]+u[3]*conj(u[6])*u[7]+u[4]*conj(u[6])*u[6]+u[5]*conj(u[6])*u[5]+u[6]*conj(u[6])*u[4]+u[7]*conj(u[6])*u[3]+u[4]*conj(u[7])*u[7]+u[5]*conj(u[7])*u[6]+u[6]*conj(u[7])*u[5]+u[7]*conj(u[7])*u[4]))/im
  A5=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[5]-p[2]/2*(p[1])^2*conj(u[3]) - (p[2]/2)*(2*p[1]*u[2]*conj(u[1])+conj(p[1])*u[2]*u[7]+2*p[1]*u[3]*conj(u[2])+conj(p[1])*u[3]*u[6]+2*p[1]*u[4]*conj(u[3])+conj(p[1])*u[4]*u[5]+2*p[1]*u[5]*conj(u[4])+conj(p[1])*u[5]*u[4]+2*p[1]*u[6]*conj(u[5])+conj(p[1])*u[6]*u[3]+2*p[1]*u[7]*conj(u[6])+conj(p[1])*u[7]*u[2] + u[1]*conj(u[1])*u[5]+u[2]*conj(u[1])*u[4]+u[3]*conj(u[1])*u[3]+u[4]*conj(u[1])*u[2]+u[5]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[6]+u[2]*conj(u[2])*u[5]+u[3]*conj(u[2])*u[4]+u[4]*conj(u[2])*u[3]+u[5]*conj(u[2])*u[2]+u[6]*conj(u[2])*u[1]+u[1]*conj(u[3])*u[7]+u[2]*conj(u[3])*u[6]+u[3]*conj(u[3])*u[5]+u[4]*conj(u[3])*u[4]+u[5]*conj(u[3])*u[3]+u[6]*conj(u[3])*u[2]+u[7]*conj(u[3])*u[1]+u[2]*conj(u[4])*u[7]+u[3]*conj(u[4])*u[6]+u[4]*conj(u[4])*u[5]+u[5]*conj(u[4])*u[4]+u[6]*conj(u[4])*u[3]+u[7]*conj(u[4])*u[2]+u[3]*conj(u[5])*u[7]+u[4]*conj(u[5])*u[6]+u[5]*conj(u[5])*u[5]+u[6]*conj(u[5])*u[4]+u[7]*conj(u[5])*u[3]+u[4]*conj(u[6])*u[7]+u[5]*conj(u[6])*u[6]+u[6]*conj(u[6])*u[5]+u[7]*conj(u[6])*u[4]+u[5]*conj(u[7])*u[7]+u[6]*conj(u[7])*u[6]+u[7]*conj(u[7])*u[5]))/im
  A6=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[6]-p[2]/2*(p[1])^2*conj(u[2]) - (p[2]/2)*(2*p[1]*u[3]*conj(u[1])+conj(p[1])*u[3]*u[7]+2*p[1]*u[4]*conj(u[2])+conj(p[1])*u[4]*u[6]+2*p[1]*u[5]*conj(u[3])+conj(p[1])*u[5]*u[5]+2*p[1]*u[6]*conj(u[4])+conj(p[1])*u[6]*u[4]+2*p[1]*u[7]*conj(u[5])+conj(p[1])*u[7]*u[3] + u[1]*conj(u[1])*u[6]+u[2]*conj(u[1])*u[5]+u[3]*conj(u[1])*u[4]+u[4]*conj(u[1])*u[3]+u[5]*conj(u[1])*u[2]+u[6]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[7]+u[2]*conj(u[2])*u[6]+u[3]*conj(u[2])*u[5]+u[4]*conj(u[2])*u[4]+u[5]*conj(u[2])*u[3]+u[6]*conj(u[2])*u[2]+u[7]*conj(u[2])*u[1]+u[2]*conj(u[3])*u[7]+u[3]*conj(u[3])*u[6]+u[4]*conj(u[3])*u[5]+u[5]*conj(u[3])*u[4]+u[6]*conj(u[3])*u[3]+u[7]*conj(u[3])*u[2]+u[3]*conj(u[4])*u[7]+u[4]*conj(u[4])*u[6]+u[5]*conj(u[4])*u[5]+u[6]*conj(u[4])*u[4]+u[7]*conj(u[4])*u[3]+u[4]*conj(u[5])*u[7]+u[5]*conj(u[5])*u[6]+u[6]*conj(u[5])*u[5]+u[7]*conj(u[5])*u[4]+u[5]*conj(u[6])*u[7]+u[6]*conj(u[6])*u[6]+u[7]*conj(u[6])*u[5]+u[6]*conj(u[7])*u[7]+u[7]*conj(u[7])*u[6]))/im
  A7=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[7]-p[2]/2*(p[1])^2*conj(u[1]) - (p[2]/2)*(2*p[1]*u[4]*conj(u[1])+conj(p[1])*u[4]*u[7]+2*p[1]*u[5]*conj(u[2])+conj(p[1])*u[5]*u[6]+2*p[1]*u[6]*conj(u[3])+conj(p[1])*u[6]*u[5]+2*p[1]*u[7]*conj(u[4])+conj(p[1])*u[7]*u[4] + u[1]*conj(u[1])*u[7]+u[2]*conj(u[1])*u[6]+u[3]*conj(u[1])*u[5]+u[4]*conj(u[1])*u[4]+u[5]*conj(u[1])*u[3]+u[6]*conj(u[1])*u[2]+u[7]*conj(u[1])*u[1]+u[2]*conj(u[2])*u[7]+u[3]*conj(u[2])*u[6]+u[4]*conj(u[2])*u[5]+u[5]*conj(u[2])*u[4]+u[6]*conj(u[2])*u[3]+u[7]*conj(u[2])*u[2]+u[3]*conj(u[3])*u[7]+u[4]*conj(u[3])*u[6]+u[5]*conj(u[3])*u[5]+u[6]*conj(u[3])*u[4]+u[7]*conj(u[3])*u[3]+u[4]*conj(u[4])*u[7]+u[5]*conj(u[4])*u[6]+u[6]*conj(u[4])*u[5]+u[7]*conj(u[4])*u[4]+u[5]*conj(u[5])*u[7]+u[6]*conj(u[5])*u[6]+u[7]*conj(u[5])*u[5]+u[6]*conj(u[6])*u[7]+u[7]*conj(u[6])*u[6]+u[7]*conj(u[7])*u[7]))/im
  @SVector [A1, A2, A3, A4, A5, A6, A7]
end
rhs = Nodes(u,p,t)
eqs = @. D(u) ~ rhs
sys = ODESystem(Array(eqs),t,u,p)
Nodes2 = generate_function(sys,expression=Val{false})[1]

N = 7; g = 0.1; psi = 0.2
Analytic_N=Int64(floor((N-1)/2))

k=sqrt(g*psi) #defining parameters
p = [psi,g,k]

u1=[0.0+0.0*im for i in 1:N]  # simple initial conditions
u1[Int64((N-1)/2)]=1.0+0.0*im
u1[Int64((N-1)/2)+1]=0.0+0.0*im
u1[Int64((N-1)/2)+2]=1.0+0.0*im
u0=SVector{N}(u1)

using BenchmarkTools
@btime Nodes(u0,p,0.0)
@btime Nodes2(u0,p,0.0)

5.920 μs (313 allocations: 14.27 KiB)
578.571 ns (2 allocations: 272 bytes)

That’s a whole order of magnitude faster. FWIW, it seems like you might be programmatically generating these equations, and ModelingToolkit might be a nice way to do that which will also automatically parallelize if you use the parallelism options. This may be worth thinking about as the complexity grows.

5 Likes
3

Thank you for your response Chris! I do indeed generate equations programmatically.
I tried putting parenthesis around every term, this did not improve performance by the same order of magnitude as what your benchmarking is showing. Will placing them around every ~ 4th term instead make the desired difference?
Also, when trying to build your solution in my editor, I get the following error:

ERROR: LoadError: MethodError: no method matching _build_function(::ModelingToolkit.JuliaTarget, ::Array{Operation,1}, ::Array{Variable{Number},1}, ::Array{Variable{Number},1}, ::Tuple{Symbol}, ::ModelingToolkit.ODEToExpr, ::Type{Val{true}}; expression=Val{false})
Closest candidates are:
  _build_function(::ModelingToolkit.JuliaTarget, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; checkbounds, constructor, linenumbers, multithread) at C:\Users\visenkon\.julia\packages\ModelingToolkit\N4mo2\src\build_function.jl:51 got unsupported keyword argument "expression"
  _build_function(::ModelingToolkit.JuliaTarget, !Matched::Operation, ::Any, ::Any, ::Any, ::Any, ::Any; checkbounds, constructor, linenumbers) at C:\Users\visenkon\.julia\packages\ModelingToolkit\N4mo2\src\build_function.jl:16 got unsupported keyword argument "expression"
  _build_function(::ModelingToolkit.JuliaTarget, ::Any, ::Any, ::Any, ::Any, ::Any) at C:\Users\visenkon\.julia\packages\ModelingToolkit\N4mo2\src\build_function.jl:51 got unsupported keyword argument "expression"
  ...
Stacktrace:
 [1] kwerr(::NamedTuple{(:expression,),Tuple{DataType}}, ::Function, ::ModelingToolkit.JuliaTarget, ::Array{Operation,1}, ::Array{Variable{Number},1}, ::Array{Variable{Number},1}, ::Tuple{Symbol}, ::ModelingToolkit.ODEToExpr, ::Type) at .\error.jl:125
 [2] (::ModelingToolkit.var"#kw##_build_function")(::NamedTuple{(:expression,),Tuple{DataType}}, ::typeof(ModelingToolkit._build_function), ::ModelingToolkit.JuliaTarget, ::Array{Operation,1}, ::Array{Variable{Number},1}, ::Array{Variable{Number},1}, ::Tuple{Symbol}, ::ModelingToolkit.ODEToExpr, ::Type) at .\none:0
 [3] #build_function#156(::ModelingToolkit.JuliaTarget, ::Base.Iterators.Pairs{Symbol,DataType,Tuple{Symbol},NamedTuple{(:expression,),Tuple{DataType}}}, ::typeof(ModelingToolkit.build_function), ::Array{Operation,1}, ::Vararg{Any,N} where N) at C:\Users\visenkon\.julia\packages\ModelingToolkit\N4mo2\src\build_function.jl:8
 [4] (::ModelingToolkit.var"#kw##build_function")(::NamedTuple{(:expression,),Tuple{DataType}}, ::typeof(ModelingToolkit.build_function), ::Array{Operation,1}, ::Vararg{Any,N} where N) at .\none:0
 [5] #generate_function#73(::Base.Iterators.Pairs{Symbol,DataType,Tuple{Symbol},NamedTuple{(:expression,),Tuple{DataType}}}, ::typeof(generate_function), ::ODESystem, ::Array{Variable,1}, ::Array{Variable,1}, ::Type) at C:\Users\visenkon\.julia\packages\ModelingToolkit\N4mo2\src\systems\diffeqs\abstractodesystem.jl:51
 [6] (::ModelingToolkit.var"#kw##generate_function")(::NamedTuple{(:expression,),Tuple{DataType}}, ::typeof(generate_function), ::ODESystem, ::Array{Variable,1}, ::Array{Variable,1}, ::Type) at .\none:0 (repeats 2 times)
 [7] top-level scope at C:\Users\visenkon\Desktop\chris_sol.jl:17
 [8] include at .\boot.jl:328 [inlined]
 [9] include_relative(::Module, ::String) at .\loading.jl:1105
 [10] include(::Module, ::String) at .\Base.jl:31
 [11] exec_options(::Base.JLOptions) at .\client.jl:287
 [12] _start() at .\client.jl:460
in expression starting at C:\Users\visenkon\Desktop\chris_sol.jl:17
[Finished in 54.9s]

Thank you!

4

Are you on the latest ModelingToolkit?

5

Hi Chris,

Thank you for your reply!

  1. I tried putting parenthesis around (roughly) every 4-th term. I include equations:
A1=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[1]-p[2]/2*(p[1])^2*conj(u[7]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[7])+conj(p[1])*u[1]*u[1] + u[1]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[2]+u[2]*conj(u[2])*u[1]+(u[1]*conj(u[3])*u[3])+(u[2]*conj(u[3])*u[2])+u[3]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[4]+u[2]*conj(u[4])*u[3]+(u[3]*conj(u[4])*u[2])+u[4]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[5]+u[2]*conj(u[5])*u[4]+u[3]*conj(u[5])*u[3]+u[4]*conj(u[5])*u[2]+u[5]*conj(u[5])*u[1]+(u[1]*conj(u[6])*u[6])+(u[2]*conj(u[6])*u[5])+u[3]*conj(u[6])*u[4]+(u[4]*conj(u[6])*u[3])+(u[5]*conj(u[6])*u[2])+u[6]*conj(u[6])*u[1]+u[1]*conj(u[7])*u[7]+u[2]*conj(u[7])*u[6]+(u[3]*conj(u[7])*u[5])+u[4]*conj(u[7])*u[4]+u[5]*conj(u[7])*u[3]+(u[6]*conj(u[7])*u[2])+u[7]*conj(u[7])*u[1]))/im
A2=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[2]-p[2]/2*(p[1])^2*conj(u[6]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[6])+conj(p[1])*u[1]*u[2]+(2*p[1]*u[2]*conj(u[7])+conj(p[1])*u[2]*u[1]) + u[1]*conj(u[1])*u[2]+u[2]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[3]+u[2]*conj(u[2])*u[2]+u[3]*conj(u[2])*u[1]+(u[1]*conj(u[3])*u[4])+(u[2]*conj(u[3])*u[3])+(u[3]*conj(u[3])*u[2])+u[4]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[5]+u[2]*conj(u[4])*u[4]+u[3]*conj(u[4])*u[3]+(u[4]*conj(u[4])*u[2])+u[5]*conj(u[4])*u[1]+u[1]*conj(u[5])*u[6]+u[2]*conj(u[5])*u[5]+u[3]*conj(u[5])*u[4]+u[4]*conj(u[5])*u[3]+u[5]*conj(u[5])*u[2]+u[6]*conj(u[5])*u[1]+(u[1]*conj(u[6])*u[7])+(u[2]*conj(u[6])*u[6])+(u[3]*conj(u[6])*u[5])+u[4]*conj(u[6])*u[4]+(u[5]*conj(u[6])*u[3])+(u[6]*conj(u[6])*u[2])+u[7]*conj(u[6])*u[1]+u[2]*conj(u[7])*u[7]+u[3]*conj(u[7])*u[6]+(u[4]*conj(u[7])*u[5])+u[5]*conj(u[7])*u[4]+u[6]*conj(u[7])*u[3]+(u[7]*conj(u[7])*u[2])))/im
A3=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[3]-p[2]/2*(p[1])^2*conj(u[5]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[5])+conj(p[1])*u[1]*u[3]+(2*p[1]*u[2]*conj(u[6])+conj(p[1])*u[2]*u[2])+2*p[1]*u[3]*conj(u[7])+conj(p[1])*u[3]*u[1] + u[1]*conj(u[1])*u[3]+(u[2]*conj(u[1])*u[2])+(u[3]*conj(u[1])*u[1])+u[1]*conj(u[2])*u[4]+u[2]*conj(u[2])*u[3]+u[3]*conj(u[2])*u[2]+u[4]*conj(u[2])*u[1]+(u[1]*conj(u[3])*u[5])+u[2]*conj(u[3])*u[4]+u[3]*conj(u[3])*u[3]+u[4]*conj(u[3])*u[2]+u[5]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[6]+(u[2]*conj(u[4])*u[5])+(u[3]*conj(u[4])*u[4])+(u[4]*conj(u[4])*u[3])+u[5]*conj(u[4])*u[2]+(u[6]*conj(u[4])*u[1])+u[1]*conj(u[5])*u[7]+u[2]*conj(u[5])*u[6]+u[3]*conj(u[5])*u[5]+u[4]*conj(u[5])*u[4]+(u[5]*conj(u[5])*u[3])+u[6]*conj(u[5])*u[2]+u[7]*conj(u[5])*u[1]+u[2]*conj(u[6])*u[7]+u[3]*conj(u[6])*u[6]+u[4]*conj(u[6])*u[5]+u[5]*conj(u[6])*u[4]+u[6]*conj(u[6])*u[3]+u[7]*conj(u[6])*u[2]+(u[3]*conj(u[7])*u[7])+(u[4]*conj(u[7])*u[6])+u[5]*conj(u[7])*u[5]+(u[6]*conj(u[7])*u[4])+(u[7]*conj(u[7])*u[3])))/im
A4=((-p[2]/2*abs2(p[1]))*u[4]-p[2]/2*(p[1])^2*conj(u[4]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[4])+conj(p[1])*u[1]*u[4]+(2*p[1]*u[2]*conj(u[5])+conj(p[1])*u[2]*u[3])+2*p[1]*u[3]*conj(u[6])+conj(p[1])*u[3]*u[2]+(2*p[1]*u[4]*conj(u[7])+conj(p[1])*u[4]*u[1]) + u[1]*conj(u[1])*u[4]+u[2]*conj(u[1])*u[3]+u[3]*conj(u[1])*u[2]+u[4]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[5]+(u[2]*conj(u[2])*u[4])+u[3]*conj(u[2])*u[3]+u[4]*conj(u[2])*u[2]+u[5]*conj(u[2])*u[1]+(u[1]*conj(u[3])*u[6])+u[2]*conj(u[3])*u[5]+(u[3]*conj(u[3])*u[4])+(u[4]*conj(u[3])*u[3])+(u[5]*conj(u[3])*u[2])+u[6]*conj(u[3])*u[1]+u[1]*conj(u[4])*u[7]+u[2]*conj(u[4])*u[6]+u[3]*conj(u[4])*u[5]+u[4]*conj(u[4])*u[4]+u[5]*conj(u[4])*u[3]+(u[6]*conj(u[4])*u[2])+u[7]*conj(u[4])*u[1]+(u[2]*conj(u[5])*u[7])+u[3]*conj(u[5])*u[6]+u[4]*conj(u[5])*u[5]+u[5]*conj(u[5])*u[4]+u[6]*conj(u[5])*u[3]+u[7]*conj(u[5])*u[2]+(u[3]*conj(u[6])*u[7])+(u[4]*conj(u[6])*u[6])+(u[5]*conj(u[6])*u[5])+u[6]*conj(u[6])*u[4]+(u[7]*conj(u[6])*u[3])+u[4]*conj(u[7])*u[7]+u[5]*conj(u[7])*u[6]+(u[6]*conj(u[7])*u[5])+u[7]*conj(u[7])*u[4]))/im
A5=((1*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[5]-p[2]/2*(p[1])^2*conj(u[3]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[3])+conj(p[1])*u[1]*u[5]+(2*p[1]*u[2]*conj(u[4])+conj(p[1])*u[2]*u[4])+2*p[1]*u[3]*conj(u[5])+conj(p[1])*u[3]*u[3]+(2*p[1]*u[4]*conj(u[6])+conj(p[1])*u[4]*u[2])+2*p[1]*u[5]*conj(u[7])+conj(p[1])*u[5]*u[1] + u[1]*conj(u[1])*u[5]+u[2]*conj(u[1])*u[4]+u[3]*conj(u[1])*u[3]+u[4]*conj(u[1])*u[2]+u[5]*conj(u[1])*u[1]+u[1]*conj(u[2])*u[6]+u[2]*conj(u[2])*u[5]+(u[3]*conj(u[2])*u[4])+u[4]*conj(u[2])*u[3]+u[5]*conj(u[2])*u[2]+u[6]*conj(u[2])*u[1]+(u[1]*conj(u[3])*u[7])+(u[2]*conj(u[3])*u[6])+u[3]*conj(u[3])*u[5]+(u[4]*conj(u[3])*u[4])+(u[5]*conj(u[3])*u[3])+(u[6]*conj(u[3])*u[2])+u[7]*conj(u[3])*u[1]+u[2]*conj(u[4])*u[7]+u[3]*conj(u[4])*u[6]+u[4]*conj(u[4])*u[5]+u[5]*conj(u[4])*u[4]+u[6]*conj(u[4])*u[3]+(u[7]*conj(u[4])*u[2])+(u[3]*conj(u[5])*u[7])+u[4]*conj(u[5])*u[6]+u[5]*conj(u[5])*u[5]+u[6]*conj(u[5])*u[4]+u[7]*conj(u[5])*u[3]+(u[4]*conj(u[6])*u[7])+(u[5]*conj(u[6])*u[6])+(u[6]*conj(u[6])*u[5])+u[7]*conj(u[6])*u[4]+u[5]*conj(u[7])*u[7]+u[6]*conj(u[7])*u[6]+(u[7]*conj(u[7])*u[5])))/im
A6=((4*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[6]-p[2]/2*(p[1])^2*conj(u[2]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[2])+conj(p[1])*u[1]*u[6]+(2*p[1]*u[2]*conj(u[3])+conj(p[1])*u[2]*u[5])+2*p[1]*u[3]*conj(u[4])+conj(p[1])*u[3]*u[4]+(2*p[1]*u[4]*conj(u[5])+conj(p[1])*u[4]*u[3])+2*p[1]*u[5]*conj(u[6])+conj(p[1])*u[5]*u[2]+(2*p[1]*u[6]*conj(u[7])+conj(p[1])*u[6]*u[1]) + u[1]*conj(u[1])*u[6]+(u[2]*conj(u[1])*u[5])+(u[3]*conj(u[1])*u[4])+u[4]*conj(u[1])*u[3]+(u[5]*conj(u[1])*u[2])+(u[6]*conj(u[1])*u[1])+u[1]*conj(u[2])*u[7]+u[2]*conj(u[2])*u[6]+u[3]*conj(u[2])*u[5]+u[4]*conj(u[2])*u[4]+u[5]*conj(u[2])*u[3]+u[6]*conj(u[2])*u[2]+u[7]*conj(u[2])*u[1]+u[2]*conj(u[3])*u[7]+u[3]*conj(u[3])*u[6]+(u[4]*conj(u[3])*u[5])+u[5]*conj(u[3])*u[4]+u[6]*conj(u[3])*u[3]+u[7]*conj(u[3])*u[2]+(u[3]*conj(u[4])*u[7])+u[4]*conj(u[4])*u[6]+(u[5]*conj(u[4])*u[5])+(u[6]*conj(u[4])*u[4])+(u[7]*conj(u[4])*u[3])+u[4]*conj(u[5])*u[7]+u[5]*conj(u[5])*u[6]+u[6]*conj(u[5])*u[5]+u[7]*conj(u[5])*u[4]+u[5]*conj(u[6])*u[7]+u[6]*conj(u[6])*u[6]+u[7]*conj(u[6])*u[5]+(u[6]*conj(u[7])*u[7])+(u[7]*conj(u[7])*u[6])))/im
A7=((9*((p[3])^2/2)-p[2]/2*abs2(p[1]))*u[7]-p[2]/2*(p[1])^2*conj(u[1]) - (p[2]/2)*(2*p[1]*u[1]*conj(u[1])+conj(p[1])*u[1]*u[7]+(2*p[1]*u[2]*conj(u[2])+conj(p[1])*u[2]*u[6])+2*p[1]*u[3]*conj(u[3])+conj(p[1])*u[3]*u[5]+(2*p[1]*u[4]*conj(u[4])+conj(p[1])*u[4]*u[4])+2*p[1]*u[5]*conj(u[5])+conj(p[1])*u[5]*u[3]+(2*p[1]*u[6]*conj(u[6])+conj(p[1])*u[6]*u[2])+2*p[1]*u[7]*conj(u[7])+conj(p[1])*u[7]*u[1] + u[1]*conj(u[1])*u[7]+u[2]*conj(u[1])*u[6]+u[3]*conj(u[1])*u[5]+u[4]*conj(u[1])*u[4]+u[5]*conj(u[1])*u[3]+u[6]*conj(u[1])*u[2]+u[7]*conj(u[1])*u[1]+(u[2]*conj(u[2])*u[7])+u[3]*conj(u[2])*u[6]+u[4]*conj(u[2])*u[5]+(u[5]*conj(u[2])*u[4])+u[6]*conj(u[2])*u[3]+u[7]*conj(u[2])*u[2]+(u[3]*conj(u[3])*u[7])+(u[4]*conj(u[3])*u[6])+u[5]*conj(u[3])*u[5]+(u[6]*conj(u[3])*u[4])+(u[7]*conj(u[3])*u[3])+u[4]*conj(u[4])*u[7]+u[5]*conj(u[4])*u[6]+u[6]*conj(u[4])*u[5]+u[7]*conj(u[4])*u[4]+(u[5]*conj(u[5])*u[7])+u[6]*conj(u[5])*u[6]+u[7]*conj(u[5])*u[5]+(u[6]*conj(u[6])*u[7])+(u[7]*conj(u[6])*u[6])+u[7]*conj(u[7])*u[7]))/im

After benchmarking this, I found out that it allocates a little more memory with additional parenthesis and solver gives

┌ Warning: Interrupted. Larger maxiters is needed. 
└ @ DiffEqBase C:\Users\visenkon\.julia\packages\DiffEqBase\wBR51\src\integrator_interface.jl:329

which does not happen when I solve original equations.
2. I have updated ModellingToolkit and everything seems to be working, however when I take larger N, for example, N=81, it took 12 minutes in MT to build, and when I tried to load these equations in a separate module to solve, they just would not load.

Is there anything else I can try? Is there a way of programmatically generating my equations in MT?

Thank you!

6

You can raise the maxiters, or make RK4 non-adaptive and set a dt, but most likely you’ll want to use a method other than RK4 since it’s not very efficient.

Yes, if you define your variables as MT symbols and directly write the equations then it should work, but it sounds like you’re hitting some bad behavior. Could you open an issue with some example code of what you’re seeing? That way we can take a look.

1 Like
7

Hi Chris,

Thank you for your reply. I have created separate issue here https://discourse.julialang.org/t/generating-equations-in-modellingtoolkit-from-infinite-series/40244
Thank you!