I am having some troubles with discretizing a 1D Schrödinger equation via finite differences.
The equation that I am trying to solve numerically is of the following form
$$
-u^{\prime \prime} +\left (\sin(x) – \frac{40}{1 + x^2} \right)u = \lambda u, \quad u \in (0, \infty),\\
u(0) = 0\,.
$$
The spectrum of this problem has band gap structure and I am aware that there will appear spurious eigenvalues when discretizing it using a simple finite difference scheme. To numerically solve the problem, I truncated it to a finite interval $(0,L)$ with $L$ sufficiently large and impose the Neumann condition $u^\prime(L) = 0$. And that is where my problem is. I am not quite sure how to incorporate this condition into the finite difference matrix. My current approach is to define the tridiagonal matrix $A$ as
$$
A = \frac{1}{h^2} \text{tridiag}(-1,2,-1) + h^2 \text{diag}(q(0), q(h), \dotsc, q((N-1)h))\,,
$$
where $q(x) = \sin(x) – 40/(1+x^2)$ is the potential from above, $h$ is the mesh size and $N$ is the number of mesh points. The original problem is then approximately solved by solving the algebraic eigenvalue problem
$$
Au_h = \lambda_h u_h.
$$
I do have high-accuracy approximations of the bands and gaps from two papers that compute the eigenvalues using different techniques but of course, since I am in fact solving a different problem, my results do not match theirs.
How can I add the Neumann boundary conditions to this discretized poblem? Is what I'm done up until now even right? Thanks in advance.
Best Answer
As $u(0)=0$ is fixed, your first equation in the PDE discretization is centered at
x=h, so your second diagonal matrix should contain $[q(h),q(2h)...,q(Nh)]$.You can use $u(L+h)=u(L-h)+O(h^3)$. This means that with the last equation centered at $x=L=Nh$ $$ -\frac{u(L+h)-2u(L)+u(L-h)}{h^2}+q(L)u(L)=\frac{2u(L)-2u(L-h)}{h^2}+q(L)u(L) $$ you have to double the entry $A_{N,N-1}$.