Saving variable in a nested for loop

2 views (last 30 days)
Nicholas Davis
Nicholas Davis on 17 Aug 2021
Answered: Walter Roberson on 18 Aug 2021
Hi,
The code I have written demonstrates the time-dependent Schrodinger equation evolving through time. My goal with the way it is written now is to demonstrate the effect of the offset, b, on the barrier transmission probability on the final iteration. I want to show this by plotting the transmission probability versus the offset, b. The end result is expected to be sinusoidal.
The issue I'm facing is that I'm unsure of how to save the transmissions as a vector that I can use to plot against my offset since my for loop iterates through different values of the offset, b, rather than through indeces. As it stands now, I only get the very last value for transmission probability for the very last offset value.
I've scoured the forums for help but I haven't seen any for loops that dont use indeces. Any help would be appreciated. Thanks!
clear all
% Initialize Parameters
N = 500; % Number of grid points
L = 100; % Width of axis
x = linspace(-L,L,N); % Grid points
h = (2*L) / N; % Grid point increments
h_bar = 1; mass = 1; w = 1; % Natural Units
tau = 1;
% Potential Barrier
V0 = 0.1;
V = zeros(1,N);
a = 5;
for b = linspace(-8,8,N/4) % N divided by 4 for speed purposes
for i = 1:N
if x(i) <= a+b && x(i) >= -a+b
V(i) = V0;
end
end
% Set up the Hamiltonian Operator Matrix
ham = zeros(N); % Set all elements to zero
coeff = -h_bar^2/(2*mass*h^2);
for i = 2:(N-1)
ham(i,i-1) = coeff;
ham(i,i) = -2*coeff + V(i);
ham(i,i+1) = coeff;
end
% First and last rows for periodic boundary conditions
ham(1,N) = coeff; ham(1,1) = -2*coeff; ham(1,2) = coeff;
ham(N,N-1) = coeff; ham(N,N) = -2*coeff; ham(N,1) = coeff;
% Compute the Crank-Nicolson matrix
dCN = ( inv(eye(N) + .5*j*tau/h_bar*ham) * (eye(N) - .5*j*tau/h_bar*ham));
% Initialize the wavefunction
x0 = 0; % Location of the center of the wavepacket
velocity = 0.5;
k0 = mass*velocity/h_bar; % Wavenumber
sigma0 = 5; % Width of wavefunction
Norm_psi = 1/(sqrt(sigma0*sqrt(pi))); % Normalization
psi = Norm_psi * cos(k0*x') .* exp(-(x'-x0).^2/(2*sigma0^2));
% Time Evolution of Psi
max_iter = L/(2*velocity*tau);
for iter = 1:max_iter
psi = dCN*psi;
end
% Transmission
xR = x;
xL = x;
for i = 1:length(x)
if x(i)>= a+b
xL(i) = 0;
end
if x(i)<= -a+b
xR(i) = 0;
end
end
pFinal = psi.*conj(psi);
tR = trapz(xR,pFinal);
tL = trapz(xL,pFinal);
transmission = abs(tR - tL)/2;
end
% Plot probability versus position at various times
figure(1)
plot(x,pFinal);
xlabel('x'); ylabel('P(x,t)');
title('Probability density at various times');
r = rectangle('Position',[(-a+b) 0 2*(a) max(max(pFinal))]);
r.LineStyle = '-.';

Answers (1)

Walter Roberson
Walter Roberson on 18 Aug 2021
% Initialize Parameters
N = 500; % Number of grid points
L = 100; % Width of axis
x = linspace(-L,L,N); % Grid points
h = (2*L) / N; % Grid point increments
h_bar = 1; mass = 1; w = 1; % Natural Units
tau = 1;
% Potential Barrier
V0 = 0.1;
V = zeros(1,N);
a = 5;
bvals = linspace(-8,8,N/4);
num_b = length(bvals);
for bidx = 1 : num_b % N divided by 4 for speed purposes
b = bvals(bidx);
for i = 1:N
if x(i) <= a+b && x(i) >= -a+b
V(i) = V0;
end
end
% Set up the Hamiltonian Operator Matrix
ham = zeros(N); % Set all elements to zero
coeff = -h_bar^2/(2*mass*h^2);
for i = 2:(N-1)
ham(i,i-1) = coeff;
ham(i,i) = -2*coeff + V(i);
ham(i,i+1) = coeff;
end
% First and last rows for periodic boundary conditions
ham(1,N) = coeff; ham(1,1) = -2*coeff; ham(1,2) = coeff;
ham(N,N-1) = coeff; ham(N,N) = -2*coeff; ham(N,1) = coeff;
% Compute the Crank-Nicolson matrix
dCN = ( inv(eye(N) + .5*j*tau/h_bar*ham) * (eye(N) - .5*j*tau/h_bar*ham));
% Initialize the wavefunction
x0 = 0; % Location of the center of the wavepacket
velocity = 0.5;
k0 = mass*velocity/h_bar; % Wavenumber
sigma0 = 5; % Width of wavefunction
Norm_psi = 1/(sqrt(sigma0*sqrt(pi))); % Normalization
psi = Norm_psi * cos(k0*x') .* exp(-(x'-x0).^2/(2*sigma0^2));
% Time Evolution of Psi
max_iter = L/(2*velocity*tau);
for iter = 1:max_iter
psi = dCN*psi;
end
% Transmission
xR = x;
xL = x;
for i = 1:length(x)
if x(i)>= a+b
xL(i) = 0;
end
if x(i)<= -a+b
xR(i) = 0;
end
end
pFinal = psi.*conj(psi);
tR = trapz(xR,pFinal);
tL = trapz(xL,pFinal);
transmission(bidx) = abs(tR - tL)/2;
end
% Plot probability versus position at various times
figure(1)
plot(x,pFinal);
xlabel('x'); ylabel('P(x,t)');
title('Probability density at various times');
r = rectangle('Position',[(-a+b) 0 2*(a) max(max(pFinal))]);
r.LineStyle = '-.';
plot(bvals, transmission); xlabel('b'); ylabel('transmission')

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!