Closed form expression for position as function of time of object falling directly into black hole from infinity Given a Schwarzschild radius $r_s=2 G M/c^2$, the escape velocity (equal to speed if falling from infinity) will be $\sqrt{2 G M/r}=\sqrt{r_s c^2/r}$
where the radial distance "r" is the point at which the measured circumference is $2 \pi r$. Lets assume that we set its speed equal to the local escape velocity when throwing our probe into the black hole from a relatively safe distance of $100 r_s$.
As you get closer to the event horizon, space is expanded in the radial direction by $1/\sqrt{1−r_s/r}$ so even as its speed gets closer to $c$ relative to an outside stationary observer, it will slow down asymptotically relative to the Schwarzschild radial distance.
What I would like to know is if there exists an exact formula for r(t) (either  for radial distance or total measured distance fallen in frame of observer with same velocity as black hole), or if you have to numerically integrate the formula for speed (divided by space expansion) as a function of radial distance.
 A: When you solve for a radial timelike geodesic, outside the horizon of a Schwarzschild black hole, for a particle falling in from rest at infinity, you can't get a nice formula for $r(t)$, but you can get one for $t(r)$:
$$t(r)-t_0=\frac{\sqrt{2}}{3}M\left[\sqrt\frac{r_0}{M}\left(6+\frac{r_0}{M}\right)-\sqrt\frac{r}{M}\left(6+\frac{r}{M}\right)\right]-4M\left[\tanh^{-1}\sqrt\frac{2M}{r_0}-\tanh^{-1}\sqrt\frac{2M}{r}\right] \tag{1}$$
in units with $c=G=1$. The integration constants have been chosen to make $r=r_0$ at $t=t_0$. The inverse hyperbolic tangent makes $t\rightarrow\infty$ as $r\rightarrow 2M$.
Addendum: The OP clarified that he is more interested in $r(\tau)$ for radial in-fall from rest at infinity, where $r$ is the radial Schwarzschild coordinate and $\tau$ is the proper time along the geodesic (i.e., the time as measured by the infalling probe). This is given by the simpler formula
$$r(\tau)=\left[\frac{9M(\tau_0-\tau)^2}{2}\right]^{1/3} \tag{2}$$
where $\tau_0$ is the proper time at which the probe reaches the singularity at $r=0$.
This is the result of integrating
$$\frac{d^2r}{d\tau^2}+\frac{M}{r^2}=0 \tag{3}$$
or
$$\frac{1}{2}\left(\frac{dr}{d\tau}\right)^2-\frac{M}{r}=0. \tag{4}$$
Remarkably, (3) and (4) are identical in form to the Newtonian versions, with the Euclidean radial coordinate replaced by the Schwarzschild radial coordinate, and the absolute time replaced by the proper time.
To derive (3), one combines
$$\frac{d^2r}{d\tau^2}+\frac{M}{r^2}\left(1-\frac{2M}{r}\right)\left(\frac{dt}{d\tau}\right)^2-\frac{M}{r^2}\left(1-\frac{2M}{r}\right)^{-1}\left(\frac{dr}{d\tau}\right)^2=0, \tag{5}$$
which is the $r$-component of the geodesic equation
$$\frac{d^2x^\mu}{d\tau^2}+\Gamma^\mu_{\alpha\beta}\frac{dx^\alpha}{d\tau}\frac{dx^\beta}{d\tau}=0 \tag{6}$$
in the case of a purely radial geodesic, with the equation
$$1=\left(1-\frac{2M}{r}\right)\left(\frac{dt}{d\tau}\right)^2-\left(1-\frac{2M}{r}\right)^{-1}\left(\frac{dr}{d\tau}\right)^2, \tag{7}$$
which comes from the expression for the proper time as given by the metric,
$$d\tau^2=\left(1-\frac{2M}{r}\right)dt^2-\left(1-\frac{2M}{r}\right)^{-1}dr^2. \tag{8}$$
Eliminating $dt/d\tau$ between the (5) and (7) gives (3).
With the $G$'s and $c$'s restored, $r(\tau)$ is
$$r(\tau)=\left[\frac{9GM(\tau_0-\tau)^2}{2}\right]^{1/3}. \tag{9}$$
The derivation of (1) is only a bit more complicated. From (4) one has
$$\left(\frac{dr}{d\tau}\right)^2=\frac{2M}{r}. \tag{10}$$
Substituting this into (7) gives
$$\left(\frac{dt}{d\tau}\right)^2=\left(1-\frac{2M}{r}\right)^{-2} \tag{11}$$
Thus
$$\frac{dr}{dt}=\frac{\frac{dr}{d\tau}}{\frac{dt}{d\tau}}=-\left(\frac{2M}{r}\right)^{1/2}\left(1-\frac{2M}{r}\right). \tag{12}$$
Integrating this gives
$$t=-\int\left(\frac{2M}{r}\right)^{-1/2}\left(1-\frac{2M}{r}\right)^{-1}dr \tag{13}$$
or
$$\frac{t}{2M}=\int u^{-5/2}(1-u)^{-1}\;du \tag{14}$$
where $u=2M/r$. This integral gives
$$\frac{t}{2M}=-2u^{-1/2}-\frac{2}{3}u^{-3/2}+\log(1+u^{1/2})-\log(1-u^{1/2})+C \tag{15}$$
or
$$\frac{t}{2M}=-2\left(\frac{r}{2M}\right)^{1/2}-\frac{2}{3}\left(\frac{r}{2M}\right)^{3/2}+\log\frac{(\frac{r}{2M})^{1/2}+1}{(\frac{r}{2M})^{1/2}-1}+C. \tag{16}$$
This is just (1) in a different form, if one uses the identity
$$\tanh^{-1}z=\frac{1}{2}\log\frac{1+z}{1-z}. \tag{17}$$
In dimensionless coordinates $r^\prime=r/2M$ and $t^\prime=t/2M$ scaled by the Schwarzschild radius $2M$, we have, dropping the primes,
$$t=-2r^{1/2}-\frac{2}{3}r^{3/2}+\log\frac{r^{1/2}+1}{r^{1/2}-1}+C. \tag{18}$$
Unfortunately, this nice expression for $t(r)$ can't be inverted to get a nice expression for $r(t)$.
