Time dependant Schrodinger Equation Leapfrog Method

[rshalloo]

Hi,
I'm trying to set up a programme to compute the numerical solution to the time dependant schrodinger equation of a ground state wave packet in a harmonic oscillator using the leapfrog method. I've been at it for two weeks trying different methods and I'm starting to get extremely stressed about it because I can't seem to get the loop for the leapfrog correct. The leapfrog method is supposed to calculate the Real values of the solution at time integer time steps dt,2dt,3dt,...,ndt and the imaginary solution at a half time steps dt/2,3dt/2,5dt/2,...
I'm not sure how to incorporate this into the programme although I think that it would be in the initial conditions and not in the loop

The time dependant schrodinger equation is a first order differential equation in time and so should only require one initial value of the wavefunction i.e x when t=0 but I can't seem to get my programme running without three initial conditions.

Please can someone point me in the right direction, for the sake of my sanity;



Here's my code

% 1. First, define the variables and constants to be used in this programme.
%The starting position of the wave packet is the middle of the initial wave
%function
x_start= input('\nWhat is the initial position of the Wave Packet\n>>');
%Size of position step used and size of time step used.
x_step= input('\nWhat position step would you like to use\n>>>');
t_step= input('\nWhat time step would you like to use\n>>>');
%The number of position steps that you will take and the number of time
%steps you will take
num_x_steps= input('\nWhat number of position steps would you like to take\n>>>');
num_t_steps= input('\nWhat number of time steps would you like to take\n>>>');



% 2. Create arrays and Matrices that you will use
%First create the Position array. With 1st element x_first and step x_step
%the length of which is the number of position steps to be taken. To do
%this first calculate the last and first element of the array which will be
%((num_x_steps*x_step)/2)+x_start and -((num_x_steps*x_step)/2)+x_start
%respectiveley
x_last= ((num_x_steps*x_step)/2)+x_start;
x_first= -((num_x_steps*x_step)/2)+x_start;
position= x_first:x_step:x_last;
%Next create the time array for the imaginary solution with first element 0 and step t_step the
%length of which is the number of time steps to be taken.To do
%this first calculate the last element of the array which will be
%(num_t_steps*t_step)+(t_step/2)
t_last= num_t_steps*t_step;
timere= 0:t_step:t_last;
%Next create the time array for the imaginary solution with first element t_step/2 and step (t_step/2) the
%length of which is the number of time steps to be taken.
timeim= t_step/2:t_step:t_last+(t_step/2);
%Next Create the EMPTY array for both realsol and imagsol. They will both have
%num_x_steps rows and 2num_t_steps collumns
realsol= zeros(num_x_steps,num_t_steps);
imagsol= zeros(num_x_steps,num_t_steps);

% 3. Define the Potential Function within the programme itself.
%To start we will use a very simple Harmonic oscillator
potentialfn= 2*(position).^2;

% 4. Set the initial Conditions of the wave function:
%To start this will be a simple gaussian function
initial_wavefn =exp(-(position-x_start).^2);




realsol([1:num_x_steps+1],1)=initial_wavefn;
imagsol([1:num_x_steps+1],1)=0;
realsol([1:num_x_steps+1],2)=initial_wavefn;
imagsol([1:num_x_steps+1],2)=0;
realsol([1:num_x_steps+1],3)=initial_wavefn;
imagsol([1:num_x_steps+1],3)=0;





% 5. Calculate the imagsol and realsol using leapfrog integration.
%Start with imagsol, using the formula from the pseudocode within a for
%loop.

for (tctr=1:num_t_steps)
for (xctr=2:num_x_steps)
realsol(xctr,tctr+1)=realsol(xctr,tctr)-...
(t_step/(2*(x_step)^2))*[imagsol(xctr+1,tctr)-2*imagsol(xctr,tctr)+imagsol(xctr-1,tctr)]...
-t_step*potentialfn(xctr)*imagsol(xctr,tctr);
imagsol(xctr,tctr+1)=imagsol(xctr,tctr)+...
(t_step/(2*(x_step)^2))*[realsol(xctr+1,tctr)-2*realsol(xctr,tctr)+realsol(xctr-1,tctr)]...
-t_step*potentialfn(xctr)*realsol(xctr,tctr);
end
end

 

[NateD]

% 6. Create Plots%First plot the potential functionfigure(1)plot(position,potentialfn)title('Potential Function');xlabel('Position');ylabel('Potential energy');%Next plot the wavefunction as a function of timefigure(2)for (tctr=1:num_t_steps) plot(position,realsol(:,tctr)) hold onendtitle('Ground State Wave Packet');xlabel('Position');ylabel('Real Part of Wave Function');hold off% Finally plot the probability density as a function of timefigure(3)for (tctr=1:num_t_steps) plot(position,(realsol(:,tctr).^2+imagsol(:,tctr).^2)) hold onendtitle('Probability Density');xlabel('Position');ylabel('Probability density');hold offTo incorporate the leapfrog method into your code, you need to modify the for loop you have created for calculating the imaginary and real solutions. The loop should look something like this: for (tctr=1:num_t_steps) %Calculate the Real Solution for (xctr=2:num_x_steps) realsol(xctr,tctr+1)=realsol(xctr,tctr)-... (t_step/(2*(x_step)^2))*[imagsol(xctr+1,tctr)-2*imagsol(xctr,tctr)+imagsol(xctr-1,tctr)]... -t_step*potentialfn(xctr)*imagsol(xctr,tctr); end %Calculate the Imaginary Solution for (xctr=2:num_x_steps) imagsol(xctr,tctr+1)=imagsol(x