Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

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

Search the Archive

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



  • Prev by Date: Re: Export
  • Next by Date: Re: Plot function with two arguments
  • Previous by thread: Re: Export
  • Next by thread: Re: Limitation on vector length in NDSolve?