# The order of accuracy for my Crank-Nicolson scheme is wrong.

10 views (last 30 days)
David Gillcrist on 23 Jul 2021
I've been trying to figure out what the issue is in my numerical scheme I've written to solve the particle in the box problem for the Schrodinger equation. I've used a Crank-Nicolson scheme to implicitly solve for the wave function numerically, the issue I'm running into is that the Crank-Nicolson scheme is second order accurate in space, but I'm only finding that my numerical solution is on first order accurate. I am near 100% positive that the Crank-Nicolson scheme is order (2,2) for the Schrodinger equation, just like the heat equation, and that there should be no special exceptions here. Because of this, I'm assuming there's an error somewhere in my code for solving of the numerical solution.
My function for solving the particle in a box is here:
function [wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N)
m = ceil((xR - xL)/h); % Number of x-axis intervals
h = (xR - xL)/m;
n = ceil(T/k); % Number of t-axis intervals
k = T/n;
mu = k/h^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% natural units ---> MASS = 1, hbar = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DiagTerm = -2*(1 - 1i*2/mu);
A = diag(DiagTerm*ones(1,m-1)) + diag(ones(1,m-2),1) + diag(ones(1,m-2),-1);
B = -conj(A);
wave_function = zeros(n+1,m-1);
x_val = xL + h*(1:m-1);
wave_function(1,:) = sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (x_val - xL));
for tstep = 2:(n+1)
wave_function(tstep,:) = ((A'*A)\A'*B * wave_function(tstep-1,:)')';
end
wave_function = [zeros(n+1,1), wave_function, zeros(n+1,1)]; % <-- Concatenates Columns to the left and right of wave_function for boundary conditions
wave_function_final = wave_function(end,:);
En = (N*pi)^2/8;
error = sqrt(sum( abs( exp(-1i*En*T)*wave_function(1,2:end-1) - wave_function_final(2:end-1) ).^2));
end
And my live script I use for running is here:
Ms = [5 10 20 40 80];
T = 16/pi;
xL = -1;
xR = 1;
N = 1;
orders = zeros(1,length(Ms)-3);
Errors = zeros(1,length(Ms));
xx = -1:0.05:1;
En = (N*pi)^2/8;
p = plot(xx,real(exp(-1i*En*T)*sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (xx - xL)) ), 'Color', [0 0 0], 'LineStyle',"--", 'LineWidth',2);
p.Color(4) = 0.2;
for j = length(Ms):-1:1
disp(j)
h = 2/Ms(j);
k = 1/Ms(j)^2;
[wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N);
hold on
xx = -1:h:1;
plot(xx,real(wave_function_final),'LineWidth',1.5)
Errors(j) = error;
end
for j = 4:length(Ms)
orders(j-3) = log(abs( (Errors(j) - Errors(j-1))/(Errors(j-1) - Errors(j-2)) )) / log( abs( (Errors(j-1) - Errors(j-2))/(Errors(j-2) - Errors(j-3)) ) );
end
grid on
hold off
disp(Errors)
disp(orders)
figure
plot(Errors)
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

### Categories

Find more on Quantum Mechanics in Help Center and File Exchange

R2020b

### Community Treasure Hunt

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

Start Hunting!