# How can this IBVP with regularity boundary conditions be solved?

I have a radial Schrödinger equation for a particle in Coulomb potential:

$$i\partial_t f(r,t)=-\frac1{r^2}\partial_r\left(r^2\partial_r f(r,t)\right)-\frac2rf(r,t)$$

with initial condition

$$f(r,0)=e^{-r^2}$$

and boundary conditions

$$\begin{cases} |f(0,t)|<\infty\\ |f(\infty,t)|<\infty. \end{cases}$$

Trying to solve it, I couldn't come up to any analytical solution, so I tried to solve it numerically.

Simplest what I thought of was finite differences method. But while I can limit domain to make it finite, imposing zero Dirichlet condition at $r=A$ for some $A$, I would need to somehow impose the condition of regularity at $r=0$, where $f$ doesn't have to be neither zero, nor non-zero in general. I don't know how to do this with finite differences.

Another approach I tried was to expand the initial condition in eigenfunctions of the Hamiltonian operator (the RHS of the PDE) and then approximate the solution taking finite number of eigenstates into account. But the expansion appeared to converge extremely slowly, and after taking 200 bound eigenstates (simple because of closed-form solution in terms of Laguerre polynomials) and 70 eigenstates of continuous spectrum (taking those of them which vanish at $r=A$, so as to impose zero boundary condition at $A$; had $A=10$), I still got nothing similar to initial function.

So, the question is: is there any explicit solution for this IBVP (be it closed-form one or in some sort of series, but necessarily explicit and fast convergent)? If no, how can it be solved numerically?

• Just a cue: have you tried looking at this? – derpy Apr 29 '14 at 15:32
• @derpy it's just what I used as the basis of expansion. It's not a solution to initial value problem — just to a time-independent Schrödinger equation boundary value problem. – Ruslan Apr 29 '14 at 15:35
• Got it. I just have this feeling that PDE should be analytically solvable, but right now I'm not in the mood to dig into this much further, sorry. But I think it should have an analytic solution. – derpy Apr 29 '14 at 15:39

Regarding the regularity condition when applying finite differences, some authors make the following approach:

Consider the RHS of the equation at $r \to 0$ and compute the limit with the help of L'Hôpital's rule:

$$\lim_{r\to 0} \left(- \frac{1}{r^2} \partial_r (r^2 f_r) - \frac{2}{r} f \right) = \lim_{r\to 0} \left( - f_{rr} - \frac{2}{r} f_r - \frac{2}{r} f \right) = \lim_{r\to 0} \left( -f_{rr} - 2 f_{rr} - 2 f_r \right),$$

and hence you have:

$$- \frac{1}{r^2} \partial_r (r^2 f_r) - \frac{2}{r} f \approx -3 f_{rr} - 2 f_r.$$

Now, you may discretize your equation at $i = 0$, if $\, r_i = r_0 + i \Delta r$, using this formula, avoiding the geometrical singularity.

Although this is a pretty "rough" way to approximate both the Laplacian and the source term, it's quite well known when applying a finite difference scheme.

I hope this helps you.

Cheers.

• Hmm, this is quite a nice trick, thanks. I'll wait a few days if someone will point to analytical solution though before I accept the answer. – Ruslan Apr 30 '14 at 9:36
• Glad to help, Ruslan! Another issue regarding the numerical approach for this problem is that of imposing Dirichlet boundary conditions on $r = A$. This is quite 'restrictive' and as far as I know, better results are obtained setting $f_r = 0$ in $r = A$, i.e, imposing homogeneous boundary conditions at the truncation distance. Of course, non-dependence on the choice of sufficiently large values of $A$ must be checked. Hope this helps! – Dmoreno Apr 30 '14 at 9:54
• Do Neumann conditions really give better approximation? The wave will similarly reflect off this boundary, so it seems to give no practical difference. – Ruslan Apr 30 '14 at 9:56
• I don't know if they will in this particular case but, at least in my field, which is fluid mechanics, homogeneous Neumann bc (also known as free-stress) give better results than Dirichlet. – Dmoreno Apr 30 '14 at 9:59

Another approach to set a regularity condition is to make a change of variables to switch the type of boundary condition. For this equation, if $f(r)$ is regular at $r=0$, then $\frac{\text{d}}{\text{d}\rho}f(\rho^2)=0$ at $\rho=0$.

We can switch to $g(\rho)=f(\rho^2)$ and rewrite the equation for $g(\rho)$:

$$i\partial_tg(\rho,t)=-\frac1{4\rho^2}\partial_{\rho\rho}g(\rho,t)-\frac3{4\rho^3}\partial_\rho g(\rho,t)-\frac2{\rho^2}g(\rho,t)$$

with boundary conditions:

\begin{cases} g^{(1,0)}(0,t)=0\\ |g(\infty,t)|<\infty, \end{cases}

i.e. we switched from regularity condition at $r=0$ to homogeneous Neumann boundary condition at $\rho=0$. Now this is trivial to setup for finite differences.

Yet simpler approach is to use $q(r,t)=rf(r,t)$. Then the equation transforms to:

$$i\partial_t q(r,t)=-\partial_{rr}q(r)-\frac2rq(r,t),$$

and the regularity condition becomes homogeneous Dirichlet condition (because finite $f$ is multiplied by zero $r$). Now we can't set good boundary condition at infinity though, but for numerical treatment this is not a problem. For finite domain $r\in[0,L]$ we have the following boundary conditions:

\begin{cases} q(0,t)=0\\ q(L,t)=0, \end{cases}

and this problem is also trivial to setup for finite differences.

This may be a wild goose chase, and at the very least it will be a mess, but here are my thoughts. The time-independent Schrodinger equation (TISE) with Coulomb potential can be solved exactly -- the negative-energy solutions are in any quantum mechanics book and the positive-energy solutions are confluent hypergeometric functions of some sort (I believe they are given in Schiff, for instance). So you can write a general solution of the time-dependent equation (TDSE) as a linear combination of all solutions of the TISE each multiplied by exp(-iEt); the coefficients are determined by the initial condition. (You only need the l=0 solutions, of course, since your initial condition is rotationally invariant.) I have not mentioned boundary conditions, and indeed I can only say that I think the above solutions of the TISE are well-behaved at the origin. The positive-energy solutions surely don't go to zero as r->infinity but I think the sum must do so.

• Well the eigenfunctions are certainly well-behaved at origin, but this is not the problem with using this method. The problem is that positive energy spectrum is continuous, and one would have to integrate over corresponding functions to represent the solution. I'd not say it's really simple, a numeric approach like finite differences in truncated space would be much easier to implement. – Ruslan Mar 18 '15 at 17:57