# Speed up NDSolveValue

I'm trying to solve this differential equation to obtain the density matrix elements $$\rho$$.

$$\frac{d\rho}{dt}=-i\Big(H_{0}(t)\rho(t)-\rho(t)H_{0}(t)\Big)$$,

for Hermitian Hamiltonian $$H_{0}(t)=\begin{pmatrix}vt-J\cos(k) & -{\it i}\sin(k)\\ {\it i}\sin(k) & -vt+J\cos(k) \\\end{pmatrix},$$

with initial time $$t_1=h_1/v$$ and final time $$t_2=h_2/v$$, $$v=\{1/100, 1/98,1/96,1/94,....\}$$, $$h_1=-50$$, $$h_2=50$$, $$J=1$$ and $$k=(2m-1)\pi/L$$, $$m=1,2,3,....,L/2$$, $$L=200$$.

my code works properly but it takes a long time to get data even for single $$k$$. I was wondering if you would be able to help me to speed up my code.

Please note that although in this code $$f(i)$$ and $$g(x)$$ do nothing but I need $$f(i)$$ and $$g(x)$$ in the second part of my code.

$PreRead = (# /. s_String /; StringMatchQ[s, NumberString] && Precision@ToExpression@s == MachinePrecision :> s <> "50." &); SetDirectory[NotebookDirectory[]]; Clear["Global*"]; L = 200; J = 1; h1 = -50.; h2 = 50; lstV = {1/100, 1/98, 1/96, 1/94, 1/92, 1/90, 1/88, 1/86, 1/84, 1/82, 1/80, 1/78, 1/76, 1/74, 1/72, 1/70, 1/68, 1/66, 1/64, 1/62, 1/60, 1/58, 1/56, 1/54, 1/52, 1/50, 1/48, 1/46, 1/44, 1/42, 1/40, 1/38, 1/36, 1/34, 1/32, 1/30, 1/28, 1/26, 1/24, 1/22, 1/20, 1/18, 1/16, 1/15, 1/14, 1/13, 1/12, 1/11, 1/10, 1/9, 1/8, 1/7, 1/6, 1/5, 1/4, 1/3, 1/2, 1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 40, 80, 100}; Do[v = lstV[[u]]; Concfile = OpenWrite[ "Concurrence_VS_V,h1=" <> ToString[N[h1, 4]] <> ", h2=" <> ToString[N[h2, 3]] <> ",N=" <> ToString[N[L, 4]] <> ".txt", FormatType -> ScientificForm]; Concfile2 = OpenWrite[ "Concurrence_VS_Inverse_V,h1=" <> ToString[N[h1, 4]] <> ", h2=" <> ToString[N[h2, 3]] <> ",N=" <> ToString[N[L, 4]] <> ".txt", FormatType -> ScientificForm]; Table[f[i], {i, 1, L/2, 1}]; Table[g[x], {x, 1, L/2, 1}]; For[m = 1, m <= L/2, m = m + 1, k = N[((2 m - 1)*\[Pi])/L]; t1 = N[h1/v]; t2 = N[h2/v]; H0[t_] = {{N[v*t - J*Cos[k]], N[-I*Sin[k]]},{N[I*Sin[k]], N[-v*t + J*Cos[k]]}}; Clear[\[Rho]]; Sol1 = NDSolveValue[{\[Rho]'[t] == -I*(H0[t].\[Rho][t] - \[Rho][t].H0[t]), \[Rho][t1] == {{1, 0},{0, 0}}}, \[Rho], {t, t1, t2}]; rhoa = Sol1[t2]; WriteString[Concfile,FortranForm@N[k/\[Pi]], " ", FortranForm@N[Norm[rhoa[[1, 1]]]], "\n"]; WriteString[Concfile2,FortranForm@N[k/\[Pi]], " ", FortranForm@N[Re[rhoa[[1, 2]]]], "\n"]; f[m] = Chop[rhoa[[1, 1]]]; g[m] = Chop[rhoa[[1, 2]]]; ]; Close[Concfile]; Close[Concfile2]; , {u, 1, Length[lstV]}];  • Could you reduce this down to a single instance of NDSolve, with appropriate parameter values, removing all the looping constructs that only clutter the code? Commented Jul 7, 2023 at 14:52 • Since there are three parameters which I need the density matrix elements as a function of them (h2,k,v) for the second part of the code, I think the differential equation can not be reduced to a single instance of NDSolve. Commented Jul 7, 2023 at 19:50 • I meant to say that you should reduce your code to a minimal working example without the file operations and the looping constructs. Choose a single set of values that gives slow execution and limit your question to that. If somebody can speed up that example, you will then be able to apply that solution to your loops yourself. Commented Jul 7, 2023 at 19:52 • Thanks for your comment. Commented Jul 7, 2023 at 20:12 • There is an analytic solution to your problem doi.org/10.1103/PhysRevA.82.032117. You can reduce monochromatically driven system ($f(t)=f_0+f_1\sin(\omega t)$) to linearly driven by setting$f_1=b/\omega$and taking the limit$\omega\rightarrow 0\$. Not sure if it easy, though. Commented Jul 8, 2023 at 11:59

Computational speed can be improved by

1. Using ParallelTable
2. Using ParametricNDSolveValue with "ParametricSensitivity" and "ParametricCaching" turned off
3. Obtaining only solutions actually needed (final value of ρ)
4. Using s in place of v t
5. Setting J = 1 in equations
6. Explicitly expanding the four differential equations

With these changes, the computation is performed by

h = 50;
psol = ParametricNDSolveValue[{v ρ11'[s] == -Sin[k] (ρ12[s] + ρ21[s]),
Sin[k] ρ11[s] + 2 I (-s + Cos[k]) ρ12[s] == Sin[k] ρ22[s] + v ρ12'[s],
Sin[k] ρ11[s] + 2 I (s - Cos[k]) ρ21[s] == Sin[k] ρ22[s] + v ρ21'[s],
v ρ22'[s] == Sin[k] (ρ12[s] + ρ21[s]),
ρ11[-h] == 1, ρ12[-h] == 0, ρ21[-h] == 0, ρ22[-h] == 0},
{ρ11[h], ρ12[h], ρ21[h], ρ22[h]}, {s, -h, h},
{Element[v, Reals], Element[k, Reals]},
Method -> {"ParametricSensitivity" -> None},
Method -> {"ParametricCaching" -> None}]


Then, with k = 1 (for instance)

ParallelTable[Chop@psol[vv, 1], {vv, lstV}]
(* {{0.000095559, 0.00103702 + 0.00579662 I, 0.00103702 - 0.00579662 I, 0.999904},
{0.000116982, 0.00256075 + 0.00705583 I, 0.00256075 - 0.00705583 I, 0.999883}, ... *)


takes about two minutes on my six-processor PC. Note that the computation is inherently slow for small v, for which the solution oscillates rapidly over {s, -h, h}.

The numerical results above suggest that ρ11[s] + ρ22[s] == 1 and that ρ21[s] == Conjugate[ρ12[s]]. And, a brief examination of the four ODEs above demonstrate that this is so. On this basis, the numerical computation can be performed in about half the time.

h = 50;
psol1 = ParametricNDSolveValue[{v ρ11'[s] == -2 Sin[k] Re[ρ12[s]],
Sin[k] (2 ρ11[s] - 1) + 2 I (-s + Cos[k]) ρ12[s] == v ρ12'[s],
ρ11[-h] == 1, ρ12[-h] == 0}, {ρ11[h], ρ12[h]}, {s, -h, h},
{Element[v, Reals], Element[k, Reals]},
Method -> {"ParametricSensitivity" -> None},
Method -> {"ParametricCaching" -> None}]

ParallelTable[Chop@psol1[vv, 1], {vv, lstV}]
(* {{0.0000956634, 0.00103934 + 0.00580275 I},
{0.000116937, 0.00256215 + 0.00705313 I}, ...} *)

• This is nice answer (+1). Commented Jul 8, 2023 at 11:41
• @bbgodfrey Thank yo so much for your valuable comments and code. Commented Jul 8, 2023 at 12:49

We can reduce time to 90s using les dimension by reducing range $$0.1\le v \le 20$$ as follows

Clear["Global*"];
L = 202;
J = 1;
h1 = -50.;
h2 = 50;

lstV = Join[Table[i/10, {i, 1, 10}], Table[i, {i, 2, 20, 2}]]; nmax =
Length[lstV];

Do[v = lstV[[u]];
For[m = 1, m <= L/2, m = m + 5, k = N[((2 m - 1)*\[Pi])/L];
H0[t_] = {{t - J*Cos[k], -I*Sin[k]}, {I*Sin[k], -t + J*Cos[k]}};
Clear[\[Rho]];
Sol[u, m] =
NDSolveValue[{\[Rho]'[
t] == -I/v*(H0[t] . \[Rho][t] - \[Rho][t] . H0[t]), \[Rho][
h1] == {{1, 0}, {0, 0}}}, \[Rho][h2], {t, h1, h2}];

];
, {u, 1, nmax}]; // AbsoluteTiming


Visualization

rho11 = Table[{lstV[[i]] // N, N[((2 m - 1)*\[Pi])/L],
Sol[i, m][[1, 1]] // Abs}, {i, nmax}, {m, 1, L/2, 5}]; rho12 =
Table[{lstV[[i]] // N, N[((2 m - 1)*\[Pi])/L],
Sol[i, m][[1, 2]] // Abs}, {i, nmax}, {m, 1, L/2, 5}]; rho21 =
Table[{lstV[[i]] // N, N[((2 m - 1)*\[Pi])/L],
Sol[i, m][[2, 1]] // Abs}, {i, nmax}, {m, 1, L/2, 5}]; rho22 =
Table[{lstV[[i]] // N, N[((2 m - 1)*\[Pi])/L],
Sol[i, m][[2, 2]] // Abs}, {i, nmax}, {m, 1, L/2, 5}];
{ListPlot3D[Flatten[rho11, 1], Mesh -> None, ColorFunction -> Hue,
PlotRange -> All, InterpolationOrder -> 3,
AxesLabel -> {"v", "k", ""},
PlotLabel -> "|\!$$\*SubscriptBox[\(\[Rho]$$, $$11$$]\)|"],
ListPlot3D[Flatten[rho12, 1], Mesh -> None, ColorFunction -> Hue,
PlotRange -> All, InterpolationOrder -> 3,
AxesLabel -> {"v", "k", ""},
PlotLabel -> "|\!$$\*SubscriptBox[\(\[Rho]$$, $$12$$]\)|"],
ListPlot3D[Flatten[rho21, 1], Mesh -> None, ColorFunction -> Hue,
PlotRange -> All, InterpolationOrder -> 3,
AxesLabel -> {"v", "k", ""},
PlotLabel -> "|\!$$\*SubscriptBox[\(\[Rho]$$, $$21$$]\)|"],
ListPlot3D[Flatten[rho22, 1], Mesh -> None, ColorFunction -> Hue,
PlotRange -> All, InterpolationOrder -> 3,
AxesLabel -> {"v", "k", ""},
PlotLabel -> "|\!$$\*SubscriptBox[\(\[Rho]$$, $$22$$]\)|"]}
`

• Dear Alex thank yo so much for your valuable comments and code. Commented Jul 8, 2023 at 12:48
• @Radmehr You are welcome! Commented Jul 8, 2023 at 13:20