Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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

1 Answer 1

Expanding my comment, you could try something like this. First we define your functions (this is not 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 avoide the 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] := Abs[InverseFourier[Fourier[kernel*inState]*kernel]];

(are you sure you want Abs here? Perhaps you want that only when you plot your function) This will compute ten subsequent steps

newStateList = NestList[nextState, initialState, 10];

You can then attach the xdomain list, but let's see if this is what you want before doing that (I am placing Abs here too to show the alternative position):

(ListPlot[Evaluate[Abs[#]], PlotRange -> All]) & /@ newStateList
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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