# NDSolve for a large system of coupled (complex valued) ODEs

The system under consideration here is a generalization of the one I discussed in this problem. Qualitatively, I am trying to solve a system of coupled nonlinear (up to order 4 in the dep variables) ODEs which can take on complex values.

The governing equation for the $n-th$ ODE is of the form:

$\sum_{k=-M}^M (Q_{nk}+Q_{kn}) \ddot{a}_k -\sum_{\ell=-M}^M\sum_{k=-M}^M S_n(\ell,k) \dot{a}_{\ell} \dot{a}_k + U_n = 0$ where

$Q_{nk}=\frac{1}{2}\sum_{j=1}^MP_{jn}P_{(-j)k}$,

with

$P_{mn}=\frac{\sqrt{|m|}}{|n|}(a_{m-n}-\frac{a_ma_{-n}}{2})$,

(and the contsraints that $a_{m-n}=0$ for $|m-n|>M$ and $a_{m-n}=0$ when $Sign(m-n)=Sign(n)$),

$S_n(\ell,k)=\left(\frac{\partial Q_{k\ell}}{\partial a_n}-\frac{\partial Q_{n\ell}}{\partial a_k}-\frac{\partial Q_{\ell n}}{\partial a_{k}}\right)$ and

$U_n = \frac{\partial V}{\partial a_n}$ where

$V = \frac{1}{8}\left(\sum_{n=1}^M\frac{a_na_{-n}}{n}\right)^2 +\frac{1}{2}\sum_{n=1}^M\frac{a_na_{-n}}{n^2} -\frac{1}{2} \left(\sum_{n=1}^M\frac{a_na_{-n}}{n})\right)\left(\sum_{n=1}^M\frac{a_na_{-n}}{n}\right) + \sum_{i+j+k=0,(i,j,k)\neq 0}\frac{a_ia_ja_k}{4ij}$

and $M$ is the highest mode under consideration.

In the limit of small $a_i$, these governing equations reduce to a system of harmonic oscillators of the form $\ddot{a}_i=-ia_i$, so I've been exploring the solutions with initial conditions that should reproduce this and give solutions of the form $a_n(t)= e^{int}$.

When I set this system up in Mathematica, I can find coherent solutions for the $M=1,M=2$ cases. However, for $M=3$ the integration takes a very long time (so far more than 6 hours) so something is awry. The bottleneck seems to be occurring at the NDSolve portion of the code, as I think I've executed all HoldAll's that could be causing time issues.

For the simpler system linked above (where all coefficients are real) I was able to integrate the $M=50$ cases in about 1 hour, so I'm trying to see why NDSolve would be taking longer in this case. I note that in that case the system was "treated as a system of differential-algebraic equations". The equations in this example are not treated as DAE, since this complicates finding appropriate ICs and I experience a "mconly" error. This could be one reason for large discrepancies in run time.

Finally, there is a relationship between the coefficients such that $a_{-k} =a^*_k$ which I have not exploited and could potentially cut the number of dependent variables in half.

Any insights into getting larger M cases to integrate, or why this is taking so long in its current state, would be greatly appreciated.

Thanks,

Nick

M = 3; (*Number of modes*)
cond1[r_, s_] := If[Sign[s] == Sign[r - s], 0, 1];
cond2[r_, s_] := If[Abs[r - s] > M, 0, 1];
F[n_] := If[n == 0, 1, Abs[n]];
a[0][t_] := 1;
P[m_, k_, t_] :=
Sqrt[F[m]]/
F[k]*(cond1[m, k]*cond2[m, k]*a[m - k][t] -
1/2*a[m][t]*a[-k][t]) ;
Q[k_, l_, t_] = 1/4*Sum[P[m, k, t]*P[-m, l, t], {m, 1, M}];
q[m_, k_, n_, t_] :=
Sqrt[F[m]]/
F[k]*(cond1[m, k]*cond2[m, k]* Boole[m - k == n] -
a[m][t]/2*Boole[-k == n] - a[-k][t]/2*Boole[m == n]);
SQ[k_, l_, n_, t_] :=
1/4*Sum[P[-m, l, t]*q[m, k, n, t] + P[m, k, t]*q[-m, l, n, t], {m,
1, M}];
SS[k_, l_, n_, t_] := (SQ[k, l, n, t] - SQ[n, l, k, t] - SQ[l, n, k, t]);
Fu[a_, b_, c_] := If[a + b + c == 0, 1, 0];
Fu2[a_] := If[a == 0, 0, 1];
V[t_] := Evaluate[
1/2*(Sum[-a[n][t]*a[-n][t]/(2*n), {n, 1, M}]^2 +
2*Sum[a[n][t] a[-n][t]/(4*n^2), {n, 1, M}] +
Sum[-a[n][t]*a[-n][t]/(2 n), {n, 1, M}]*
Sum[a[n][t] a[-n][t]/n, {n, 1, M}] +
Sum[Fu[n, m, o]*Fu2[n]*Fu2[m]*a[n][t]*a[m][t]*a[o][t]*
Abs[o]/(4 Abs[F[n]*F[m]*F[o]]), {n, -M, M}, {m, -M,
M}, {o, -M, M}])];
U[n_, t_] := Evaluate[D[V[t], a[n][t]]];
Gov[n_, t_] :=
Sum[a[l]''[t]*(Q[n, l, t] + Q[l, n, t]), {l, -M, M}] -
Sum[SS[k, l, n, t] a[k]'[t] a[l]'[t], {k, -M, M}, {l, -M, M}] +
U[n, t];
ao = 0.01;
eqns2 = {
Evaluate[Gov[-3, t] == 0],
Evaluate[Gov[-2, t] == 0],
Evaluate[Gov[-1, t] == 0],
Evaluate[Gov[1, t] == 0],
Evaluate[Gov[2, t] == 0],
Evaluate[Gov[3, t] == 0],
a[1][0] == a[-1][0] == ao,
a[1]'[0] == I*ao, a[-1]'[0] == -I*ao,
a[2][0] == a[-2][0] == ao^2,
a[2]'[0] == ao^2*I*2,
a[-2]'[0] == -ao^2*I*2 ,
a[3][0] == a[-3][0] == ao^3,
a[3]'[0] == 3*I*ao^3,
a[-3]'[0]==-3*I*ao^3
};
MMode = NDSolve[
eqns2, {a[-3][t],a[-2][t], a[-1][t], a[1][t], a[2][t], a[3][t]}, {t, 0, 10}, Method ->      {"EquationSimplification" -> "Solve"}]
• there is a syntax error or something else. copying/pasting your code gives errors : !Mathematica graphics Commented Aug 19, 2013 at 2:19
• Hi Nasser, there were a few commas missing when I was copy/pasting to the forum. Thanks for pointing this out. Commented Aug 19, 2013 at 5:01
• @NickP, unfortunately there is no IDA code for complex values, this is why the message suggests that should use Solve as an equation simplification method, that however, takes a long time to finish. It is calling Solve to modify your system of equations. I do not know of a workaround for this. You could send this to the support and ask for a recommendation.
– user21
Commented Aug 19, 2013 at 9:06
• @ruebenko thanks for the reply. Perhaps I can split the governing equations into real and imaginary parts and see if that helps out. Commented Aug 19, 2013 at 19:20
• I have indeed split the equations into real and imaginary parts, allowing the IDA code to run and solve these equations in reasonable time. Commented Aug 23, 2013 at 22:15