Tips for efficiently solving large system coupled (nonlinear) ODEs

Question

I'm trying to solve a system of nonlinear, coupled ODEs, where the governing equation for the $n$-th ODE is of the form:

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

where

$$Q_{nk}=\frac{1}{2}\sum_{j=1}^MP_{jn}P_{jk},$$

with

$$P_{jn}=\frac{1}{\sqrt{m}n}(a_{n-j}+a_{j+n}-a_ja_n),\quad \mathrm{and} \quad a_{n-j}=0 \quad \mathrm{for}\quad j>n,$$

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

and

$$U_n=\frac{\partial}{\partial a_n} \frac{1}{2\pi} \int_0^{2\pi}\left\{ \left(-\sum_{j=1}^M \left(\frac{a_j^2}{2j}-\frac{a_j}{j}\cos j\theta \right)\right)^2\,\left(1+\sum_{j=1}^M a_j \cos j\theta \right) \right\}d\theta.$$

Note, $(Q_{k\ell}, S_N(k,\ell), U_n)$ are all polynomials (of maximum order 4) in the dependent variables $a_i(t)$ and $M$ is the total number of dependent variables under consideration (i.e. $a_j\equiv 0$ for $j>M$). I have been able to solve for low $M$ cases (i.e. $M$ less than 10) using NDSolve. For fixed initial conditions, I have been increasing $M$ and finding a significant increase in computation time. I'm wondering if there are any obvious ways to decrease computation time that I have overlooked. (Ideally, for the system I'm trying to model, I'd like to increase the number of dependent variables to around 200.)

I'm new to Mathematica, and I'm sure I'm making some 'rookie mistakes' that are costing me in computation time. Also, if I'm barking up the wrong tree, that is I should be using a different platform, I'd like to find out before I invest too much effort into this project.

Some questions I have on basic problem set-up:

1. Is it more efficient to isolate the dependent variables in question. That is, have equations of the form $\ddot{a}_i = \ldots$ versus $\sum_k^M Q_{nk} \ddot{a}_k = \ldots$?

2. Is it worth fully simplifying the governing equations (i.e. using FullSimplify) before executing the NDSolve command?

3. Should I rewrite each second order ODEs as a pair of first order ODEs?

Any advice on this would be greatly appreciated.

Here is the code:

M = 10; (*Number of modes*)
f[r_, s_] := If[r > s, 0, 1]; (*Will ensure we have no negative modes*)
α[i_, t_] := If[i <= M, A[i, t], 0]; (*Truncate series at Mth mode*)
A[0, t_] := 1; (*By definition*)

P[m_, n_, t_] := 
1/Sqrt[m*n^2]*(f[m, n]*α[n - m, t] + α[n + m, 
t] - α[m, t]*α[n, t]) ;  

Q[n_, l_, t_] := 1/2 Sum[P[m, n, t]* P[m, l, t], {m,1,M}] ; 

S[n_, k_, l_, t_] := 
1/2 (D[Q[k, l, t], A[n, t]] - D[Q[n, l, t], A[k, t]] - 
D[Q[n, k, t], A[l, t]]); 

F[t_] := 1/(2 π) Integrate[
Sum[((A[j, t]*Cos[j*ξ])/j - 
       A[j, t]^2/(2j)),{j,1,M}]^2* 
(1+Sum[1+A[j, t]*
      Cos[j*ξ],{j,1,M}]), {ξ, 0, 2π}]; 
U[n_, t_] := D[F[t], A[n, t]];
Table[A[j, t_] = B[j][t], {j, 1, M}]; (*Redefine dep variables so that independent     variable is distinguished from dep var index*)
Gov[n_, t_] := Sum[B[l]''[t]*Q[n,l,t],{l,1,M}]-Sum[Sum[S[n,k,l,t]B[k]'[t]B[l]'[t],      {k,1,M}],{l,1,M}]+1/2U[n,t];

eqns = {Table[Gov[j, t] == 0, {j, 1, M}], B[1][0] == 0.3, B[2][0] == 0.1, Table[B[j][0] ==   0, {j, 3, M}], Table[B[j]'[0] == 0, {j, 1, M}]};

TenMode = NDSolve[eqns, Table[B[n][t], {n, 1, M}], {t, 0, 10}]

I apologize for any errors in the code since I had some issues parsing from the symbolic notation present in my editor to a code that would be easily copied into this forum.

Answer 1

The bottleneck in this example is not NDSolve but the fact that for every time you call F a symbolic integration needs to be done. F is defined via SetDelayed (:=) and that function has Attributes HoldAll. What you can do is to force a symbolic (!) evaluation of the integral and use that evaluated integral in stead.

showStatus[status_] := 
  LinkWrite[$ParentLink, 
   SetNotebookStatusLine[FrontEnd`EvaluationNotebook[], 
    ToString[status]]];
clearStatus[] := showStatus[""];
clearStatus[]

M = 10;(*Number of modes*)

f[r_, s_] := 
 If[r > s, 0, 1];(*Will ensure we have no negative modes*)α[i_,
   t_] := If[i <= M, A[i, t], 0];(*Truncate series at Mth mode*)

A[0, t_] := 1;(*By definition*)
P[m_, n_, t_] := 
 1/Sqrt[m*n^2]*(f[m, n]*α[n - m, t] + α[n + m, 
     t] - α[m, t]*α[n, t]);
Q[n_, l_, t_] := 1/2 Sum[P[m, n, t]*P[m, l, t], {m, 1, M}];
S[n_, k_, l_, t_] := 
  1/2 (D[Q[k, l, t], A[n, t]] - D[Q[n, l, t], A[k, t]] - 
     D[Q[n, k, t], A[l, t]]);

(* pre-evaluate the integral and define F to be that *)
F[t_] := Evaluate[
   1/(2 π) Integrate[
      Sum[((A[j, t]*Cos[j*ξ])/j - A[j, t]^2/(2 j)), {j, 1, M}]^2*
       Sum[1 + A[j, t]*Cos[j*ξ]/j, {j, 1, M}], {ξ, 0, 
       2 π}] // Simplify];

U[n_, t_] := D[F[t], A[n, t]];
Table[A[j, t_] = B[j][t], {j, 1, M}];

Gov[n_, t_] := 
  Evaluate[Sum[B[l]''[t]*Q[n, l, t], {l, 1, M}] - 
    Sum[Sum[S[n, k, l, t] B[k]'[t] B[l]'[t], {k, 1, M}], {l, 1, M}] + 
    1/2 U[n, t]];
eqns = {Table[Gov[j, t] == 0, {j, 1, M}], B[1][0] == 0.3, 
   B[2][0] == 0.1, Table[B[j][0] == 0, {j, 3, M}], 
   Table[B[j]'[0] == 0, {j, 1, M}]};
TenMode = 
 NDSolve[eqns, Table[B[n][t], {n, 1, M}], {t, 0, 10}, 
  EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]
  (*,Method->{"EquationSimplification"->"Residual"}*)
  ]

This gives a message (NDSolve::ndsdtc). If you comment in the (*Method->..*) this will go away. Also have a look at the documentation.

Concerning your questions:

1. Not really, NDSolve does that for you.
2. Not really, the overhead in this example is not the numerical function evaluation. In some cases using e.g. Simplify will make the equations smaller and thus less function evaluations are necessary.
3. No, NDSolve does that for you.

Hope this helps.