The probability density of the ground state is time independent, so there is no motion in this sense. Yet the expectation value of the kinetic energy is non-zero, so there is movement in this sense. How are these notions of movement reconciled?
First off, classically, if we had a particle in a $1/r$ potential and released it from rest, it would indeed bob back and forth like a pendulum as you describe. But in quantum mechanics we can't say the electron is taking any specific path around the proton. As there is no specific path, we cannot fully reconcile these notions of movement with any classical preconceptions.
Let's discuss some various notions of motion in quantum mechanics, that may help you here.
The ground state of the hydrogen atom is
$$ \psi(r,\theta,\phi,t) = A\, \exp\left(-\frac{r}{a}-i\frac{E t}{\hbar}\right)$$
Where $A =\sqrt{\frac{1}{\pi a^3}},a=\frac{\hbar^2}{me^2},E=\frac{-m e^4}{8 h^2 \varepsilon_0^2}$
The radial momentum operator in this basis is:
$$\vec{p}_r = -i\hbar\hat{r}\frac{\partial}{\partial r}$$
where $\hat{r}$ is the radial unit vector (not an operator).
In calculating the expectation value of this:
$$\langle \psi|\vec{p}_r|\psi\rangle
= \int \psi^* (- i\hbar\hat{r}\frac{\partial}{\partial r} \psi) \sin(\theta) r^2 dr\,d\theta\,d\phi
= \int \psi^* (- i\hbar\hat{r}\frac{-1}{a} \psi) \sin(\theta) r^2 dr\,d\theta\,d\phi
$$
Due to symmetry, this will of course be zero. But the density term in the integral is
$$\psi^*(- i\hbar\hat{r}\frac{-1}{a} \psi) = i\hbar\frac{1}{a}(\psi^*\psi)\hat{r}$$
This might be what you want to interpret as 'motion', but since $(\psi^*\psi)\ge 0$ this is purely imaginary and doesn't have a directly physical interpretation as motion. As its imaginary, it's neither towards or away from the center.
Another notion of motion is the probability current:
$$ \vec{j} = \frac{\hbar}{2mi}\left(\psi^* \vec\nabla \psi - \psi \vec\nabla \psi^{*} \right)$$
This is related to conservation of probability by:
$$ \rho = \left|\psi\right|^2,\quad \frac{\partial \rho}{\partial t} + \vec\nabla \cdot \vec{j} = 0 $$
For the hydrogen ground state we have:
$$ \vec{j} = \frac{\hbar}{2mi}\hat{r}\left(\psi^* \left(\frac{\partial}{\partial r}\psi\right) - \psi \left(\frac{\partial}{\partial r}\psi^{*}\right) \right)$$
$$ = \frac{\hbar}{2mi}\hat{r}\left(\psi^* \left(\frac{-1}{a}\psi\right) - \psi \left(\frac{-1}{a}\psi^{*}\right) \right) = 0$$
There is no probability current at any point. So any sense in which there is motion at some location, the net current into/out of this point is still zero. Which takes me to the only remaining way I know to discuss "motion" here. We are writing the state in the position basis, let me make this more clear, and also use the cartesian basis for a bit:
$$ |\psi\rangle = \int \phi(x,y,z,t)|x,y,z\rangle $$
The state $|\phi\rangle$ is but one vector in the infinite vector space, that is the Hilbert space for the electron here. When we write $\phi(x,y,z,t)$ these are really time dependent components for each basis element $|x,y,x\rangle$ in the chosen basis for this vector space. The state $|1,0,0\rangle$ by itself is the closest interpretation of your idea of starting with the electron at say x=1,y=0,z=0 and dropping it to see how it moves.
We can start with this pure position state and watch how it evolves according to the Hamiltonian operator. Since this state is not an energy eigenstate, it will spread (evolve to a state that now needs to be written as a super-position of many of our $|x,y,x\rangle$ basis states). It however will not swing like a pendulum though the origin like you imagine. It will spread out in all directions (since by the uncertainty principle, a pure position state is completely spread in momentum space).
The magic of the ground state is that if we consider this special weighted super position of an enormous (infinite) number of position states individually spreading, they spread exactly such that the super position of states remains the same and the net current is zero at every point. You could view this a bit like equilibrium with the principle of detailed balance: the position states will evolve into each other, but the amount that is "leaving" a pure position state must be replaced with exactly the same amount "entering" that state from other position states in this super position.
So in a sense, there is movement (kinetic energy is non-zero, the time evolution operator (Hamiltonian) is constantly evolving pure position states at each point to spread out), but the "net movement" of the wavefunction is zero (probability current is zero) and the probability density is time independent.
Consider this section an extended comment:
Akrasia suggested another way of looking at motion here: momentum decomposition.
Basically we can also write the state in terms of momentum basis in Hilbert space.
$$ |\psi\rangle = \int \phi(k_x,k_y,k_z,t)|k_x,k_y,k_z\rangle $$
These basis states are spread over all space (uniformly so). So they can't tell us about motion in some region. But we can get a probability density on this space, giving a notion of motion for parts of the state. And for the hydrogen ground state, it will be built up as standing waves of opposing momentum basis states. Since these cover all space, the momentum of a plain wave state is not just in the radial direction. So in this sense, the "motion" is not just in the radial direction.
Here is a related stack exchange question:
Hydrogen wave function in momentum space
And here is a paper which claims to work out the hydrogen wave function in spherical momentum space.
They find the ground state to be:
$$\phi_{1,0,0} = \frac{1}{(p_r - i p_0)^2} \frac{1}{p_\theta^{1/2}}J_{1/2}(p_\theta)\delta(p_\phi)$$
Which means there is no contribution from basis elements with non-zero $p_\phi$ momentum, but there is for non-zero $p_\theta$. That is surprising to me, and I don't have time to read through the paper right now. So it would be best if someone else wrote up an answer covering this portion. If that paper works, then that seems to be a very nice notion of motion for which to answer the question here "Is there only radial motion in the Hydrogen ground state?".
A nice overview of the problem is given in arXiv:1205.3740. I'll summarise the most important points here.
Let $d$ be the number of space dimensions. Then the Laplace operator is given by
$$
\Delta=\frac{\partial^2}{\partial r^2}+\frac{d-1}{r}\frac{\partial}{\partial r}+\frac{1}{r^2}\Delta_S\tag{1}
$$
where $\Delta_S$ is the Laplace operator on the $d-1$ sphere.
The Coulomb potential is given by the solution of
$$
-\Delta V=e\delta(\boldsymbol r)\tag{2}
$$
which is solved by
$$
V(r)=\frac{2(d/2-1)!}{(d-2)\pi^{(d-2)/2}}\frac{e}{r^{d-2}}\tag{3}
$$
One way to show the expression above is to consider the Gauss Law in $d$ dimensions, that is, $E(r)=e/\int_r\mathrm ds$ where the area of the $d-1$ sphere can be found here.
With this, the Schrödinger equation reads
$$
\left[\frac{1}{2m}(-\Delta)+V(r)-E\right]\psi(\boldsymbol r)=0\tag{4}
$$
The problem has spherical symmetry so that, as usual, we may write
$$
\psi(r,\boldsymbol \theta)=\frac{1}{r^{(d-1)/2}}u(r)Y_\ell(\boldsymbol\theta)\tag{5}
$$
where $\boldsymbol \theta=(\theta,\phi_1,\phi_2,\cdots,\phi_{d-2})$ and the power $(d-1)/2$ of $r$ is chosen to remove the linear term in $(1)$. The superspherical harmonics ($\sim$ Gegenbauer polynomials) are the generalisation of the usual spherical harmonics to $d$ dimensions:
$$
\Delta_S Y_\ell+\ell(\ell+d-2)Y_\ell=0\tag{6}
$$
Using this form for $\psi$, the Schrödinger equation becomes
$$
u''(r)+2m[E-V_\ell(r)]u(r)=0\tag{7}
$$
where the effective potential is
$$
V_\ell=V(r)+\frac{1}{2m}\frac{\ell_d(\ell_d+1)}{r^2}\tag{8}
$$
with $\ell_d=\ell+(d-3)/2$.
Now, this eigenvalue problem has no known analytical solution, so that we must resort to numerical methods. You can find the numerical values of the energies on the arxiv article. An important point that is addressed in that article is that there are no negative eigenvalues for $d\ge 4$, that is, there are no bound states in more than three dimensions. But for $d\ge 5$ there are stable orbits, with positive energy, with well behaved wave functions.
To answer to some of your questions:
In general you need $d$ quantum numbers for $d$ dimensions, modulo degeneracies. In the case of the hydrogen atom in 3D, there are the spherical symmetry and the accidental symmetry$^1$, so that you only need one quantum number. In $d$ dimensions the spherical symmetry remains, but I think that the accidental does not, which would mean that you need $d-1$ quantum numbers for $d\neq 3$. One would need to check whether the Runge-Lenz vector is conserved in $d$ dimensions or not (left as an exercise). If this vector is actually conserved, then the energies would depend on $d-2$ quantum numbers.
As there is no analytical solution to the radial Schrödinger equation, we don't know. In the case of $d=3$ dimensions, the Bohr-Sommerfeld quantisation rule turns out to be exact. We could check what this scheme predicts for $E_n$ (though we couldn't know whether it would be exact or not: we must compare to the numerical results).
These are well-known for mathematicians. You can read about them on the wikipedia article.
In closed form, no. I don't know of an asymptotic formula, but it should be rather easy to derive from the radial Schrödinger equation, where the centrifugal term $r^{-2}$ dominates for $r\to\infty$, so that we can neglect the Coulomb term.
$^1$ In spherical symmetric systems, the potential depends on neither $\theta$ nor $\phi$. In these systems energies don't depend on $m$, the azimuthal quantum number. Therefore, in general for radial potentials the energies depend on two quantum numbers, $n,\ell$. In the specific case of $V=1/r$, there is another symmetry that is kind of unexpected (or at least it is not very intuitive geometrically). When $V=1/r$ the rotational symmetry $SO(3)$ is enlarged into a $SO(4)$ symmetry, and this new symmetry is known as accidental symmetry. This symmetry makes the energy levels independent of $\ell$, that is, $E=E_n$. Note that this symmetry is not present in the rest of the atomic table, that is, multielectronic atoms, which makes the energy levels depend on the angular momentum (and thus, Hund's rules).
One may illustrate the above as
\begin{aligned}
\text{If $V=V(r,\theta,\phi)$} &\longrightarrow E=E_{n\ell m}\\
\text{If $V=V(r)$}&\longrightarrow E=E_{n\ell}\\
\text{If $V=1/r$}&\longrightarrow E=E_{n}
\end{aligned}
The first line is the general result for 3D systems; the second line is the result of spherical symmetry; and the third line is the result of the accidental symmetry. If you want to read more about this symmetry, you'll find some nice references in Why are hydrogen energy levels degenerate in $\ell$ and $m$? and/or here.
Best Answer
The infinitesimal probability for the electron to be in the volume $dV$ around a point $(r,\theta,\phi)\leftrightarrow (x,y,z)$ is given by $$ dP = dV\cdot |\psi(x,y,z)|^2 = dV\cdot |R(r)|^2\cdot |Y_{lm}(\theta,\phi)|^2 =\dots$$ as you can see if you substitute your Ansatz for the wave function. However, the infinitesimal volume $dV=dx\cdot dy\cdot dz$ may be rewritten in terms of differentials of the spherical coordinates as $$ dV = dr\cdot r^2 \cdot d\Omega = dr\cdot r^2 \cdot \sin\theta\cdot d\theta\cdot d\phi $$ where the small solid angle $d\Omega$ was rewritten in terms of the spherical coordinates. You see that for dimensional reasons (or because the surface of a sphere scales like $r^2$), there is an extra factor of $r^2$ in $dV$ and therefore also in $dP$ which suppresses the probability. There is simply not enough volume for small values of $r$.
So $|R(r)|^2$ may still go like $1/r^2$ for small $r$ and in that case, $dV$ will be proportional to $dr$ times a function that is finite for $r\to 0$. Such $dP$ may be integrated and there's no divergence at all near $r=0$.
That's why one should allow the wave function to go like $1/r$ near $r=0$ which is the true counterpart of one-dimensional wave function's being finite near a point. However, Nature doesn't use this particular loophole because the wave function $\psi$ for small $r$ actually scales like $r^l$ where $l$ is the orbital quantum number and the wave function actually never diverges even though it could.
Update 2016: I should have and could have written it four years ago but I didn't. While the normalizability allows $1/r$ around $r=0$, such singular functions ultimately can't be in stationary or nearly stationary states, for the following reason which differs from various reasons above and those in the comments.
For example, someone mentioned that $1/r$ could lead to a continuous spectrum or some surprising degeneracies. But if the correct wave functions predicted a continuous or degenerate spectrum in a box, then it would be how Nature works. The actual reason why $1/r$ is not finally allowed as a stationary wave function near $r=0$ is that the Laplacian of this wave function (and Schrödinger's equation contains such a Laplacian) is proportional to a delta-function at the origin (or contains such a term) and no other term in Schrödinger's equation can cancel this delta-function, so the Schrödinger equation must be violated.