Numerical solution of Schrodinger equation script

[aaaa202]

I am right now working on a script that solves the Schrodinger equation numerically for arbitrary potentials using the finite difference method. The idea is that I diagonalize the Hamiltonian with elements:
H(i,i+1)=1/dx^2 * constants
H(i,i-1)=1/dx^2 * constants
H(i,i) = -2/dx^2 * constants
and zero elsewhere.
Now what I am not so sure about is what boundary conditions this procedure assumes? I think it assumes automatically that the wave functions tend to zero at the boundaries of your grid.
In any case the script works for 1d, but I am now trying to make it work for 2d so I can solve the circular well problem. However, it has not worked so far:
I have used that the radial equation (polar coordinates) is given as that shown here
http://www.physics.csbsju.edu/QM/square.08.html
which has led me to a hamiltonian of the form (dx is mesh size)
H(i,i) = -2/dx^2 - 1/r(i)^2
H(i,i+1)=1/dx^2 + 1/(dx*r(i))
H(i,i-1)=1/dx^2 - 1/(dx*r(i))
But the solutions I am getting oscillate more and more and don't look anything like the analytic solution.
Can anyone see what the problem is? I suspect it may be the fact that my method forces the wave function to go to zero at the boundaries, which does not make sense for r=0. If this is indeed the problem, what could I do to fix it?

 

[Geofleur]

You need to make sure that your numerical scheme preserves the normalization of the wave function. I am not sure if straight finite differences will do that. You might consider using the method of lines instead to solve the equation. The idea is this: use finite differences, but only on the spatial derivatives. You will end up with a system of first-order ODEs in time, which you can then stuff into your favorite ODE system solver.

 

[cgk]

@OP: Regarding the boundary conditions: If you do not put in any additional information constraints, the (implicit) boundary conditions will be that the wave function is zero everywhere off the grid. That comes simply from the fact that there is no basis there to expand it, and thus your solution will be equivalent to psi(x)=0 in these regions. But you can change that; for example, you can make periodic boundary conditions (even with non-trivial phase angles phi) by simply identifying points off the grid with points on the other side.
E.g., if you have a grid with indices 0, 1, 2, ..., N-1, you'd hack up the Hamiltonian elements:
H(N-1,0) = exp(i phi) <put here what you would normally put into H(N-1,N) if the matrix element existed), e.g., the same as in H(N-2,N-1)>
H(0,N-1) = exp(-i phi) <put here what you would normally put into H(0,-1) if the matrix element existed), e.g., the same as in H(1,0)>
(in place of H(N-1,N), where N is identified with 0 now). By varying the phase angles you would get what corresponds to the band structure of the periodic system.

Regarding the cylindrical coordinates: Your representation of the Laplace operator in cylindrical coordinates does not look right. Could you explain how you got these finite difference formulas? From the first equation of the link you posted I would expect that this is supposed to be
$$\left[\frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r}\right]$$
I see how would take these terms:
H(i,i) = -2/dx^2 - ...
H(i,i+1)=1/dx^2 +
H(i,i-1)=1/dx^2 - ...

but these terms look like they are off by factor of two (the distance between those is 2dx, not dx. Central difference)
H(i,i+1)=...+ 1/(dx*r(i))
H(i,i-1)=... - 1/(dx*r(i))

and this one:
H(i,i) = ...- 1/r(i)^2

I don't see it in the equation at all. Is this supposed to be for the angle derivative? I think if you don't have angles, the angle derivative will be zero, not one.

UPDATE: it also looks like the first order derivative terms could have a wrong sign.

 

[cgk]

I missed the comment about the wave function being 0 on the inside of the cylindrical coordinate systems. So.. some extension.

I probably should have expanded on how the boundary conditions come in: In the end, you just want to solve a linear algebra system---but only in the space you consider as "free" (i.e., without the boundary condition constraints). This means that you need to take into account the boundary conditions already when translating the formal equations into matrix elements, such that your linear algebra system respects them implicitly.

For example, in your original post, you said nothing about what happens for the Hamiltonian matrix elements H(i,i+1) when i is still in grid, but i+1 is not. Normally this matrix element would connect the wave function at point i to the wave function at point i+1:
$$(H \psi)(i) = ... + H(i,i+1) \psi(i+1)$$
but if you simply leave out these matrix elements if i=N-1 is in the grid, but i+1=N is not, that is exactly the same as what you would get if you had set psi(N) = 0 as a boundary condition, since in this case the this contribution would also vanish (just as it does if you simply drop it). You could force other boundary conditions. For example, if you identify
$$\psi(N+i) := \exp(i \phi) \psi(i),$$
where $\phi$ is some arbitrary, but fixed phase angle, you get what is called "twisted boundary conditions" (a form of periodic boundary conditions. You get pure periodic boundary conditions if you just set $\psi(N+i) = \psi(i)$). These boundary conditions on the wave function would translate into how you make the matrix equations, by leading to additional terms in H which connect the two sides of the grid to each other.

You need to do something similar if you want to impose other boundary conditions. For example, if you want to force that your wave function has a value of $\psi(-1)=174$ just off the left side of the grid, you need to put this value of 174 in all equations which otherwise would have contained a reference to the value of $\psi(-1)$ if it had been inside the grid. Now I am not sure what exactly the right boundary conditions in your cylindrical case would be, but you could do things like insert $\psi(0)$ for the value of $\psi(-1)$ (which would be equivalent to asserting that the first radial derivative of the wave function at r=0 is zero, not the wave function itself) and similar things.

 

[aaaa202]

You are saying the term from the angle derivative shouldn't appear. But if you separate variables doesn't it appear in the radial equation? Like in the link in my original post.

 

[cgk]

If you have the phi angle, you need the angle term. But the equations you wrote is how the angle term would look if $\frac{\partial^2}{\partial\phi^2}$ applied to the wave function would give a value of "1" universally. I could understand if you set this term to zero (i.e., looking only for spherically symmetric solutions, which do not depend on phi and thus have phi derivatives zero). But this does not look right to me.

 

[aaaa202]

oh yes sorry. Its not supposed to be 1 universally. It should be m^2, where m is a quantum number.

 

[aaaa202]

But I guess that the problem is that my finite difference scheme doesn't enforce the right boundary conditions. Does anyone have an idea of how to fix this? The boundary conditions are I guess that phi_r(R)=0, where R is the radius of the circular well and phi_r is the radial solution, and since this should hold for all angles the more correct boundary condition is phi(R,Φ)=0. But can I implement this is in my finite difference scheme?