Animated Wave Propagation using Fourier & InverseFourier

Question

This is a continuation off of previous help on the first part of my project: fourier issue arising from input miscommunication Now I want to go one step further in the current code. Here's the code from earlier with some comments to explain the thought-process:

xdomain = Table[i, {i, -10, 10, .1}];
initialState= E^-2#^2 & /@ xdomain;

f1 = E^-I #^2*Δt & /@ xdomain;

(* For the first step, I want to multiply my initial function by a second function and then fourier transform it. So I'm breaking down both function into a table of points and then fourier-transforming the product of the multiplied-points, because (I think) this is the most efficient way of doing this numerically.*)

f2 = E^-I #^2*Δt & /@ xdomain;

(* The second step is the same principle. We then take this fourier'd result, multiply it by another function, and then Inverse-fourier transform back *)

finalstate=InverseFourier[Fourier[f1*initialState]*f2]; 

(* ^This line is doing the two steps mentioned in the above two comments. *)

My result represents a wave that is propagated 1-unit-delta-t forward in time. To find out what this wave looks like 2-units in time, we need to take our "finalstate" function and use that instead of our intial state. (so to find something 10-units in time forward, I would have to run this 10 times, each time replacing the output with the input.)

Eventually I'd want to have a series of plots that I would assemble to form a movie of the wave's behavior.

This is what I have so far:

xdomain = Table[i, {i, -50, 50, .05}];
initialState = E^-#^2 & /@ xdomain;

cache = {initialState};
(*For 1\[Rule]10 *)
Δt = .001;

f1 = E^-I #^2*Δt & /@ xdomain;
f2 = E^-I #^2*Δt & /@ xdomain;

For[i = 0, i < 15, i++,
 ftot = Abs@InverseFourier[Fourier[f1*initialState]*f2];
 initialState = ftot;
 cache = Join[cache, {ftot}];
 ListLinePlot[
   Partition[Riffle[xdomain, cache[[i + 1]]]/Norm[cache[[i + 1]]], 
    2]] // Print]

So there are two things that I'm looking for answers for:

1. If my method numerically fourier-transforming these functions by turning them into points and multiplying the points individually, makes sense and is the best way of handling it numerically in Mathematica.

2. How can I export these printed for-loop pictures into an animation. (I eventually will decrease the value of delta_t, and increase the amount of pictures taken, so I won't be able to export the animation by hand.)

Thanks again everyone for your help so far! Let me know what you think.

P.S. In case anyone's interested, what I'm constructing is a program that numerically performs the Split-Operator method for solving the wave equation in an arbitrary potential (that commutes with p).

Answer 1

Expanding my comment, you could try something like this. First we define your functions (this is not strictly necessary, but I prefer to see what I am dealing with in a more natural form. Also, you could make Δt a global parameter to avoid carrying around a second argument)

f[x_, Δt_] := Exp[-x^2]Δt
k[x_, Δt_] := Exp[-I x^2]Δt

Let's get discrete:

xdomain = Table[i, {i, -50., 50., .05}];
Δt = .001;
initialState = f[xdomain, Δt];
kernel = k[xdomain, Δt];

The main code I was talking about in my comment is here.

nextState[inState_List] := InverseFourier[Fourier[kernel*inState]*kernel];

(EDIT: I removed the extra Abs[] from this line, and added it later, when we plot our function). This will compute ten subsequent steps:

newStateList = NestList[nextState, initialState, 10];

And this plots the solution:

(ListPlot[Evaluate[Abs[#]], PlotRange -> All]) & /@ newStateList

You can easily modify the code to attach the xdomain list, if you wish to plot the result on a scale with actual x values.

EDIT - I will now expand this answer to show the general case and another example, since the answer I wanted to answer to has been put on hold.

Basically, when you can compute the next state in term of a given state

nextState[previousState_] := functionOf[previousState]

you can use Nest to iterate this computation a number num of times. If you desire a list with all the subsequent states that originated from a given initialState, you can use NestList.

stateEvol[initialState_, num_] := NestList[nextState,initialState, num]

This approach is fairly general and can be used to iterate complex computations. Suppose you start with an initial state ψ[0] given by:

initialψ[hbar_, b_, x_, y_] = b/(2 Sqrt[π]) (x + I y) E^(-
    b (x^2 + y^2)/(4hbar))*{{1}, {1}};

and your next state ψ[0+Δt] can be computed by nextState[previousState_]=M2[Δt].M1[Δt].M2[Δt].previousState, or more precisely with

nextψ[previousψ_] := M1[Δt, v, b, y, hbar].M2[Δt, v, b, y, hbar, kx,
       ky].M1[Δt, v, b, y, hbar].previousψ

where

M1[Δt_, v_, b_, y_, hbar_] := {{1, 0}, {0, 1}}Cos[(Δt v b 
    y)/(2hbar)] - Sin[(
          Δt v b y)/(2hbar)]/((Δt v b y)/(2hbar)){{0, -b y}, {-b y, 0}};
M2[Δt_, v_, b_, y_, hbar_, kx_, ky_] := {{1, 0}, {
          0, 1}} Cos[(Δt v Sqrt[kx^2 + ky^2])/(2hbar)] - Sin[
            Δt v Sqrt[kx^2 + ky^2]]/(Δt v Sqrt[kx^2 + ky^2]){{0, 
      kx - I*ky}, {kx + I*ky, 0}};

then, once you have defined

ψEvolution[initialState_, num_] := NestList[nextψ, initialState, num]

the list of the evolving states {ψ[0],ψ[Δt],ψ[2Δt],...,ψ[nΔt]} is given by

ψEvolution[initialψ, n]

where n is the number of successive iterations. Now, the tricky part is that this compounding of computations can create big symbolical monsters, so one has better make the result of the computation the more compact as possible - a numeric expression if possible (as this should be the goal of this iteration method).

Here's a dirty way to do that when all the parameters specified are constants (I am using the first numerical values that came to mind - and after all in natural units you can see h or hbar equal to one, right? Anyway to get meaningful results one has to specify suitable values for all parameters. The GIGO rules applies here):

listofψ = Block[
    {hbar = 1, b = 1, v = 1, x = 1, y = 1, kx = 1, ky = 1, Δt = .01},
    ψEvolution[initialψ[hbar, b, x, y], 10]
]

This will return a list of numerical state vectors. the first being initialψ, that is ψ[0], and the last being ψ[0 + n*Δt]. It is now trivial to compute the square module of these ψs and plot them on a graph where the points are separated by equal intervals of width Δt.