Hi all!
I was advised to open this issue as a follow-up of this thread: https://discourse.julialang.org/t/speeding-up-solution-to-the-large-system-of-non-linear-complex-valued-odes/39075
My issue: I am generating nonlinear complex-valued ODES from the infinite series, by essentially summing out the series for a given N (number of terms). Series have linear, quadratic and cubic terms in it (single, double and triple sums), making equations (summations) very long for large N as the number of terms in the entire system grows approximately as a quadratic function of N. Therefore when I call function containing my equations in the solver, code allocates large amounts of memory, making it impossible for me to solve for realistic N (I would be looking at N=100 as a minimum).
Examples of generated equations, which I give to the solver:
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
I tried putting parenthesis around every 4-th term or so to avoid splatting of +() but this did not solve the issue.
Implementation:
Currently, I generate equations as string arrays, write them to a text file, and then load them to the solver. Example code I use for cubic term:
using OffsetArrays, ..dnDot_String
function C(N)
Cubic=Vector{String}(undef,6N+1)
oCubic=OffsetArray(Cubic,-3N:3N)
dummy=Array{String}(undef,2N+1,2N+1,6N+1)
odummy=OffsetArray(dummy,-N:N,-N:N,-3N:3N)
fun=Vector{String}(undef,3)
cTuple=NTuple{3,Int64}
for q = -3N:1:3N #q, n and m are indices appearing in triple summation
for n = -N:1:N
for m = -N:1:N
if abs(q-n+m) <= abs(N) #condition on indices to select only the terms I need
cTuple = (n,m,q-n+m)
fun = dnDot_String.boolMatch(cTuple,N)
@inbounds odummy[n,m,q] = fun[1]*"*conj("*fun[2]*")*"*fun[3]*"+" #storing term into the array.
elseif abs(q-n+m) > abs(N)
@inbounds odummy[n,m,q] = ""
end
end
end
@inbounds oCubic[q] = chop(string(odummy1[:,:,q]...))::SubString # concantenating the array for each q. This gives vector where q-th component contains all the cubic terms for q-th equation
end
and function dnDot_String.boolMatch
function boolMatch(n::NTuple{T},N::Int64) where T<:Any
x=lastindex(n)
dummy=Vector{String}(undef,x)
for i=1:x
dummy[i]="u[$(n[i]+N+1)]"
end
return dummy
end
This essentially generates u[i] terms recognised by DifferentialEquations package base on the values of indices in cTuple.
So my question is, is there a more efficient way of generating equations in ModellingToolkit or otherwise that would avoid memory overhead I get when trying to solve the system?
Thank you!