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




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



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_] := 
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_] := 
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,
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"}]
  • $\begingroup$ there is a syntax error or something else. copying/pasting your code gives errors : !Mathematica graphics $\endgroup$
    – Nasser
    Commented Aug 19, 2013 at 2:19
  • $\begingroup$ Hi Nasser, there were a few commas missing when I was copy/pasting to the forum. Thanks for pointing this out. $\endgroup$
    – Nick P
    Commented Aug 19, 2013 at 5:01
  • $\begingroup$ @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. $\endgroup$
    – user21
    Commented Aug 19, 2013 at 9:06
  • $\begingroup$ @ruebenko thanks for the reply. Perhaps I can split the governing equations into real and imaginary parts and see if that helps out. $\endgroup$
    – Nick P
    Commented Aug 19, 2013 at 19:20
  • $\begingroup$ I have indeed split the equations into real and imaginary parts, allowing the IDA code to run and solve these equations in reasonable time. $\endgroup$
    – Nick P
    Commented Aug 23, 2013 at 22:15


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.