MATLAB Wavefunction with fft and ifft in matlab

[JohanL]

I start with calculating psi(x,0). Then i use fft(psi(x,0)) to get the Fourier coefficients. and with fftshift i get the terms with negative index in the beginning. Then i multiply with the timedependent factor and use ifftshift and ifft to get the wavefunction psi(x,t).
Here is the code:

N=4096;
L=30;
k0=5.0;
x=linspace(-L,L,N);
t=linspace(0,5,N);

psi0=f.*exp(im.*k0.*x); % where f is a function that gives the initial condition

c=fft(psi0);
cs=fftshift(c);
n=-N/2:(N/2-1);

fftPSI=cs.*exp(im.*pi^2.*n.^2.*t);

fftPSIs=ifftshift(fftPSI);
PSI=ifft(fftPSIs);

But then when i plot i only get nonsense.
Can somebody please help me with this?

 

[J77]

What's your initial function f?

And what plot commands are you using?

 

[JohanL]

What's your initial function f?

And what plot commands are you using?

Im scaling the varibles with
x'=x/a and t'=2ma^2/h-bar

fftPSI was wrong in my first post...this is the new one:
fftPSI=cs.*exp(-im.*pi^2.*n.^2.*t/900);

This all the code with the plot commands in the end and also the initial function f on row 8.

N=4096;
L=30
dt=0.25;
im=0+1i;
k0=5.0;
x=linspace(-L,L,N);

f(1:N)=0;
I=find(abs(x)>=2 & abs(x)<2.5);
f(I)=1.0;
psi0=f.*exp(im.*k0.*x);

c=fft(psi0,N);
cs=fftshift(c);
n=-N/2:(N/2-1);

t=linspace(0,5,N);

fftPSI=cs.*exp(-im.*pi^2.*n.^2.*t/900);

fftPSIs=ifftshift(fftPSI);
PSI=ifft(fftPSIs,N);



clf
subplot(3,1,1);
plot(x,f,'b')
hold on
plot(x,real(psi0),'g')
TITLE('real(psi0) och f')
axis([-30 30 -1.5 1.5])

subplot(3,1,2);
plot(x,imag(psi0),'r')
hold on
plot(x,f,'b')
TITLE('imag(psi0) och f')
axis([-30 30 -1.5 1.5])

subplot(3,1,3);
A=abs(psi0);
plot(x,A,'k')
hold on
plot(x,f,'b')
TITLE('abs(psi0) och f')
axis([-30 30 -1.5 1.5])

pause

clf
hold on
s=[0.1 0.25 0.5 1.0];
for t=s
subplot(3,1,1);
plot(x,f,'b')
hold on
plot(x,real(PSI),'g')
TITLE('real(PSI) och f')
axis([-30 30 -1.5 1.5])

subplot(3,1,2);
plot(x,imag(PSI),'r')
hold on
plot(x,f,'b')
TITLE('imag(PSI) och f')
axis([-30 30 -1.5 1.5])

subplot(3,1,3);
A2=abs(PSI);
plot(x,A2,'k')
hold on
plot(x,f,'b')
TITLE('abs(PSI) och f')
axis([-30 30 -1.5 1.5])

pause

end

 

[dimensionless]

I'm not totally sure what your doing but I noticed a couple of things. The first, you want an out of psi(x,t). This is two dimensions, so rather than having a 1 by x vector you should have a m by n matrix. I don't know if this is where the porblem is, but replacing the line

fftPSI=cs.*exp(-im.*pi^2.*n.^2.*t/900);

with

fftPSI=cs'*exp(-im.*pi^2.*n.^2.*t/900);

will get you a matrix.

The second, the plots do not look like total nonsense. It all follows from the code

Also if you type

plot(real(cs))

you should definitely be able to see some order.

 

[dimensionless]

I'm just curious. What are you trying to do?