# Wavefunction with fft and ifft in matlab

• MATLAB
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?

## Answers and Replies

What's your initial function f?

And what plot commands are you using?

J77 said:
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

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.

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