Geodesics on spheres are great circles How does one prove that on $S^n$ (with the standard connection) any geodesic between two fixed points is part of a great circle?
For the special case of $S^2$ I tried an naive approach of just writing down the geodesic equations (by writing the Euler-Lagrange equations of the length function) and solving them to gain some insights but even if the equations are solvable I can't see how to show that they are great circles. (the solutions are some pretty complicated functions which don't give me much insight)
I checked the article on Great Circles on Wolfram Mathworld for a coordinate geometry approach to it but that article looked quite cryptic to me! 
One knows that on compact semi-simple lie groups any one-parameter subgroup generates a geodesic and $S^n$ is the quotient of 2 compact semi-simple lie groups $SO(n+1)/SO(n)$. Is this line of thought useful for this question?  
=================================================================================
After some of the responses came let me put in "a" way of seeing the above for $S^2$ (wonder if it is correct). If $\theta$ and $\phi$ are the standard coordinates on $S^2$ then the equations for the curve are
$$\ddot{\theta} = \dot{\phi}^2 sin(\theta)cos(\theta)$$
$$\dot{\phi}sin^2{\theta} = k$$ 
where $k$ is some constant set by the initial data of the curve. 
Now given the initial point I can choose my coordinate system such that the the initial data looks like $\dot{\phi}=0$, $\theta = \text{some constant}$, $\dot{\theta}=\text{some constant}$, $\phi = \text{some constant}$. Then the differential equations tell me that the $k=0$ and the only way it can happen for times is by having  , 
$$\dot{\phi} = 0$$ 
Which clearly gives me a longitude in this coordinate system. 
Hence the geodesic equation gives as a solution a great circle.
Surely not an elegant proof like Bar's reference. 
But I hope this is correct. 
{As a friend of mine pointed out that this set of coordinates is motivated by the fact that the way the "energy" of the curve is being parametrized the z-component of the angular momentum is conserved which is in fact my second Euler-Lagrange equations} 
 A: For $S^2$, I would consider taking two points on the equator, and seeing if that is easier to solve. Since you can use isometries to move two arbitrary points onto the equator, that should prove it for $S^2$. For $S^n$ the same should also work.
