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 RHS
inside 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];
sysdim
andA
andmatrix
? $\endgroup$ – xzczd May 3 '14 at 13:19A Matrix
I meant thatRHS
is a matrix. I thought it would make the question simpler. I will upload all the code soon. $\endgroup$ – MOON May 3 '14 at 22:56Table
, rather thanParallelTable
. If you really need to parallelize this code, you'll probably need to take a close look at techniques for resource sharing. $\endgroup$ – Mark McClure May 6 '14 at 13:50