# Tips for efficiently solving large system coupled (nonlinear) ODEs

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 == 0.3, B == 0.1, Table[B[j] ==   0, {j, 3, M}], Table[B[j]' == 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.

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[FrontEndEvaluationNotebook[], 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 == 0.3, B == 0.1, Table[B[j] == 0, {j, 3, M}], Table[B[j]' == 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. • ruebenko, thanks for the suggestions. They are appreciated. I understand the rationale (and I think the implementation) of evaluating the integral for$F[t]\$ explicitly before integrating. However, I note that I have been running the code for over an hour (on a macbook pro) and so far it has not yielded a solution for the M=10 case. Aug 5, 2013 at 21:15
• @NickP, once the Integrate is done (about 30 sec.) NDSolve finishes almost immediately. What version of M- do you have?
• @NickP, it is impossible to say what exactly might be the issue. It could be an improved algorithm or fixing of a bug. And yes, I use 9.0.1. I just checked V8 and it seems slower. For M=4 it takes 21s. vs. 17s on v9 (different machines though), for M=5` 27s. vs. 18s.