It's the same trick that is used to find the geodesics in the upper half-plane with the hyperbolic metric, really. You compute the geodesics between two points lying on a vertical line, then use isometries to find the rest of them.
A: It is enough to show that the geodesics through a specific point are all of that form, and we can do this just for the north pole $N=(1,0,\dots,0)$, as the isometry group of the sphere is transitive. Moreover, we need only consider one unit tangent vector at $N$, for the stabilizer of $N$ in the isometry group of the sphere acts transitively on the unit vectors in the tangent space at $N$. So let $v=(0,1,0,\dots,0)$ be the initial speed of a geodesic $\gamma$ starting at $N$.
Since the map $(x_1,\dots,x_n)\mapsto (x_1,x_2,-x_3,\dots,-x_n)$ preserves both the point and the vector, the geodesic $\gamma$ must be also preserved by it. It follows immediately that the curve is contained in a great circle.
A: Actually it is pretty easy even to solve the Euler-Lagrange equations in the case of S^n. Please check example 7.3 in the following reference, where the sphere S^n is parameterized by unit vectors in R^(n+1).
A: Although Jose has made the essentially the same point, I just want to elaborate (this is really just a comment, but I always run out of room in the comment box).
What nobody else has mentioned explicitly is that you should have trouble solving the Euler-Lagrange equation for the length functional. The Euler-Lagrange equation is a second order ODE, but a highly degenerate one. And you know this before you even start. Why? Well, suppose you have a solution. Then if you reparameterize that curve using any arbitrary parameterization (i.e., any monotone function of the original parameter), the newly parameterized curve is still a solution to the Euler-Lagrange equation. That means that the ODE has an infinite dimensional space of solutions and it is nothing like any ODE we learned about in our ODE courses or textbooks.
A trick is needed to get around this, namely to use the so-called energy functional $E[\gamma] = \int_0^1 |\gamma'(t)|^2 dt$ (which is not invariant under reparameterization of the curve) instead of the length functional (which is). The Holder inequality shows that a minimum of the energy functional is necessarily a minimum of the length functional that is parameterized by a constant times arclength, i.e. a constant speed geodesic.
The Euler-Lagrange equation for the energy functional is a nice nondegenerate 2nd order ODE that can be handled by standard ODE techniques and theorems.
As for the standard sphere, there are many different ways to solve for the geodesics. To review the ways already suggested in other answers:
1) I recommend that you first do it without the machinery of Riemannian geometry and using only the Euclidean structure of $R^{n+1}$. Using the discussion above, you should be able to show that a curve on the unit sphere is a constant speed geodesic if and only if its acceleration vector is always normal to the sphere. You should then be able to work out the solutions to this ODE. The suggestion that you assume one point is the north pole and the other lies in a co-ordinate plane is a good one and makes the ODE easy to solve.
2) The other way is to do it all intrinsically. Here, I recommend using stereographic co-ordinates and assuming one point is the origin in those co-ordinates. Again, everything becomes very easy in that situation.
3) And the third way is to view the sphere as a homogeneous space and use formulas for that situation. I don't remember the details myself, but I learned them from the book by Cheeger and Ebin.
I recommend that you work through all 3 different ways, as well as any other way you can find.
As others have noted, the calculations for geodesics on hyperbolic space are identical, except that you are working with a "unit sphere" in Minkowski instead of Euclidean space. There is even a notion of stereographic projection (but onto what?). This is also fun to work out carefully.
Finally, I do want to note that after you work this all out and have it all in your head, it's a really beautiful picture and story. And if you find the right angle, it's all very simple, so you can work out the details yourself and not rely on reading a book line-by-line or having someone else show you all the details. Try to get the essential ideas and necessary tricks (like using the energy functional) from books, lectures, or teachers, but try to work everything else out from scratch (i.e., minimal reliance on theorems you can't prove yourself).
A: A simple argument that shows that great circles are geodesics is that if you parametrise them in such a way that they have unit speed, their acceleration is normal to the surface.  To show that all geodesics are great circles, just use uniqueness of the initial value problem for a geodesic after noticing that through every point and in any direction there is a great circle.
A: Here is a much simpler proof. We know that the great circles of a sphere $S^{2}$ are geodesics. Let $p$ and $q$ be two points on $S^{2}$. Now find a plane that contains the center of the sphere, $p$ and $q$. The intersection of the plane and the sphere is a great circle with $p$  and $q$ being points on the great circle. Hence, the geodesic joining $p$ and $q$ is part of the great circle.
A: *

*Notice that a geodesic (between points p, q) must be a subset of a plane - to see this, if the geodesic is 3 dimensional, then, looking at any plane through p, q, the projection of the path onto it can be reparametrised to become a valid path on the sphere with less length.

*Any plane through p, -p must contain the center of the sphere. Hence it cuts a great circle on the sphere.

*Let p, q be two points. If the plane that contains the geodesic between the two points does not go through the center, then taking the great segment between q,-p and the path through plane not containing center between p, q, we get a smaller length path between p, -p that isn't confined to a single plane, contradiction to 1. Hence, the plane goes through center, and the path between p and q is a great circle segment.

A: Here is another argument.
Given a great circle in $S^n$, there is an isometry of the sphere $S^n$ whose fixed points are precisely that great circle, and on the normal bundle (fibers) the symmetry is the antipodal map.
So whatever the covariant 2nd derivative of the great circle is, the component normal to the great circle is invariant under this antipodal map.  So it must be zero, and the great circle is a geodesic.
This argument works equally well for geodesics in $\Bbb R^n$ or in hyperbolic space.
