Updated by Michael Weiss, 2017.
Original by Michael Weiss and John Baez.
In special cases, yes. In general, it depends on what you mean by "energy", and what you mean by "conserved".
In flat spacetime (the backdrop for special relativity), you can phrase energy conservation in two ways: as a differential equation, or as an equation involving integrals (gory details below). The two formulations are mathematically equivalent. But when you try to generalize this to curved spacetimes (the arena for general relativity), this equivalence breaks down. The differential form extends with nary a hiccup; not so the integral form.
The differential form says, loosely speaking, that no energy is created in any infinitesimal piece of spacetime. The integral form says the same for a non-infinitesimal piece. (This may remind you of the "divergence" and "flux" forms of Gauss's law in electrostatics, or the equation of continuity in fluid dynamics. Hold on to that thought!)
An infinitesimal piece of spacetime "looks flat", meaning that even a non-zero curvature is hard to detect on an infinitesimal piece of spacetime; but the effects of curvature do become more evident in a non-infinitesimal piece. (The same holds for curved surfaces in space, of course.) In GR, the curvature of spacetime is "felt" by us as gravity. Now, even in newtonian physics, you must include gravitational potential energy to get energy conservation. And GR introduces the new phenomenon of gravitational waves; perhaps these carry energy as well? Perhaps we need to include gravitational energy in some fashion, to arrive at a law of energy conservation for non-infinitesimal pieces of spacetime?
Casting about for a mathematical expression of these ideas, physicists came up with something called an energy pseudo-tensor. (In fact, several of 'em!) Now, GR takes pride in treating all coordinate systems equally. Mathematicians invented tensors precisely to meet this sort of demand: if a tensor equation holds in one coordinate system, it holds in all. Pseudo-tensors are not tensors (surprise!), and this alone raises eyebrows in some circles. In GR, one must always guard against mistaking artifacts of a particular coordinate system for real physical effects. (See the FAQ entry What is a black hole, really? for some examples.)
These pseudo-tensors have some rather strange properties. If you choose the "wrong" coordinates, they are non-zero even in flat empty spacetime. By another choice of coordinates, they can be made zero at any chosen point, even in a spacetime full of gravitational radiation. For these reasons, most physicists who work in general relativity do not believe the pseudo-tensors give a good local definition of energy density, although their integrals are sometimes useful as a measure of total energy.
One other complaint about the pseudo-tensors deserves mention. Einstein argued that all energy has mass, and all mass acts gravitationally. Does "gravitational energy" itself act as a source of gravity? Now, the Einstein field equations are $$ G_{\mu\nu} = 8\pi T_{\mu\nu} \;. $$ Here $G_{\mu\nu}$ is the Einstein curvature tensor, which encodes information about the curvature of spacetime, and $T_{\mu\nu}$ is the so-called stress–energy tensor, which we will meet again below. $T_{\mu\nu}$ represents the energy due to matter and electromagnetic fields, but includes no contribution from "gravitational energy". So one can argue that "gravitational energy" does not act as a source of gravity. On the other hand, the Einstein field equations are non-linear; this implies that gravitational waves interact with each other (unlike light waves in Maxwell's (linear) theory). So one can argue that "gravitational energy" is a source of gravity.
In certain special cases, energy conservation works out with fewer caveats. The two main examples are static spacetimes and asymptotically flat spacetimes.
Let's look at three examples before plunging deeper into the mathematics.
Light from the Sun appears redshifted to an Earth bound astronomer. In quasi-newtonian terms, we might say that light loses kinetic energy as it climbs out of the gravitational well of the Sun, but gains potential energy. General relativity looks at it differently. In GR, gravity is described not by a potential but by the metric of spacetime. No problem: the Schwarzschild metric describes spacetime around a massive object, if the object is spherically symmetrical, uncharged, and "alone in the universe". The Schwarzschild metric is both static and asymptotically flat, and energy conservation holds without major pitfalls. For more details, consult MTW, chapter 25.
A binary pulsar emits gravitational waves, according to GR, and one expects (innocent word!) that these waves will carry away energy. So its orbital period should change. Einstein derived a formula for the rate of change (known as the quadrapole formula), and in the centenary of Einstein's birth, Russell Hulse and Joseph Taylor reported that the binary pulsar PSR1913+16 bore out Einstein's predictions within a few percent. Hulse and Taylor were awarded the Nobel prize in 1993.
Despite this success, Einstein's formula remained controversial for many years, partly because of the subtleties surrounding energy conservation in GR. The need to understand this situation better has kept GR theoreticians busy over the last few years. Einstein's formula now seems well-established, both theoretically and observationally.
The Cosmic Background Radiation (CBR) has red-shifted over thousands of millions of years. Each photon gets redder and redder. What happens to this energy? Cosmologists model the expanding universe with Friedmann-Robertson-Walker (FRW) spacetimes. (The familiar "expanding balloon speckled with galaxies" belongs to this class of models.) The FRW spacetimes are neither static nor asymptotically flat. Those who harbor no qualms about pseudo-tensors will say that radiant energy becomes gravitational energy. Others will say that the energy is simply lost.
It's time to look at mathematical fine points. There are many to choose from! The definition of asymptotically flat, for example, calls for some care (see Stewart); one worries about "boundary conditions at infinity". (In fact, both spatial infinity and "null infinity" clamor for attention—leading to different kinds of total energy.) The static case has close connections with Noether's theorem (see Goldstein or Arnold). If the catch-phrase "time translation symmetry implies conservation of energy" rings a bell (perhaps from quantum mechanics), then you're on the right track. (Check out "Killing vector" in the index of MTW, Wald, or Sachs and Wu.)
But two issues call for more discussion. Why does the equivalence between the two forms of energy conservation break down? How do the pseudo-tensors slide around this difficulty?
First point: in relativity, energy is not a scalar, but the time component of the energy–momentum 4-vector. If you already know what that means, skip the next paragraph.
In ordinary high-school analytic geometry, a (two-dimensional) vector $\boldsymbol{v}$ has components, say $(v_1,v_2)$. While $\boldsymbol{v}$ is a geometrical object, existing free from the confines of any coordinate system, the same is not true for the components $v_1$ and $v_2$—they depend on the choice of the axes and their scales. If you change the coordinate system—say you rotate the axes—then the components change according to a standard formula. We say that $\boldsymbol{v}$ has an invariant meaning, while the components are coordinate dependent. Now in relativity, neither energy nor momentum by themselves have invariant meaning, just like time and space. But if we weld the energy and momentum together, we get a geometric object called a 4-vector that is invariant. That is, if $E$ is the energy and $\boldsymbol{p}=(p_1,p_2,p_3)$ is the momentum, then $(E,p_1,p_2,p_3)$ is the energy–momentum 4-vector. The energy–momentum 4-vector basks in celebrity, being the second most famous 4-vector; the top spot is held by the spacetime-displacement 4-vector, $(\Delta t,\Delta x,\Delta y,\Delta z)$. That's why we say that the energy is the time component of the energy–momentum 4-vector: it occupies the same slot that time does in the spacetime-displacement 4-vector.
If we want to conserve something, it must have an invariant meaning. The energy–momentum 4-vector $p$ fits the bill. Let's consider first the case of flat Minkowski spacetime. Recall from your favorite text on relativity that that the notion of "inertial frame" corresponds to a special kind of coordinate system (Minkowski coordinates, a.k.a. Lorentz coordinates), in which the metric takes the form $-\mathrm{d}t^2+\mathrm{d}x^2+\mathrm{d}y^2+\mathrm{d}z^2$.
Pick an inertial reference frame. Pick a volume $V$ in this frame, and pick two times $t=t_0$ and
$t=t_1$. One formulation of energy–momentum conservation says that the energy–momentum inside $V$
changes only because of energy–momentum flowing across the boundary surface (call it $S$). It is
"conceptually difficult, mathematically easy" to define a quantity $T$ so that the captions on equation
\eqref{equation-1} below are correct. (The quoted phrase comes from Sachs and Wu.)
Valid in flat Minkowski spacetime, when Minkowskian coordinates are used (where $p$ is the energy–momentum
4-vector),
\begin{equation}\label{equation-1}
\begin{array}{ccccl}
\int_{V,\,t\,=\,t_0} T\,\mathrm{d}V & - & \int_{V,\,t\,=\,t_1} T\,\mathrm{d}V & = & \int_{t_0}^{t_1} T\,\mathrm{d}t\,\mathrm{d}S
\;.\\
\\
p\text{ in vol }V & - & p\text{ in vol }V & = & p\text{ flowing out through }\\
\text{at time }t_0 & &\text{at time }t_1 & & \text{boundary }S\text{ of }V\\
& & & &\text{from }t_0\text{ to }t_1\;.
\end{array}
\end{equation}
$T$ is called the stress–energy tensor. You don't need to know what that means!—just that you can
integrate $T$, as shown, to get 4-vectors. Equation \eqref{equation-1} may remind you of Gauss's theorem, which
deals with flux across a boundary. If you look at \eqref{equation-1} in the right 4-dimensional frame of mind,
you'll discover it really says that the flux across the boundary of a certain 4-dimensional hypervolume is zero.
(The hypervolume is swept out by $V$ during the interval $t=t_0$ to $t=t_1$.) Misner, Thorne, and Wheeler,
chapter 7, explains this with pictures galore. (See also Wheeler.)
A 4-dimensional analogue to Gauss's theorem shows that \eqref{equation-1} is equivalent to the following:
Valid in flat Minkowski spacetime, with Minkowskian coordinates,
\begin{equation}\label{equation-2}
\nabla_\text{coord}\cdot T = \sum_\mu {\partial T\over\partial x_\mu} = 0 \;.
\end{equation}
We write $\nabla_\text{coord}\cdot T$ for the divergence, for we will meet another divergence in a moment.
Proof? Quite similar to Gauss's theorem: if the divergence is zero throughout the hypervolume, then the flux
across the boundary must also be zero. On the other hand, the flux out of an infinitesimally small
hypervolume turns out to be the divergence times the measure of the hypervolume.
Pass now to the general case of any spacetime satisfying Einstein's field equation. It is easy to
generalize the differential form of energy–momentum conservation, \eqref{equation-2}:
Valid in any GR spacetime,
\begin{equation}\label{equation-3}
\nabla_\text{cov}\cdot T = \sum_\mu \nabla_\mu T= 0 \;.
\end{equation}
$\nabla_\mu$ is the covariant derivative, and $\nabla_\text{cov}\cdot T$ is the covariant divergence.
[Side comment: \eqref{equation-3} is the correct generalization of \eqref{equation-1} for SR when non-Minkowski
coordinates are used.]
GR relies heavily on the covariant derivative, because the covariant derivative of a tensor is a tensor, and, as we've seen, GR loves tensors. Equation \eqref{equation-3} follows from Einstein's field equation (because something called Bianchi's identity says that $\nabla_\text{cov}\cdot G=0$). But \eqref{equation-3} is no longer equivalent to \eqref{equation-1}!
Why not? Well, the familiar form of Gauss's theorem (from electrostatics) holds for any spacetime, because essentially you are summing fluxes over a partition of the volume into infinitesimally small pieces. The sum over the faces of one infinitesimal piece is a divergence. But the total contribution from an interior face is zero, since what flows out of one piece flows into its neighbor. So the integral of the divergence over the volume equals the flux through the boundary. "QED".
But for the equivalence of \eqref{equation-1} and \eqref{equation-3}, we would need an extension of Gauss's theorem. Now the flux through a face is not a scalar, but a vector (the flux of energy–momentum through the face). The argument just sketched involves adding these vectors, which are defined at different points in spacetime. Such "remote vector comparison" runs into trouble precisely for curved spacetimes.
The mathematician Levi-Civita invented the standard solution to this problem, and dubbed it "parallel transport". It's easy to picture parallel transport: just move the vector along a path, keeping its direction "as constant as possible". (Naturally, some non-trivial mathematics lurks behind the phrase in quotation marks. But even pop-science expositions of GR do a good job explaining parallel transport.) The parallel transport of a vector depends on the transportation path; for the canonical example, imagine parallel transporting a vector on a sphere. But parallel transportation over an "infinitesimal distance" suffers no such ambiguity. (It's not hard to see the connection with curvature.)
To compute a divergence, we need to compare quantities (here vectors) on opposite faces. Using parallel transport for this leads to the covariant divergence. This is well defined, because we're dealing with an infinitesimal hypervolume. But to add up fluxes all over a non-infinitesimal hypervolume (as in the contemplated extension of Gauss's theorem) runs smack into the dependence on transportation path. So the flux integral is not well defined, and we have no analogue for Gauss's theorem.
One way to get round this is to pick one coordinate system, and transport vectors so their components stay constant. Partial derivatives replace covariant derivatives, and Gauss's theorem is restored. The energy pseudo-tensors take this approach (at least some of them do). If you can wrangle \eqref{equation-3}, $\nabla_\text{cov}\cdot T= 0$, into the form: $$ \nabla_\text{coord}\cdot\Theta = 0 \;, $$ then you can get an "energy conservation law" in integral form. Einstein was the first to do this; Dirac, Landau and Lifshitz, and Weinberg all came up with variations on this theme. We've said enough already on the pros and cons of this approach.
We will not delve into definitions of energy in general relativity such as the hamiltonian (amusingly, the energy of a closed universe always works out to be zero according to this definition), various kinds of energy one hopes to obtain by "deparametrizing" Einstein's equations, or "quasilocal energy". There's quite a bit to say about this sort of thing! Indeed, the issue of energy in general relativity has a lot to do with the notorious "problem of time" in quantum gravity... but that's another can of worms.
References (vaguely in order of difficulty):