# NDSolve for a large system of simple ODEs

I am solving a system of many (more than 100) ODEs. It is the kind of standard rate equation encountered in semiconductor physics. Here is the system:

eqns := Table[
P[k]'[t] ==
Sum[-P[k][t] P[q][t] a[k,q,t] - P[k][t] (1 - P[q][t]) b[k,q,t] +
(1 - P[k][t]) P[q][t] c[k,q,t] + (1 - P[k][t]) (1 - P[q][t]) d[k,q,t],
{q, Drop[Table[i, {i, 1, M}], {k}]}]
- e[k,t] (P[k][t])^2 + f[k,t] (1 - P[k][t])^2,
{k, 1, M}];


where P[k][t] represent probabilities (bounded between 0 and 1).

The form of the solution for every P[k][t] is very close to an analytical one, and it is very smooth (basically has just one maximum, quite like a gaussian). I am interested in gettign Sum[P[k][t], {k, 1, M}].

I solve the system with

functionlist = Table[P[k], {k, 1, M}];
initialconds = Table[P[k][t0] == initialvalue[[k]], {k, 1, M}];
eqnstosolve = Join[eqns, initialconds];

NDSolve[eqnstosolve, functionlist, {t, t0, tfin},
Method -> {"EquationSimplification" -> "Solve"}];


Question: is there a better way? I tried several method and this seems the faster. However, for my large system of equation evaluation takes hours.

I would like, for example, to use a very fast and sloppily method at the beginning and then to refine the peak. Or I would just like to solve faster than this method by adjusting the steps.

So far everything I've tried just make things worse. I not even sure what EquationSimplification->Solve really does.

Any help is appreciated.

EDIT: since I have suggested to paste a running code, here is my actual code. The problem is what I stated previously in a more polished way:

Ns = 40;
\[Alpha] = 0.01;
T=0.1;
v=0.0001;

kEven[n_, N_] := 2 \[Pi] n/N  - (N - 1) \[Pi]/N;
kPos[Ns_] := Table[kEven[n, Ns], {n, Ns/2, Ns - 1}];
\[Epsilon][g_, k_] := Sqrt[1 + g^2 - 2 g Cos[k]];

Jf = 1.0;
\[Theta][k_, g_] := ArcTan[g - Cos[k], Sin[k]];
J[\[Omega]_, \[Alpha]_, s_, \[Omega]c_] :=
2 \[Alpha] \[Omega]^
s Exp[-\[Omega]/\[Omega]c] HeavisideTheta[\[Omega]];
gOhm[E_, T_] :=
If[E == 0,
2 N[\[Pi]] \[Alpha] T, \[Pi] ((J[(E + 10^(-20)), \[Alpha], 1,
10^6] -
J[-(E + 10^(-20)), 0.01, 1, 10^6]) Exp[(E + 10^(-20))/
T] )/(Exp[(E + 10^(-20))/T] - 1)]

h[t_] := HeavisideTheta[1 - v t] (1 - v t);
tin = -2/v;
relaxfact = 1;
tfin = relaxfact/v;

eqnsPOS := Table[P[k]'[t] ==((1/(2 Ns)) Sum[
-P[k][t] P[q][t] (W[+1, +1, k, q] + W[+1, +1, k, -q])
- P[k][t] (1 - P[q][t]) (W[+1, -1, k, q] + W[+1, -1, k, -q])
+ (1 - P[k][t]) P[q][t] (W[-1, +1, k, q] + W[-1, +1, k, -q])
+ (1 - P[k][t]) (1 - P[q][t]) (W[-1, -1, k, q] +
W[-1, -1, k, -q])
, {q, Drop[Table[i, {i, 1, Ns}], {k}]}]
- (1/(2 Ns)) W[+1, +1, k, k] (P[k][t])^2
+ (1/(2 Ns)) W[-1, -1, k, k] (1 - P[k][t])^2)
, {k, 1, Ns}];

kPos[2 Ns];
Clear[W];
W[s1_, s2_, k_, q_] :=
W[s1, s2, k,
q] = (1 -
s1 s2 Cos[
Sign[k] \[Theta][kPos[2 Ns][[Abs[k]]], h[t]] +
Sign[q] \[Theta][kPos[2 Ns][[Abs[q]]],
h[t]]]) (gOhm[(s1 \[Epsilon][h[t], kPos[2 Ns][[Abs[k]]]] +
s2 \[Epsilon][h[t], kPos[2 Ns][[Abs[q]]]]), T]);

functionlist = Table[P[k], {k, 1, Ns}];

thermalPsFD[j_, g_] := 1/(Exp[\[Epsilon][g, kPos[2 Ns][[j]]]/T] + 1);

thermalPsInit =
Table[1/(Exp[\[Epsilon][h[t], kPos[2 Ns][[k]]]/T] + 1), {k, 1,
Ns}] /. t -> tin;

initialconds = Table[P[k][tin] == thermalPsInit[[k]], {k, 1, Ns}];

eqnstosolve = Join[eqnsPOS, initialconds];

time = AbsoluteTiming[
spos = NDSolve[eqnstosolve, functionlist, {t, tin, tfin},
Method -> {"EquationSimplification" -> "Solve"}];
];
Print["time NDSolve" <> ToString[time[]/60.] <> ToString[" min"]];

• What are values of the intialvalues?
– user21
Apr 27, 2013 at 5:50
• they are actually taken from a fermi-dirac distribution at a given temperature. However for the sake of the example let's fix them at zero. And the functions a,b,c,d,e which are physically scattering/absorption/excitation of phonons rates, can be taken as constant or simple sinusoidal functions of k^2+q^2 for what matters.... Apr 28, 2013 at 1:47
• your will increase the chance that someone will look into this, if you give a complete set of parameters that can easily be copied. You could edit your question to contain all information and code that actually runs.
– user21
Apr 28, 2013 at 6:34
• Thank you, ok I will do it!! I just didn't want to past a too long code. The actual definition of a,b,c,d,e,f is a little long.... I'll notify you later when I do. Thanks! Apr 28, 2013 at 18:33
• it should be a minimal consistent set of code that captures the issue you like to be improved. The easier it is to copy and paste the more people will play with.
– user21
Apr 28, 2013 at 18:36

One option is to reduce the PrecisionGoal and AccuraryGoal. Those are usually set quite high and perhaps this is a viable option here. Setting
, PrecisionGoal -> 4, AccuracyGoal -> 4

for the Ns=10 case runs in 12 seconds, where the original code runs in 20 seconds.
It were good if you removed the calls to N and replaced the real numbers with exact rational numbers and let NDSolve do the coercion.