1

When solving the hydrogen radial Schroedinger equation (with $r > 0$ the radial coordinate) for angular momentum $L=1$ and the modified radial wave function $P(r)=rR(r)$, $P(r)$ satisfies: $d^2P/dr^2=2(−1/r−E+1/r^2) P$, where $E$ is the eigenenergy to be found. One can re-write this as $1$st order system of equations (with $Q = dP/dr$):

\begin{align} dP/dr &= Q \\ dQ/dr &= 2(−1/r−E+1/r^2)P \\ \end{align} , with $P(0)=0$ and $Q(0)=0$.

Surprisingly, RK4 will output zero for all r > 0!

How to mend this? $P(r) \sim r^2$, for $r$ small. $E < 0$ is a parameter, usually of the order of unity.

Chip
  • 1,199
  • Well, $r \mapsto (0,0)$ is a solution, so it is not really a surprise! – copper.hat Feb 10 '16 at 06:55
  • @copper.hat we would like to see the non-trivial cases :) The Schroedinger equation reads $H \psi = E \psi $, yet this does not mean $\psi = 0$, which would render the quantum mechanics useless. – Chip Feb 10 '16 at 07:27
  • my question relates to the question: http://math.stackexchange.com/q/992327 – Chip Feb 10 '16 at 08:17
  • Not sure if it is too surprising, since with your boundary conditions of $P(0) = 0, \frac{dP}{dr}(0) = 0$ then the numerical solution will yield $P(r) = 0$. Perhaps you can start at some other place than $r=0$, with $P(x) = x, \frac{dP}{dr}(x) = 2 \sqrt{x}$ and then normalise the resulting solution? Or what about defining a function $F(r)$ such that $P(r) = r^2 F(r)$? – jim Feb 10 '16 at 09:47
  • @jim: 1) any particular reason for your recommendations $P(x)=x$ and $dP/dr(x)=2 \sqrt x$? In the article I cite below (from ArXiv), the authors do start from some small $r$, not $r=0$. 2) indeed, i can confirm that using your function $F(r)$ solves this problem (taking $F(0)$ different from $0$ and finding somehow $dF/dr(r=0)$ ). – Chip Feb 10 '16 at 10:18
  • @Chip I was thinking that, if for small $r$, $P(r) \alpha r^2$, say $P(r) = Ar^2$, and put $A = 1$, then $\frac{d P(r)}{d r} = 2 r$. If you start off the numerical integration with $P(x) = X$, where $x$ is some initial starting point and $X$ is some number, then for consistency shouldn't you initialise $\frac{d P(r)}{d r} = 2 \sqrt{X}$? I think the confusion is that I see I previously used $x$ instead of $X$, poor notation on my part. I think the following paper may be of interest to people: Cooper and Kermode Journal of Physics A-mathematical and General - vol. 19, no. 6, pp. 859-863, 1986 – jim Feb 10 '16 at 18:57
  • thanks @jim! I will have a look at the reference. – Chip Feb 11 '16 at 01:20

1 Answers1

1

$P( r ) = 0$ and $Q( r ) = 0$ are the solutions of the given equation for the initial values $P(0) = 0$ and $Q(0) = 0$. The problem is that the rate function is unbound for any solution starting anywhere but from zero (because of $1/r$ and $1/r^2$ terms). So, one cannot even invoke Cauchy theorem here, which means the problem is not well poised.

mobiuseng
  • 251
  • The $1/r$ is unbound, but $P(r)/r$ need not be so. Same goes for $1/r^2$ term. I can assure you these equations (that include the eigenvalue problem for the hydrogen atom) have very well studied and some analytical solutions. In fact, for the eigenvalue $E=-1/8$ (the energy of the 2S level of hydrogen, in atomic units of energy), the solution reads $P=r^2 exp(-r/2)$. – Chip Feb 10 '16 at 07:25
  • @Chip If $P(0)\neq 0$, the rate will be unbound. If, however, $P(0)=0$ then there is a solution $P(r)=0$. You do have a different solution $P(r)=r^2\exp(-r/2)$ starting from the same point, hence my point of the problem being not well poised as a Cauchy problem. Is it BVP, perhaps? – mobiuseng Feb 10 '16 at 07:52
  • it is an eigenvalue problem. The procedure is: for a given $E$ one propagates from zero to a small $r_{ctp}$ value (outward propagation). Secondly, starting from large $r$ (where we assume $P \approx exp(-\sqrt(-2E)r)$, one propagates backward to meet the inward solution at $r_{ctp}$, where $P$ and $Q$ should be the same. Then, one eigenvalue $E$ is determined. Hope this clarifies better the issue. Alternatively, you can say is a BVP (looking at the second order ODE, not in the form I wrote) which requires the function P to be zero at $r=0$ and some very large $r$ (order of $100$). – Chip Feb 10 '16 at 08:08
  • However, the question here is: how come RK4 predicts zero, while one can write down a Taylor series around $r=0$ that is not identically zero? – Chip Feb 10 '16 at 08:12
  • @Chip RK4 (and I suspect any IVP solvers) will predict zero as it is a valid solution: all the rates are 0 at $r=0$ thus $P(h)=0$ for small $h$, and so on. If you are trying to solve BVP, however, the case is different: RK4 won't work - you'll need something like collocations method. – mobiuseng Feb 10 '16 at 08:34
  • to prove that RK4 can solve this, albeit in a somewhat `cheating' way (probably just to avoid the $r=0$ as starting point issue I pointed out, and use instead something like $r \approx 10^{-7}$), please have a look at this article http://arxiv.org/abs/1209.1752 , Eqns. ($69$) onwards and their references in the text before. – Chip Feb 10 '16 at 09:01
  • Close to $r=0$, the ODE is stiff, thus explicit methods will fail. Use implicit methods designed for stiff problems. But even then you will have to avoid the singularity. – Lutz Lehmann Feb 10 '16 at 09:23
  • @LutzL: please see the article cited in my previous comment and you will be convinced that explicit RK4 works just fine (though they mention some other method in addition to RK4 used only as starting method). One ODE solution is given also in one of my previous comments and it is well behaved. The singularity is removable (ie, P/r and the like are finite). – Chip Feb 10 '16 at 09:38
  • Then you will need a dynamic step size obeying $h<2·r^2$. Or you need to get to $r=0.03$ using alternative means. – Lutz Lehmann Feb 10 '16 at 09:44
  • @LutzL: this should be tested numerically. Anyone interested? In the article from ArXiv (link above), you will see that there is no mentioning of stiffness in relation to this ODE. This ODE and RK4 application is puzzling. – Chip Feb 10 '16 at 10:07
  • Quote from the article: In practice, we use RK4 to compute the first 4 points needed by implicit Adams, then switch to the more accurate implicit Adams method to compute the remaining points. I guess even this is not the whole story, it would be interesting what they use for $f(0,0)$ in their RK4. – gammatester Feb 10 '16 at 10:12
  • 1
    @gammatester: indeed, I agree. From what I read, they use a change of coordinates (see Section $3.2$ in the article - "Mesh") such that $r=0$ is never reached, and instead start from very close to $0$. – Chip Feb 10 '16 at 10:22