# Solving coupled ODEs

I lately faced some problems trying to solve systems of coupled nonlinear ODEs with NDSolve. Cause I couldn't find a solution on my own, here I am with a MWE. Consider a system of $n$ coupled linear ODEs with coupling radius r.

$$\dot{x}_i=\frac{a}{2r}\sum_{j=i-r}^{i+r}(x_j-x_i)$$

where $x_i\in \mathbb{R}$, $i\in \{1,\dots,n\}$. Here the sum is used in a modulo kind of way. Every $j<1$, so $j=0,-1,-2,\dots$ corresponds to $j=n,\,n-1,\,n-2,\dots$ and every $j>n$, so $j=n+1,n+2,n+3\dots$ corresponds to $j=1,2,3\dots$.

For somewhat large $n$ and $r$ NDSolve fails to solve the system displaying:

Cannot solve to find an explicit formula for the derivatives.
Consider using the option Method->{"EquationSimplification"->"Residual"}.


I would like to know (if possible)

1. why NDSolve isn't able to find an explicit formula for the derivative and tells me to use an DAE-solver when the system clearly isn't DAE.
2. if there is some way to work around this.

Here is the code that produces the error

n = 1000; r = 400;
a = 0.05; tend = 500.;

vars = x[#] & /@ Range[1, n];

eqns = a/(2 r) Sum[x[Mod[j, n, 1]][t] - x[#][t], {j, # - r, # + r}]&/@Range[1, n];

ics = x[#][0] & /@ Range[1, n] == RandomReal[{-2.0, 2.0}, n];

sol = NDSolve[{x[#]'[t] & /@ Range[1, n] == eqns, ics}, vars, {t, 0, tend}][[1]];


I'm using version 10.2

• Can't reproduce the warning in v9.0.1 and v11.2… anyway, your first question has been explained in this post: mathematica.stackexchange.com/a/158519/1871 – xzczd Jun 27 '18 at 11:28
• @freddy90 Additionally, you could also look at the Fermi-Pasta-Ulam-Tsingou model on Mathematica demonstrations. What you describe here is a first order coupled set of point bodies (?). FPUT have a more complex model (with periodic boundary conditions, I think) that describes n bodies connected by quadratic springs and their evolution. Best wishes. – dearN Jun 27 '18 at 13:23
• @xzczd thanks to your guidance I was able to solve my problem by setting Method->{"EquationSimplification" -> "SolveExplicitly"} – freddy90 Jun 27 '18 at 16:18
• Oh, how did you figure out this option value? (I failed to find it in the document and this site. ) I suggest writing a self-answer to elaborate. – xzczd Jun 27 '18 at 16:35

By chance I found the option value

Method -> {"EquationSimplification" -> "SolveExplicitly"}


which works for me. Another possible option instead of "SolveExplicitly" would be "Solve". I also found out by trial and error that for highly coupled ODEs it might help to use the function Expand before feeding the equations to NDSolve

You have to use Thread for the definition of initial conditions and equations.

With Version 8.0, I get a result after about 30 minutes without problems.

n = 1000; r = 400;
a = 5/100; tend = 500;

vars = x[#] & /@ Range[1, n];

eqns = a/(2 r) Sum[x[Mod[j, n, 1]][t] - x[#][t], {j, # - r, # + r}] & /@
Range[1, n];

ics = Thread[x[#][0] & /@ Range[1, n] == RandomReal[{-2.0, 2.0}, n]];

sol = NDSolve[Join[Thread[x[#]'[t] & /@ Range[1, n] == eqns], ics],
vars, {t, 0, tend}][[1]];

Plot[Evaluate[Table[x[i][t] /. sol, {i, 900, 1000}]], {t, 0, 100},
PlotRange -> All, ImageSize -> 600]


• Yeah, in v8.0.4 equations like {x'[t], y'[t]}={…, …} is not supported by NDSolve, one needs to use e.g. Thread to transform the equation. The syntax in OP's code is supported since v9. – xzczd Jun 27 '18 at 16:43
• BTW, 30 minutes?? – xzczd Jun 27 '18 at 16:43