# Animated Wave Propagation using Fourier & InverseFourier

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).

-
any reason why you write InitialState instead of initialState and F2 instead of f2 etc..? It makes look like they are Mathematica own symbols and can be confusing. – Nasser Oct 19 '13 at 2:25
sorry about that, Ill go ahead and fix it now – Steven Sagona Oct 19 '13 at 2:39
I suggest you define a function nextState that accepts initialState (and optionally the kernel) as an argument and then use NestList to create the list of successive states. This will get rid of the For loop and the cache. – Peltio Oct 19 '13 at 11:24

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.

-
What an elegant solution! Thanks a lot! Do you have any suggestions for reducing the noise? As mentioned in the previous post (mathematica.stackexchange.com/questions/34269/…), someone suggested I use ListDeconvolve. But I feel like I'm swinging in the dark deciding what kernel or method to use to deconvolve the list. – Steven Sagona Oct 20 '13 at 4:09
You are welcome. Please make sure to use Abs only where it is... ABSolutely needed. As for what method to use, since it is not clear to me what is the purpose of the operation (I am not familiar at all with split-operatore methods) I cannot be of further assistance. – Peltio Oct 21 '13 at 13:32
@Peltio,I try to use this command for the modulus square for equal intervals of of width Δt but this is not the plot I expect. SquareNorm[x_] := Norm[x]^2; ListLinePlot[Map[SquareNorm, listof[Psi]]] – David H May 3 '15 at 11:13
I just implemented a way to compute what was asked --- oh WTF, I cannot flip back and forward from one question to the other each time. Somebody reopen the original question - marking that as a duplicate is like marking all question using NestList a duplicate!! – Peltio May 3 '15 at 17:33
@Pelito, I've asked them to reopen the question but no body cares. I believe the question was not a duplicate. – David H May 3 '15 at 21:02