Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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;

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];
    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[[1]]/60.] <> ToString[" min"]];
share|improve this question
What are values of the intialvalues? – user21 Apr 27 '13 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.... – Davide Venturelli Apr 28 '13 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 '13 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! – Davide Venturelli Apr 28 '13 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 '13 at 18:36
up vote 3 down vote accepted

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.

share|improve this answer
Thank you Ruebenko this helped a lot! – Davide Venturelli Jun 5 '13 at 0:02

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.