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:

  • $\begingroup$ How do you want to specify the sequence of operators? $\endgroup$
    – BlacKow
    Oct 4, 2016 at 3:19
  • $\begingroup$ Look up Composition, maybe? Either that, or ComposeList. $\endgroup$
    – march
    Oct 4, 2016 at 3:43
  • $\begingroup$ Or maybe G[f0, #] & /@ NestList[H, f0, 4]? $\endgroup$
    – march
    Oct 4, 2016 at 3:47
  • 1
    $\begingroup$ Or G[f0, #] & /@ ComposeList[{A, A, B, A, B, B}, x]? Which of these is closest to what you want? $\endgroup$
    – march
    Oct 4, 2016 at 3:48
  • 2
    $\begingroup$ 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. $\endgroup$
    – march
    Oct 4, 2016 at 4:54

1 Answer 1


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},
   Sow[inr[f0, ft = u[ft]]],

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},
   Sow[g[f0, ft = op[order[[j]]][ft]]],
   {j, num}

in which case, for instance,

L[f0, 5] // Column

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.