Where can I find the most direct and simplest presentation of what geodesics on a (complex) Grassmannian look like? I know how to do it from scratch, but, if I want to provide a reference to, say, a graduate student in EE who doesn't want to deal with any unnecessary abstract mathematical machinery, what should I point him to?

Grassmanians are symmetric spaces, and symmetric spaces are "geodesic orbit spaces", that is, their geodesics are orbits of their group of isometries. Your Grassmanians, in particular, are of the form $SU(p+q)/SU(p)\times SU(q)$. If $g$ is the Lie algebra of the big group and $h\subseteq g$ the Lie algebra of the subgroup, then there is a $SU(p)\times SU(q)$invariant complement $p$ to $h$ in $g$. The geodesics are the orbits of the $1$parameter subgroups of $SU(p+q)$ whose tangent vectors are in $p$. So to compute the geodesics, you need only find that complement $p$ and compute exponentials... 


This is a variation on the 1st answer, but I find it more straightforward and have explained it to EE students. Consider the map $U(n) \to Gr(k,n)$ from the unitary group to the Grassmannian by $g \mapsto gx_0$, with $x_0$ a chosen `base point'. Put the biinvariant metric on the unitary group $U(n)$. The metric on $Gr(k,n)$ is defined so that this map is a Riemannian submersion: the orthogonal complement to the fiber is linearly isometric to the base tangent space. As such, geodesics in $U(n)$ which are ORTHOGONAL TO THE FIBER project onto geodesics in the Grassmannian. And all geodesics in the Grassmannian arise this way. Now use that geodesics in the Unitary group are all of the form $g_0 exp(t \xi)$  translates of oneparameter subgroups, and work out what it means , relative to $\xi$ for the geodesic to be tangent to the fiber in the case $g_0 = Id.$ 


This answer is a little bit redundant with the other two answers given so far, but here goes anyway. It is easier to describe the real Grassmannian case. We can look at the Grassmannian $\text{Gr}(n,k)$, and suppose that $2k \le n$; if not then you can pass to the opposite Grassmannian. If $V$ and $W$ are two $k$planes in $\mathbb{R}^n$, then there a set of $k$ 2dimensional planes that are each perpendicular to each other and each intersect $V$ and $W$ in a line. Call the angles between these lines $\theta_1,\ldots,\theta_k$. Then in the connecting geodesic $V_t$ with $V_0 = W$ and $V_1 = V$, the angles are instead $t\theta_1,\ldots,t\theta_k$. This is an explicit description that is basically equivalent to Mariano's remark about invariant complements. The complex version has the same system of angles, but complexified lines, 2planes, and $k$planes. The "angle" between two lines in a 2plane can be defined as the geometric angle between their slopes plotted on the Riemann sphere. Actually these angles are twice as large as the angles in the real case in the previous paragraph, but that makes no difference. I was vague on the positions of the planes. The orthogonal projection of $V$ onto $W$ has a singular value decomposition, and the singular values are $\cos \theta_1,\ldots,\cos \theta_n$. The orthogonal projection the other way is the transpose, or Hermitian transpose in the complex case, so it has the same singular values. The corresponding singular vectors are the lines $V \cap P_k$ and $W \cap P_k$. So the actual explicit work of finding the geodesic comes from solving a singular value problem, or equivalently an eigenvalue problem. 


The complex Grassmannian SU(n)/S(U(k) * SU(nk)) being a Hermitian symmetric space enjoys the property that its geodesics (in the standard Kaehler metric) are homogeneous, i.e., generated by action of a one parameter subgroup of SU(n). In the following reference there is an explicit construction of this map in the affine coordinates. http://www.emis.de/journals/BBMS/Bulletin/bul972/berceanu1.pdf Update: Another method for the computation of the geodesics on symmetric spaces is through the solution of the radial part of the HamiltonJacobi equation. In the case of the complex Grassmannian, it depends on min(k, nk) coordinates and depends only on the restricted roots of the symmetric space and their multiplicity (see, Helgason: Groups and geometric analysis for the definitions of the radial coordinates and the radial differential operators). 


Have a look at the paper:
There are explicit formulas for geodesics and even for the geodesic distance on real Grassmannians. This ties in with Greg Kuperberg's answer. 

