4
$\begingroup$

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[[1]]/60.] <> ToString[" min"]];
$\endgroup$
  • 1
    $\begingroup$ What are values of the intialvalues? $\endgroup$ – user21 Apr 27 '13 at 5:50
  • $\begingroup$ 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.... $\endgroup$ – Davide Venturelli Apr 28 '13 at 1:47
  • 1
    $\begingroup$ 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. $\endgroup$ – user21 Apr 28 '13 at 6:34
  • $\begingroup$ 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! $\endgroup$ – Davide Venturelli Apr 28 '13 at 18:33
  • $\begingroup$ 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. $\endgroup$ – user21 Apr 28 '13 at 18:36
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you Ruebenko this helped a lot! $\endgroup$ – Davide Venturelli Jun 5 '13 at 0:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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