Wavefunction with fft and ifft in matlab

  • MATLAB
  • Thread starter JohanL
  • Start date
  • #1
158
0
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

  • #2
J77
1,083
1
What's your initial function f?

And what plot commands are you using?
 
  • #3
158
0
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
 
  • #4
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 :smile:

Also if you type

plot(real(cs))

you should definitely be able to see some order.
 
  • #5
I'm just curious. What are you trying to do?
 

Related Threads on Wavefunction with fft and ifft in matlab

  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
12
Views
5K
  • Last Post
Replies
0
Views
3K
  • Last Post
Replies
1
Views
5K
Replies
1
Views
4K
Replies
11
Views
2K
  • Last Post
Replies
2
Views
4K
  • Last Post
Replies
2
Views
2K
  • Last Post
Replies
10
Views
5K
Replies
8
Views
7K
Top