# Generating a list of values from the application of a sequence of operators

I have a module 'L' that reads:

L[f0_,n_]:=Module[{ft=f0}, a[t]={}; Do[AppendTo[a[t], Inr[f0,ft]];
ft=U@ft, {n}]];


It takes an input function f0, applies unitary operator U repeatedly n-times, i.e. it is a n-fold composition of U on f0. The resulting function after each single composition step is evaluated with Inr (Inr is the inner product here). Finally, the value is added to a list a[t], and so on.

What I am looking to achieve is instead of applying U n-times, (U[U[U...[f0]..], I want to construct an action of a repeated sequence of different operators A,B. For example, A[A[A[B[B[B[A[A[A[B...[f0]..], that likewise amounts to a total of n-fold composition. The list a[t] then should be constructed via Inr employing ft's obtained from each single application of A or B along the composition sequence.

Could anyone provide suggestions how to implement this with Do command, or in some other way! Building a list of the operator sequences and using NestList does not seem to be permissible due to memory constraints since all intermediate values ft are retained. Only a[t]'s are of interest here.

The functions are defined as:

U@ft_:=UV InverseFourier[UK Fourier[ft]];


i.e. Going between position/momentum-spaces to evaluate application of a potential energy term UV and momentum UK term, where:

UV:=Exp[-I V dt];
UK:=Exp[-I K dt];


This is the split-operator method solution of a 't'-dependent 2D Schrodinger eq. on a lattice grid. f0 -initial is any discretized 2D function, say a Gaussian. Difference between U and A and B is due to different 'V' terms of 'UV's. Inr on the spatial grid:

Inr[f0_,ft_]:=Dot[Flatten[Conjugate[f0]],Flatten[ft]]dxdy;

• How do you want to specify the sequence of operators? – BlacKow Oct 4 '16 at 3:19
• Look up Composition, maybe? Either that, or ComposeList. – march Oct 4 '16 at 3:43
• Or maybe G[f0, #] & /@ NestList[H, f0, 4]? – march Oct 4 '16 at 3:47
• Or G[f0, #] & /@ ComposeList[{A, A, B, A, B, B}, x]? Which of these is closest to what you want? – march Oct 4 '16 at 3:48
• I see what you mean. I believe this is a job for a Sow/Reap construction, because AppendTo is notoriously slow, and using Sow/Reap allows you to gather desired outputs at any time during a Do loop. However, I think it will be hard to determine if it's useful without a minimal example showing the actual format of your functions/matrices/etc. Obviously, don't post a 1024X1024 grid, but perhaps you could give an example grid (4x4?) with example A and B, and your expression for G. – march Oct 4 '16 at 4:54

The definitions for the functions are as follows:

uv = Exp[-I v dt];
uk = Exp[-I k dt];
u[ft_] := uv InverseFourier[uk Fourier[ft]]
inr[f0_, ft_] := Dot[Flatten[Conjugate[f0]], Flatten[ft]] dx dy;


where it is understood that v, k, and ft all have the same dimension.

Then, do

L[f0_, n_] := Last@Last@Reap@Module[{ft = f0},
Do[
Sow[inr[f0, ft = u[ft]]],
{n}
]
]


If instead we have a sequence of operations, and we want to Sow a function of the result at each step, I would define the operators as an indexed set of operators, then define a list that includes the order of operations, and apply it.

As a generic example, let's suppose we have three operators op[1], op[2], and op[3], all defined as functions or pure functions of some sort. Let's suppose we want to apply these operators in order, each twice, and then repeat n times. Define

order = With[{n = 3}, Join @@ Table[{1, 1, 2, 2, 3, 3}, {n}]]
(* {1, 1, 2, 2, 3, 3, 1, 1, 2, 2, 3, 3, 1, 1, 2, 2, 3, 3} *)


Then, we have some function g that we want to apply to the output of these operations at each step. Then:

L[f0_, num_] := Last@Last@Reap@Module[{ft = f0},
Do[
Sow[g[f0, ft = op[order[[j]]][ft]]],
{j, num}
]
]


in which case, for instance,

L[f0, 5] // Column


For your problem, I would define

op[1] = Fourier;
op[2] = Exp[-I k dt] #&;
op[4] = InverseFourier;
op[3] = Exp[-I v dt] #&;


and (if this makes sense, which I don't really think it does because you will be taking inner products of wave functions that are expanded in different bases, namely position and momentum bases, but you can work that out four yourself)

g[f0_, ft_] := Dot[Flatten[Conjugate[f0]], Flatten[ft]] dx dy


Then, everything should work out.