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.

[itex]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[/itex]

E = energy, V = potential, m = mass, [itex]\psi[/itex] is the field variable.

The problem is the coordinate singularity at the origin.

For zero angular momentum [itex]l[/itex], one can solve for [itex]r\psi[/itex], 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 [itex]r^l[/itex] near the origin.

I'm now trying a different approach: expanding (psi) as [itex]\psi(r) = \sum_k \psi_k r^m J_n(x_k(r/r_{max}))[/itex]
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

[tex]\frac{d\psi}{dr} \approx r \frac{d^2 \psi}{dr^2} [/tex]

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.

 

[blue_raver22]

I understand the difficulty you are facing in solving the radial part of the Schroedinger's equation for a central force. The coordinate singularity at the origin adds an extra layer of complexity to the problem. However, there are several approaches that can be used to handle this issue.

One solution is to use a transformation of variables, such as the logarithmic or square root transformation, which can remove the singularity at the origin. This allows for a more straightforward finite-difference solution to be applied.

Another approach is to use a different numerical method, such as the finite-element method, which can handle singularities in a more efficient manner. This method can also be combined with the transformation of variables mentioned above for better accuracy.

Additionally, there are specialized numerical techniques, such as the asymptotic iteration method, that have been developed specifically for solving problems with singularities at the origin.

It is also worth considering if there are any symmetries in your specific problem that can be exploited to simplify the solution. For example, if the potential is spherically symmetric, then the problem can be reduced to a one-dimensional problem, which may be easier to solve.

In conclusion, while the coordinate singularity at the origin does present a challenge in solving the radial part of the Schroedinger's equation for a central force, there are various approaches and techniques that can be used to handle it. It may require some experimentation and trial and error to find the most suitable method for your specific problem.