5

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
freddy90
  • 851
  • 4
  • 12
  • 4
    Can't reproduce the warning in v9.0.1 and _v11.2_… anyway, your first question has been explained in this post: https://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
  • 1
    @xzczd Using something like Method -> {"EquationSimplification" -> nonsense} returns an error message listing allowed options, {Automatic, Residual, MassMatrix, Solve, SolveExplicitly}. Of course, it does not tell what the options do. – bbgodfrey Nov 09 '21 at 05:05

2 Answers2

2

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

freddy90
  • 851
  • 4
  • 12
0

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]

enter image description here

Akku14
  • 17,287
  • 14
  • 32
  • 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