Parallel computation of NDSolve and a problem

Question

I used this code to solve a system of differential equations:

    ro[t_] := 
    Table[Subscript[ρ, i, j][t], {i, 1, sysdim}, {j, 1, sysdim}];

    RHS = A matrix;
    RHS2 = RHS

    ParallelTable[Flatten[NDSolve[{ro'[t] == RHS, ro[0] == initial}, Flatten[ro[t]], 
                          {t, 0, 10}]], {γ, 1, 2}]

Whenever I use RHSinside the NDSolve, I get this error:

   NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

But, when I use RHS2 which is basically RHS, I don't get any error message.

Does anybody know the problem?

---

Edit.

Here I put the whole code:

    δ[i_, j_] := KroneckerDelta[i, j]
    sysdim = 7;

    Hsystem[dimension_] := Module[{sysdim = dimension},
    V[ii_, jj_] := Module[{i = ii, j = jj}, 
        If[i > j, Return[Subscript[V, j, i]], Return[Subscript[V, i, j]]]];

    hsys[i_, j_] := δ[i, j]*Subscript[ϵ, i] + V[i, j] (1 - δ[i, j]);

    MHsys = Table[hsys[i, j], {i, 1, sysdim}, {j, 1, sysdim}]
    ]

    Hsys = Hsystem[sysdim];

    A[m_] := Table[δ[i, m]*δ[j, m], {i, 1, sysdim}, {j, 1, 
    sysdim}]
    ro[t_] := 
    Table[Subscript[ρ, i, j][t], {i, 1, sysdim}, {j, 1, sysdim}];
    (*LL=Table[(-1+δ[i,j])*ro[t][[i,j]],{i,1,sysdim},{j,1,sysdim}]*)


    L = γ*
    Sum[A[m].ro[t].A[m]\[ConjugateTranspose] - 
     1/2 (A[m].A[m]\[ConjugateTranspose].ro[t] + 
        ro[t].A[m].A[m]\[ConjugateTranspose]), {m, 1, sysdim}];

    Hrecom = Table[-I*ℏ*Γ*δ[i, j], {i, 1, 
    sysdim}, {j, 1, sysdim}];

    Htrap = Table[-I*ℏ*κ*δ[i, 3]*δ[j, 3], {i, 
    1, sysdim}, {j, 1, sysdim}];

    Hdiss = Hrecom + Htrap;

    RHS = -I/ℏ (Hsys.ro[t] - ro[t].Hsys) + -I/ℏ (Hdiss.ro[t] +
       ro[t].Hdiss) + L ;

    hcp = 6.62606957*10^-34*3*10^10*10^-12; Subscript[ϵ, 1] = 
    280*hcp; Subscript[ϵ, 2] = 
    420*hcp; Subscript[ϵ, 3] = 0; Subscript[ϵ, 4] = 
    175*hcp; Subscript[ϵ, 5] = 
    320*hcp; Subscript[ϵ, 6] = 
    360*hcp; Subscript[ϵ, 7] = 260*hcp;

    Subscript[V, 1, 2] = -106*hcp; Subscript[V, 1, 3] = 
     8*hcp; Subscript[V, 1, 4] = -5*hcp; Subscript[V, 1, 5] = 
     6*hcp; Subscript[V, 1, 6] = -8*hcp; Subscript[V, 1, 7] = -4*
     hcp;  Subscript[V, 2, 3] = 28*hcp; Subscript[V, 2, 4] = 
     6*hcp; Subscript[V, 2, 5] = 2*hcp; Subscript[V, 2, 6] = 
     13*hcp; Subscript[V, 2, 7] = 
     1*hcp; Subscript[V, 3, 4] = -62*hcp; Subscript[V, 3, 5] = -1*
     hcp; Subscript[V, 3, 6] = -9*hcp; Subscript[V, 3, 7] = 
     17*hcp; Subscript[V, 4, 5] = -70*hcp; Subscript[V, 4, 6] = -19*
     hcp; Subscript[V, 4, 7] = -57*hcp; Subscript[V, 5, 6] = 
    40*hcp; Subscript[V, 5, 7] = -2*hcp; Subscript[V, 6, 7] = 32*hcp;


    κ = 1;
    Γ = 10^-3;
    ℏ = (6.62606957*10^-34)/(2*Pi);
    initial = Table[δ[i, 6]*δ[j, 6] + δ[i, 1]*δ[j, 
      1], {i, 1, sysdim}, {j, 1, sysdim}];

    RHS2 = RHS;

    num = 20;
    upper\[TripleDot]limit = 10;
    sol = ParallelTable[NDSolve[{ro'[t] == RHS, ro[0] == initial}, Flatten[ro[t]], 
    {t, 0,upper\[TripleDot]limit}, MaxSteps -> 10^5], {γ, 1, num, 1}];//AbsoluteTiming

    dens = Flatten[Table[ro[t] /. sol[[i, All, All]], {i, 1, num}],1];

Answer 1

I decided to write this answer because the equivalence of this and the OP's question may be not immediately obvious. In short, though I'm not sure if it should be called a bug or side-effect, the culprit is the Subscript inside RHS.

I suggest not to use Subscript. List or something like v[1, 1] are more convenient. Personally I think Subscript is a compromise to the traditional form of math symbol and a somewhat pathological object in Mathematica which will easily trigger problems and look ugly when pasting here. There are many posts on this site about the side-effects of Subscript, which you can find by searching. Now let me focus on your problem.

First of all, I'd like to point out that RHS and RHS2 aren't the same. Just remove the semicolon after the definition of RHS and RHS2 and run your code, you'll see the differences:

RHS:

RHS2:

Why do the differences cause an error? To make this answer cleaner I'll use a much simpler sample reproducing your problem (btw, next time you ask a question please put in some effort to reduce your problem into a minimal, workable example to attract more attention and upvotes):

Clear[Subscript]
eqn = y'[x] == Sin[x] + Subscript[a, 1];
Subscript[a, 1] = 1;
ParallelTable[NDSolve[{eqn, y[0] == i}, y, {x, 0, 1}], {i, 1}]

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

Removing Subscript, the problem doesn't exist anymore:

Clear[subscript]
eqn = y'[x] == Sin[x] + subscript[a, 1];
subscript[a, 1] = 1;
ParallelTable[NDSolve[{eqn, y[0] == i}, y, {x, 0, 1}], {i, 1}]
(* This works well. *)

Btw, any symbol that doesn't belong to the Global` context suffers the same problem:

Clear["a`*"]
eqn = y'[x] == Sin[x] + a`a;
a`a = 1;
ParallelTable[NDSolve[{eqn, y[0] == i}, y, {x, 0, 1}], {i, 1}]

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

But in this case the problem can be fixed by the DistributedContexts option:

Clear["a`*"]
eqn = y'[x] == Sin[x] + a`a;
a`a = 1;
ParallelTable[NDSolve[{eqn, y[0] == i}, y, {x, 0, 1}], {i, 1}, DistributedContexts -> All]
(* This also works well. *)

While this won't work on Subscript:

Clear[Subscript]
eqn = y'[x] == Sin[x] + Subscript[a, 1];
Subscript[a, 1] = 1;
ParallelTable[NDSolve[{eqn, y[0] == i}, y, {x, 0, 1}], {i, 1}, DistributedContexts -> All]

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

It seems that DistributedContexts option can't distribute System`. Similar options or functions like DistributeDefinitions etc. don't work either. As I said above, I'm not sure if it's a bug or side-effect i.e. these functions are designed to work like this. Anyway, if you abandon those Subscripts in your code, the problem will disappear.