Efficient Dyson series implementation

I'm trying to implement a Dyson series \begin{array}{lcl} U(x,x_0) & = & 1 + \int_{x_0}^{x}{dy_1V(y_1)}+\int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2V(y_1)V(y_2)}}+\cdots \\ & &{} + \int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2 \cdots \int_{x_0}^{y_{n-1}}{dy_nV(y_1)V(y_2) \cdots V(y_n)}}} +\cdots \end{array} where $V(x)$ is a large matrix. (I'm aware of the solution for symmetric functions, it doesn't work in my case.)

Probably the main bottleneck is the integration but, since it involves multiple integrals, I don't see any easy way of dealing with it. Compile seems to be out of reach as well. I'll greatly appreciate any kind of advice.

Here's the code I've written so far:

myInt[exp_, x1_, rat_] :=
Module[{mulInt, len, int, sum, s1, s2, rar, x, x2 = rat x1, r2,
r1},
SetAttributes[x, NHoldFirst];
(*to prevent switching into numerical values*)

mulInt[y1_, y2_, n_Integer /; n > 0] := Module[{tab},
tab = Table[{x[i], y1, x[i - 1]}, {i, 1, n - 1, 1}];
Return@Prepend[tab, {x, y1, y2}]
];
(*mulInt creates lists of variables to be integrated*)

int[n_] :=
Integrate[Dot @@ (#1 /@ #2[[;; , 1]]), Sequence @@ #2] &[exp[#] &,
mulInt[x1, x2, n]];
(*maps first elements of mulInt's lists as variables of exp,
creates a matrix product out of all exps*)

sum[n_] := Module[{set},
set = Sum[int[m], {m, 1, n}];
len = Length[set[]];
Return[IdentityMatrix[len] + set]
];
(*Dyson series*)

s1 = sum;
rar = Total[Abs@Flatten[#]] &;
r1 = rar[s1];
(*a function for convergence check*)
Do[s2 = s1 + int[n];
r2 = rar[s2];
Print[{n, r2 - r1}];
If[Abs[r2 - r1] <= 10^-3. len^2, Break[]]; s1 = s2;
r1 = r2;, {n, 4, 15}];
(*convergence test*)

Return[s2]
]

Here's a simple example with the 3rd Pauli matrix:

MatrixExp[Integrate[PauliMatrix, {x, 1., 2.}]]
myInt[PauliMatrix &, 1., 2.]
(*results:
{{2.718282, 0.}, {0., 0.3678794}}
{{2.716667, 0}, {0, 0.3666667}}
*)

Of course, MatrixExp[Integrate[function, {x, 1., 2.}]] does not work in more general cases.

The operators I work with are sparse arrays (code below) and it take ages to get the result. Typical values of M: 10 to 20 (M=2 can be a good working example). To run the code, type for instance myInt[toMT[ #, 2]& , 10^-7, 2.].

toMT[x_, M_] :=
Module[{nMax = 8 M + 4, sw, ret, v = 0.1, t = 5 10^8., ntol, mrul},
ntol[nn_, Mn_] := Mn - (Ceiling[nn/4] - 1.);
mrul[n_, m_, y_] := {n, m} -> y;

sw[n_] :=
Switch[{Mod[n, 4], n + 12 <= nMax, n + 15 <= nMax, n - 12 >= 1},
{1, True, True, _},
{mrul[n, n, -ntol[n, M] /x],
mrul[n, n + 12, I v  t],
mrul[n, n + 15, 2 v   (ntol[n, M] - 2) /x]},

{1, True, False, _},
{mrul[n, n, -ntol[n, M] /x],
mrul[n, n + 12, I v  t]},

{1, False, _, _},
{mrul[n, n, -ntol[n, M] /x]},

{2, _, _, _},
{mrul[n, n, (ntol[n, M] - 1) /x],
mrul[n, n + 1, -I t ]},

{3, _, _, True},
{mrul[n, n, ntol[n, M] /x],
mrul[n, n - 12, I t v  ],
mrul[n, n - 13, -2 v  (ntol[n, M] + 2) /x]},

{3, _, _, False},
{mrul[n, n, ntol[n, M] /x]},

{0, _, _, _},
{mrul[n, n, -(ntol[n, M] + 1) /x],
mrul[n, n - 3, -I t ]}
];
ret = Map[sw[#] &, Range[1, 8 M + 4]];
Return[SparseArray@Flatten[ret]]]
• Closely related question: Solving a Volterra integral equation numerically. Doing the iteration to some finite order is often useful to derive perturbative (analytic) results, but numerically it's really better to avoid the iterated integrals, as @ybeltukov did. – Jens Sep 5 '14 at 18:29

One can decrease the difficulty of the problem by reducing the Dyson series to a matrix ODE. Let's start from the definition

$$U(x,x_0) = 1 + \int_{x_0}^{x}{dy_1V(y_1)}+\int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2V(y_1)V(y_2)}}+\ldots$$

and take the derivative with respect to $x$

$$\frac{\partial}{\partial x}U(x,x_0) = V(x)\left(1+\int_{x_0}^{x}{dy_2V(y_2)}+\ldots\right).$$

It has the same tail as $U(x,x_0)$, therefore $$\frac{\partial}{\partial x}U(x,x_0) = V(x)U(x,x_0).$$

Now it is very easy to solve it with NDSolve

σ3 = PauliMatrix;
u0 = IdentityMatrix@Length@σ3;
{x1, x2} = {1, 2};
NDSolve[{u'[x] == σ3.u[x], u[x1] == u0}, u, {x, x1, x2}][[1, 1, 2]][x2]

{{2.71828, 0.}, {0., 0.367879}}

m = 20;
u0 = IdentityMatrix@Length@toMT[x, m];
{x1, x2} = 10^-7 {1, 2};
res = NDSolve[{u'[x] == toMT[x, m].u[x], u[x1] == u0},
u, {x, x1, x2}][[1, 1, 2]][x2]; // AbsoluteTiming

{5.583331, Null}

It takes about 5 sec to obtain res for M=20!

• Looks too simple! :-) Checking. – Gregory Rut Sep 5 '14 at 13:54
• Works! Thanks you! – Gregory Rut Sep 5 '14 at 16:09
• @GregoryRut You are welcome! =) – ybeltukov Sep 5 '14 at 16:11