Hi,
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] -
1/2*a[m][t]*a[-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] =
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_] :=
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,
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}] +
1/2*U2[n];
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,
Nick
