Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Limitation on vector length in NDSolve?

• To: mathgroup at smc.vnet.net
• Subject: [mg122116] Limitation on vector length in NDSolve?
• From: DmitryG <einschlag at gmail.com>
• Date: Sat, 15 Oct 2011 06:02:38 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com

Dear All,

Working on vectorization of solutions of systems of ODEs, I have
stumbled upon a problem that is likely a limitation on the length of
the vectors that can be processed by NDSolve. Below a critical length
that is 845 in my case, the calculation is fast but above the critical
length it suddenly becomes very slow, although the results are
correct. Probably NDSolve switches to another and less efficient
method.

The task is to solve the system of equations

ds[t][[i]]/dt==-s[t][[i]]Sum[ f[i-j]s[t][[j]] ,{j,1,NN}],
i=1..NN,             f[i_]=Exp[-a i^2],          a=0.0001;

with some initial conditions

Below are two solutions, with vectorization (very slow!) and with
vectorization.


***************************************************************************
(* Definitions *)
tMax = 10.;
NN = 200;  (* Number of sites *)
a = 0.001;
f[i_] = Exp[-a i^2];
IniConds = s[0] == Table[1.  i/NN, {i, 1, NN}];

(* Solution with direct summation in the RHS - slow *)
Eqs = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(s[t]\)\) == -Table[
     s[t][[i]] Sum[f[i - j] s[t][[j]], {j, 1, NN}], {i, 1, NN}];
AbsoluteTiming[
 Sol = NDSolve[{Eqs, IniConds}, s, {t, 0, tMax},
   StartingStepSize -> 1/100,
   Method -> {"FixedStep", Method -> "ExplicitEuler"}]]
st[t_] := s[t] /. Sol[[1]];

ListPlot[{st[0], st[.003], st[0.01], st[0.05], st[0.5]}]


(* Vectorized solution *)
AbsoluteTiming[fMatr = Table[f[i - j], {i, 1, NN}, {j, 1, NN}];]
H[s_] := fMatr.s;
Eqs = \!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\(s[t]\)\) == -s[t]
H[s[t]];

AbsoluteTiming[
 Sol = NDSolve[{Eqs, IniConds}, s, {t, 0, tMax},
   StartingStepSize -> 1/100,
   Method -> {"FixedStep", Method -> "ExplicitEuler"}]]
st[t_] := s[t] /. Sol[[1]];
ListPlot[{st[0], st[.003], st[0.01], st[0.05], st[0.5]}]

**************************************************************************

I do not know if the problem depends on the computer because at the
moment I have Mathematica on my laptop only. It would be great if you
checked the issue. With NN around the critical value 845, the non-
vectorized solution should not be run because it never ends.

Thank you for your insights,

Dmitry

• Follow-Ups:
  • Re: Limitation on vector length in NDSolve?
    • From: Oliver Ruebenkoenig <ruebenko@wolfram.com>