Message Boards Message Boards

Efficiently tabulating system of coupled ODEs for NDSolve

Posted 11 years ago

I am solving a (large) system of second order coupled nonlinear ODEs.

I have been increasing the number of dependent variables present in the system slowly, to converge on more physically realistic results, but have experienced some pretty severe bottlenecks in the code.

I will link to some theory (is LaTex supported in this forum?) and present my code below, but my question is, what is the best (i.e. most efficient) way to tabulate a large system of indexed functions, to be passed into NDSolve?

I've tried a few things to improve efficiency which I'll list after I present the code, but first some details.

The system under consideration is described here. Qualitatively it is a system of coupled second order nonlinear ODEs. The i-th governing equation is described in terms of a series of summations, and is of the form Sum_j (Q_ij +Q_ji) d^2a_i/dt^2 + Sum_jSum_k S_i (j,k) da_j/dt da_k/dt + V = 0, where the coefficients (Q_ij, S_i (j,k), V) are of maximum degree 4 polynomials in the dependent variables a_i.
The particular formulation of these terms is given in the link but in general I am dealing with a number of indexed functions, which leads to a bit of my confusion, and will be discussed below the code.

The code is
 M = 10;
 cond1[r_, s_] :=
   cond1[r, s] =
    If[Sign[s] == Sign[r - s], 0, 1];
 cond2[r_, s_] :=
   cond2[r, s] =
    If[Abs[r - s] > M, 0, 1];
 F[n_] := F[n] = If[n == 0, 1, Abs[n]];
 a[0][t_] := a[0][t] = 1.;
P[m_, k_, t_] :=
  P[m, k, t] =
   1./(Sqrt[F[m]]*F[k])*(cond1[m, k]*cond2[m, k]*a[m - k][t] -
Q[k_, l_, t_] :=
  Q[k, l, t] = 1./4.*Sum[P[m, k, t]*P[-m, l, t], {m, 1, M}];
q[m_, k_, n_, t_] :=
  q[m, k, n, t] =
   1./(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_] :=
  SQ[k, l, n, t] =
     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_] :=
  SS[k, l, n,
    t] = (SQ[k, l, n, t] - SQ[n, l, k, t] - SQ[l, n, k, t]);
Fu[a_, b_, c_] := Fu[a, b, c] = If[a + b + c == 0, 1, 0];
Fu2[a_] := Fu2[a] = If[a == 0, 0, 1];
V[t_] := V[
    t] = (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]/(8*Abs[F[n]*F[m]*F[o]]), {n, -M, M}, {m, -M, M}, {o, -M,

Join[Table[With[{i = i}, U2[i] = D[V[t], a[i][t]]], {i, 1, M}],
  Table[With[{i = i}, U2[-i] = D[V[t], a[-i][t]]], {i, 1, M}]];
Gov[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}] +

Table[With[{i = i}, a[i][t_] := ar[i][t] + I*ai[i][t]], {i, 1, M}];
Table[With[{i = i}, a[-i][t_] := ar[i][t] - I*ai[i][t]], {i, 1, M}];
Table[ar[i][t] \[Element] Reals, {i, 1, M}];
Table[ai[i][t] \[Element] Reals, {i, 1, M}];

Table[With[{i = i}, Test[i] = Evaluate[Gov[i, t]]], {i, 1, M}];
Table[With[{i = i}, ReGov[i, t_] := ComplexExpand[Re[Test[i]]]], {i,
   1, M}];
Table[With[{i = i}, ImGov[i, t_] := ComplexExpand[Im[Test[i]]]], {i,
   1, M}];

ao = 0.1;

eqns = Join[Parallelize[Table[ReGov[i, t] == 0, {i, 1, M}]],
   Parallelize[Table[ImGov[i, t] == 0, {i, 1, M}]],
   Table[ar[i][0] == ao, {i, 1}], Table[ar[i][0] == 0, {i, 2, M}],
   Table[ai[i][0] == 0, {i, 1, M}], Table[ar[i]'[0] == 0, {i, 1, M}],
   Table[ai[i]'[0] == ao, {i, 1}], Table[ai[i]'[0] == 0, {i, 2, M}]];

MMode = NDSolve[eqns,
   Join[Table[ar[i][t], {i, 1, M}], Table[ai[i][t], {i, 1, M}]], {t,
    0, 10}];

The first bottleneck occurs when defining a variable that evaluates the governing equations:
Table[ With[{i = i}, Test = Evaluate[Gov[i, t]]], {i, 1, M}];

Note, for the M=10 case the Timing on my machine is about 1.55, while for the M=20 case, the Timing is about 28.555. In particular, the bottleneck seems to come from evaluating the function SS[n,k,l,t] in the definition of Gov[i,t]. There are quite a few terms in this function, but they all involve low order polynomials of maximum degree 4 in the dependent variables. I'm not sure if there is a way to speed these calculations up, or if this is a necessary area of intense computation.

Next, I decompose my governing equations into real and imaginary parts, to speed up the NDSolve portion of the computation (if I don't, I get an IDA error and must use "Solve" - see here). This leads to the next, more severe, bottle neck, which comes from tabulating the real and imaginary parts of the governing equations, i.e. defining the variable "eqns". I have used ComplexExpand[Re(Im) ] to get the real and imaginary parts of the governing equations, respectively, but this seems to be computationally expensive and accounts for the gross majority of the computation time. The reason for this post is to find out if there is a more efficient way to do this. Note, for the M=5 case the Timing of the eqns calculation is 0.0173 while for the M=10 case it is 12.25, i.e. 4 orders of magnitude larger!

Qualitatively, I wonder if I have not effectively constrained some of the definitions in my code, and if this generality is leading to significant slow down. Also, I feel as if I need to gain some understanding of the compile function, which seems to be very useful when defining functions, but the subtleties of using it with indexed functions to be passed into NSDSolve is beyond my very naive understanding and I cannot find any relevant examples to study. Finally, an alternative approach would be using matrices instead of term by term operations, but I'm not convinced this will lead to a significant speed up in Mathematica, but would welcome any advice on this. Any comments/suggestions are greatly appreciated,

POSTED BY: Nick Pizzo
2 Replies
Posted 11 years ago
Hi Peter,

Thanks for pointing this out. I had some trouble formatting my post, and in the confusion did not realize many terms had been dropped. I believe the code is now as it should be.

POSTED BY: Nick Pizzo
Posted 11 years ago
I'm noting that the expression Sign == Sign[r - s] in cond1 doesn't make sense, according to your reference it's a typo (Sign -> Sign[ s ] ).

Is all the remaining code as desired?
POSTED BY: Peter Fleck
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract