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 somewhere between first and second order accurate. I am nearly 100% positive that the Crank-Nicolson scheme is of 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. I'm determining the order of accuracy by checking if the following ration holds.
$$
\frac{O\left(\left(\dfrac{1}{2} \Delta x\right)^2\right)}{O(\Delta x^2)} \approx \frac{1}{4}
$$
The scheme I derived from the Schrodinger equation
$$
i \partial_t \psi(t,x) = -\frac{1}{2} \partial_x^2 \psi(t,x)
$$
is given by
$$
\textbf{A}\psi^{n+1} = -\textbf{A}^{*}\psi^{n}
$$
And $\textbf{A}$ is the tridiagonal matrix given by
$$
\textbf{A}_{jk} =
\begin{cases}
-2 + i \cdot 2/\mu, & j = k\\
1, & |j-k|=1\\
0, & \textrm{otherwise}
\end{cases}
$$
And $\mu = \Delta t/\Delta x^2$. The lack of potential term, $V(x)$, in the original Schrodinger equation is because I'm writing this for the particle in the box problem, and the potential's conditions are enforced via the boundary conditions, since there is no tunneling outside the "walls of the box".
My scheme is written in MATLAB 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,1:end) - wave_function_final(1:end) ).^2));
end
In my function $\Delta t = k$ and $\Delta x = h$. And I run it using the script
xL = -1;
xR = 1;
N = 1;
rev = 2*pi;
En = (N*pi)^2/8;
T = rev/En;
M = [10, 20, 40, 80];
Errors = zeros(1,length(M));
orders = zeros(1,length(M)-1);
tic
parfor i=1:length(M)
h = 1/M(i);
k = 1/10000;
[wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N);
Errors(i) = error;
end
for i=1:length(orders)
orders(i) = Errors(i+1)/Errors(i);
end
disp(orders)
The values in my vector orders are all on approximately equal to $0.35$. When in reality these should be approximately $0.25$. I keep my value of $k$ $(\Delta t)$ low so that the error changing should only be that of space and not time. Ultimately, I don't know what's going on as I cannot personally find any problems in the code but I know my order of accuracy should still be much better.
Scheme Background (indirectly solving the Schrodinger equation by means of the evolution of initial wave state).
The choice for the time evolution is $e^{-i H t/\hbar}$,using Cayley's approximation of an exponential gives
$$
\hat{U}(t) = e^{-i H t/\hbar} \approx \frac{1 – \frac{1}{2\hbar} i H t}{1 + \frac{1}{2\hbar} i H t}
$$
which is unitary. Solving the Schrodinger equation at this points amounts to numerically solving the ODE:
$$
\psi(t+\Delta t,\vec{x}) = \hat{U}(\Delta t) \psi(t,\vec{x}) = \frac{1 – \frac{1}{2\hbar} i H \Delta t}{1 + \frac{1}{2\hbar} i H \Delta t} \psi(t,\vec{x})
$$
which can be rewritten as
$$
\left( 1 + \frac{1}{2\hbar} i H \Delta t \right) \psi(t+\Delta t,\vec{x}) = \left(1 – \frac{1}{2\hbar} i H \Delta t \right) \psi(t,\vec{x})
$$
EDIT #1: Removed unnecessary derivations, not overly relevant.
EDIT #2: Used clearer indices so not to confuse $i$ the index with $i$ the number
Best Answer
You are computing the Euclidean norm of the vector that contains the pointwise errors. Each component is essentially $O(h^2)$ where $h$ is the spatial stepsize. In your minimal working example the time step has been reduced to the point where it makes a neglible contribution to the error. Therefore your error behaves as $$E(m,h)=\sqrt{m}h^2$$ where $m$ is the length of the vector and $h$ and $mh$ is a constant, the length of your interval. Your order is computed as the ratio $$E(2m,h/2)/E(m,h) = \frac{\sqrt{2m} \frac{1}{4} h^2}{\sqrt{m} h^2} = \sqrt{2} \frac{1}{4} \approx 0.35.$$ You are getting exactly the results that you should.
Your objective is to compute an approximation of a function $$x \rightarrow f(T,x)$$ for some fixed value $t=T$ of the time. You have approximations $a_j$ of $f(T,x_j)$ for equidistant $x_j$, i.e., $x_j = x_0 + jh$. The global error is $$e_j = f(T,x_j) - a_j = O(k^2 + h^2).$$ You have reduced the temporal stepsize $k$ to the point where it is irrelevant. I would like you to compute either $$E_1 = h \sum_{j=0}^{m-1} |e_j|$$ or $$E_2 = \sqrt{h \sum_{j=0}^{m-1} e_j^2}.$$ Both expressions are of class $O(h^2)$ because $mh = x_R - x_L$ is a constant. [Note: You should be using $m$ to compute $h$, rather than the other way around using the ceiling function. It is not a significant issue here, it will show up in a test based on Richardson's techniques.] Your original/new measure is $$\tilde{E} = \sqrt{\sum_{j=0}^{m-1} e_j^2} = O(\sqrt{mh^4}) = O(\sqrt{m}h^2) = O(h^{\frac{3}{2}}).$$