# Physics Simulation Forum

 Page 1 of 1 [ 1 post ]
 Print view Previous topic | Next topic
Author Message
 Posted: Tue May 01, 2012 6:49 pm

Joined: Tue May 01, 2012 6:43 pm
Posts: 1
Dear friends,
I have written a matlab code to solve the two-dimensional schrodinger equation in a 2D Harmonic potential.
But the program gives diverging result.
Can anybody help me to find out the bug in the program. I am attaching the Matlab program.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
m=200;
n=m;
T=50;
dt=0.2;
hbar=1;
me=1;
xstart=-2;
xend=2;
ystart=-2;
yend=2;
dx=(xend-xstart)/(m-1);
dy=(xend-xstart)/(m-1);
x=(xstart:dx:xend);
y=(ystart:dy:yend);
%%%%%%%Two-Dimensional_Potential%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xx,yy]=meshgrid(x,y);
Vxy=2*(xx.^2+yy.^2);
%figure,surf(Vxy)
%%%%%%%%%%%%%%%%%%initial wave packet%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vecini=0.01*exp(-xx.*xx-yy.*yy);
vec=vecini;
%figure,surf(vecini)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dx=sparse(n-1,n-1);
Dy=sparse(n-1,n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rx=hbar*dt/(4*me*dx*dx);
ry=hbar*dt/(4*me*dy*dy);
cv=dt/(2*hbar);
%%%%%%%%%%%%%%%%%%%Tridiagonal Matrices%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dx=diag(rx*ones(1,n))+diag((i-2*rx)*ones(1,n-1),-1)+diag(rx*ones(1,n-1),1);
Dy=diag(ry*ones(1,n))+diag((i-2*ry)*ones(1,n-1),-1)+diag(ry*ones(1,n-1),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vec(1,:)=0; %Dritchlet boundary conditions
vec(n,:)=0;
vec(:,1)=0;
vec(:,n)=0;
for tt=1:T
%starting peaceman ADI method%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Y-direction derivative%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:length(x)-1
g(1)=0;%Dritchlet boundary conditions
g(n)=0;
for k = 2:length(y)-1
g(k)=(i+2*ry)*vec(j,k)-ry*vec(j,k+1)-ry*vec(j,k-1)+dt/2*Vxy(j,k)*vec(j,k);
end
vecn(:,j)=Dx\g';
end
vecn(:,1)=zeros(n,1);
vecn(:,n)=zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%X-direction derivative%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 2:length(y)-1
g(1)=0;
g(n)=0;
for j = 2:length(x)-1
g(j)=(i+2*rx)*vecn(k,j)-rx*vecn(k,j+1)-rx*vecn(k,j-1)+dt/2*Vxy(k,j)*vecn(k,j);
end
vecn2(:,k-1)=Dy\g';
end
vecn2(:,1)=zeros(n,1);
vecn2(:,n)=zeros(n,1);
vec=vecn2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
surf(x,y,abs(vec.*conj(vec)))
F(tt)=getframe()
end
movie(F)
figure,surf(x,y,Vxy),shading flat,title('2D harmonic potential')
figure,surf(x,y,vecini),shading flat,title('initial gaussian wavepacket')

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 1 post ]

#### Who is online

Users browsing this forum: No registered users and 2 guests

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot post attachments in this forum

Search for:
 Jump to:  Select a forum ------------------ BULLET PHYSICS LIBRARY USERS    General Bullet Physics Support and Feedback    Release Announcements    Applications, Games, Demos or Movies using Bullet PHYSICS AUTHORING TOOLS, SERIALIZATION AND STANDARDS    Physics authoring tools, serialization, standards and related topics RESEARCH AND DEVELOPMENT IN COLLISION DETECTION & PHYSICS. Don't post Bullet support questions here!    Research and development discussion about Collision Detection and Physics Simulation    Links, Papers, Libraries, Demos, Movies, Comparisons       Non-technical forum and license/patent discussion    Career Opportunities