# Parallel computation of NDSolve and a problem

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];

• Where's the definition of sysdim and A and matrix? – xzczd May 3 '14 at 13:19
• sysdim is just a number. Buy A Matrix I meant that RHS is a matrix. I thought it would make the question simpler. I will upload all the code soon. – MOON May 3 '14 at 22:56
• @ Kuba. How did you enter the greek letters? Using LaTex? – MOON May 6 '14 at 13:12
• No, using this addon. P.s. do not put space after @. – Kuba May 6 '14 at 13:15
• It works just fine when you use Table, rather than ParallelTable. If you really need to parallelize this code, you'll probably need to take a close look at techniques for resource sharing. – Mark McClure May 6 '14 at 13:50

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] + aa;
aa = 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] + aa;
aa = 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.

• Thank you for your response. You are right the problem is using subscript, which is really convenient. I put the assignment of the variables of RHS after it, and when you execute those assignments, after that RHSand RHS2 are the same. I used the subscripts again but I moved the assignment of the variables of RHSto before the definition of RHS. Now, if I use RHSto solve the differential equations, there won't be any problem. – MOON May 9 '14 at 13:44
• It seems ParallelTablehave some problems with subscript and there is not any problem with subscript and Table, however, it's not clear to me why. If the assignments are after the definition of RHS, it doesn't have any variable anymore(except for gama which a table is made by it). So, the question is that why can't ParallelTable see that there is an assignment? and how does assigning RHSto another variable resolve this issue? – MOON May 9 '14 at 13:56
• @yashar It's because ParallelTable has the attribute HoldAll, so those Subscripts are passed into it before changing into numeric values. When you use RHS2 to store RHS, evaluation happens, and Subscripts are removed, it's equivalent to use Evaluate@RHS inside ParallelTable. – xzczd May 9 '14 at 14:05