Finite-Difference Solutions of Radial Equations: Handling the Origin for D > 1?

[lpetrich]

I have a problem that I've been unable to find a simple solution for.

When solving Schroedinger's equation for a central force, the nonradial part of the solution can be found with spherical harmonics, but the radial part is much more difficult.

$E\psi = - \frac{1}{2m}\left(\frac{d^2\psi}{dr^2} + \frac{2}{r}\frac{d\psi}{dr} - \frac{l(l+1)}{r^2}\psi\right) + V(r)\psi$

E = energy, V = potential, m = mass, $\psi$ is the field variable.

The problem is the coordinate singularity at the origin.

For zero angular momentum $l$, one can solve for $r\psi$, and one gets an equation without a coordinate singularity. One can easily finite-difference it and then solve it as an eigensystem.

But for nonzero angular momentum, that coordinate singularity cannot be eliminated by this means. Is there some alternative way of handling that coordinate singularity in this case?

 

[bigfooted]

The usual way is to impose a boundary condition at r=0. If it is cylindrically or spherically symmetric, then at r=0 you can impose a symmetry boundary condition, like putting the gradient at r=0 to 0.
In a finite difference scheme, the unknowns at r=0 can now be eliminated and expressed in terms of the unknowns at the neighboring nodes.

 

[lpetrich]

That I understand, but the problem is handling behavior like $r^l$ near the origin.

I'm now trying a different approach: expanding (psi) as $\psi(r) = \sum_k \psi_k r^m J_n(x_k(r/r_{max}))$
where n = dim/2 + l - 1
and m = 1 - dim/2
The x_k's make the J's zero at r = r_max.

The Schroedinger kinetic-energy term is automatically handled, but the potential term requires constructing integrals of the potentials multiplied by pairs of the J's.

I've found Mathematica's NIntegrate useful for high-quality integrals, but useless for fast ones, so I've been doing the integrals with the midpoint rule and caching the Bessel-function evaluations for additional speed.

That works for the harmonic oscillator, but fails miserably for the 1/r potential. Any better ideas?

 

[Redbelly98]

Assuming the ψ is sufficiently well-behaved so that dψ/dr=0 at the origin, then

$$\frac{d\psi}{dr} \approx r \frac{d^2 \psi}{dr^2}$$

where the 2nd derivative on the RHS is evaluated at r=0. Thus you get an r to cancel out the 1/r in the differential equation.