How does one motivate the analytic continuation of the Riemann zeta function? I saw the functional equation and its proof for the Riemann zeta function many times, but usually the books start with, e.g. tricky change of variable of Gamma function or other seemingly unmotivated things (at least for me!).
My question is, how does one really motivate the functional equation of the zeta function? Can one see there is some hidden symmetry before finding/proving it? For example, I think $\Gamma(s+1)=s\Gamma(s)$ for $s>1$ "motivates" the analytic continuation of the Gamma function.
 A: When I speak about intuition behind the proof of the functional equation.  I am talking about proofs similiar to this one http://www.math.harvard.edu/~elkies/M259.02/zeta1.pdf
As far as I can tell, this idea was originally formulated out of necessity.  Riemann needed it for the prime number theorem. However the intuition becomes more natural if one accepts both of the following facts:
1)The Cahen Mellin Integral transforms a Dirichlet series (something hard to work with) into a Fourier series/polynomial (which we usually are more familiar with). 
2) When applying $\zeta(2s)$ to this transform, we get the Jacobi theta function $\theta(y)$ which has tons of special structure to it.
So in a sense it is the classic scheme of "I can't work with A so I will transform it to B and then transform it back."  When you transform the $\zeta(2s)$ you are getting a Fourier series (which by itself has a ton of facts about it) but it also has a really useful functional equation to it.
A: Here is one motivation. By elementary algebraic manipulation, we have
$$1 - 2 + 3 - 4 + \cdots = \frac{1}{4};$$
see, e.g. [broken link - ?; search Wikipedia for 1 - 2 + 3 - 4]
But (see here) we have
$$(1 - 2 \cdot 2)(1 + 2 + 3 + 4 + \cdots) = 1 - 2 + 3 - 4 + \cdots,$$
and therefore
$$1 + 2 + 3 + 4 + \cdots = - \frac{1}{12}.$$
Of course, to do all of this, one has to ignore all those rules that pesky analysis professors tell you about. But the last identity is just so cool that one feels compelled to try to prove it rigorously.
A: You do not try to motivate it! Even Riemann didn't see a nice argument right away. 
Riemann's first proof of the functional equation used a contour integral and led him to a yucky functional equation expressing $\zeta(1-s)$ in terms of $\zeta(s)$ multiplied by things like $\Gamma(s)$ and $\cos(\pi{s}/2)$. Only after finding this functional equation did Riemann observe that the functional equation he found could be expressed more elegantly as $Z(1-s) = Z(s)$, where $Z(s) = \pi^{-s/2}\Gamma(s/2)\zeta(s)$. Then he gave a proof of this more symmetric functional equation for $Z(s)$ using a transformation formula for the theta-function $\theta(t) = \sum_{n \in {\mathbf Z}} e^{-\pi{n^2}t},$ which is $$\theta(1/t) = \sqrt{t}\theta(t).$$ In a sense that transformation formula for $\theta(t)$ is equivalent to the functional equation for $Z(s)$. The transformation formula for $\theta(t)$ is itself a consequence of the Poisson summation formula and also reflects the fact that $\theta(t)$ is a modular form of weight 1/2. 
Instead of trying to motivate the technique of analytically continuing the Riemann zeta-function I think it's more important to emphasize what concepts are needed to prove it: Poisson summation and a connection with modular forms (for $\theta(t)$). These ideas are needed for the analytic continuation of most other Dirichlet series with Euler product which generalize $\zeta(s)$, so an awareness that the method of Riemann continues to work in other cases by suitably jazzing up Poisson summation and a link to modular forms will leave the reader/student with an appreciation for what goes into the proof.
This proof is not intuitive and I think it's a good illustration of von Neumann's comment that in mathematics we don't understand things, but rather we just get used to them.
A: Euler–Maclaurin formula allows to give the analytic continuation of the Riemann zeta function to the halfplane ${\rm Re} s>-n$. It also calculates $\zeta(s)$ at integer points $s<0$. You get Bernoulli numbers as for positive $s$. This is the reason to look for a symmetry.
A: (1) Titchmarsh points out in his book on the zeta function (section 2.3) that if you blindly apply the Poisson summation formula to the function $f(s)=|x|^s$, you get the functional equation of the Riemann zeta function immediately, and gives a reference to a paper of Mordell where this procedure is justified.
(2) However, if I had to motivate this in a (introductory graduate) class, with students knowledgeable about complex analysis, I would do as follows: 
-- We want to count prime up to $x$ (after all, this is what Riemann had in mind);
-- Writing the number of primes as a sum of some function $g(n)$ for $n$ up to $x$, it is fairly natural (though perhaps with hindsight) to use harmonic analysis to go further; it is here also rather natural to use characters of the positive reals, and hence one gets a Mellin integral on a vertical line with large enough real part ($Re(s)=\sigma>1$, say), where the dependency on $x$ is entirely in $x^s$, and the logarithmic derivative of the zeta function comes out;
-- A naive attempt to estimate from this fails, since $|x^s|=x^{\sigma}$ is worse than the trivial bound for the number of primes up to $x$; it is therefore natural to try to move the contour of integration to the left to reduce the size of $\sigma$;
-- Hence whether something can be done along these lines is immediately related to the possibility of performing analytic continuation of the (logarithmic derivative) of the zeta function;
-- In Riemann's time, I have the feeling that people like him would then just shrug and say, "Ok, let's continue $\zeta(s)$ as much as we can", and would find a way of doing so.  Most attempts succeed at also getting the functional equation, and it becomes natural to try to understand the latter. But there are other problems where the basic strategy is the same, and one gets a limited range of analytic continuation (and there is no functional equation), while still being able to get an asymptotic formula for the sum of original interest. [Indeed, in some sense, this is what happened to Riemann: to actually count primes, one needs to continue the logarithmic derivative, and if the zeros of $\zeta(s)$ are badly located, the strategy fails -- and he did not prove the Prime Number Theorem.]
A: The analytic continuation of the Riemann zeta function doesn't depend on knowing the functional equation at all. It can be done, and motivated, much more simply.
Consider the Hurwitz zeta function, $\zeta(s, q) = \sum_{n = 0}^{\infty} (q + n)^{-s}$. The Riemann zeta function is of course the $q = 1$ case of this.
Actually, it will be slightly more convenient to work with the related functions $F_s(q) = -s \sum_{n = 0}^{\infty} (q + n)^{-s - 1} = -s\zeta(s + 1, q) = \frac{d}{dq} \zeta(s, q)$. Clearly, once we know all the functions $F_s$, we can recover $\zeta(s, q)$ as $-F_{s - 1}(q)/(s - 1)$.
The series defining $F_s$ converges for $\Re(s) > 0$, so there is no difficulty interpreting things in that regime. How do we extend to further $s$?
Note firstly that $\frac{d}{dq} F_s(q) = -s F_{s + 1}(q)$. This determines $F_s$ up to an additive constant once we already know $F_{s + 1}$.
But how do we determine the additive constant? Well, since $F_s = \frac{d}{dq} \zeta(s, q)$, the average value of $F_s$ on any unit interval from $q$ to $q + 1$ is to be $\zeta(s, q + 1) - \zeta(s, q) = -q^{-s}$. Thus, we simply choose the unique additive constant that makes this true.
In this way, we can determine each $F_s$ from $F_{s + 1}$. Iterating in this way starting from the $\Re(s) > 0$ regime, we find that $F_s$ has been defined for all $s$. And thus so has the Hurwitz zeta function $\zeta(s, q)$, and thus so has the Riemann zeta function $\zeta(s)$.
It's that simple.
Note also that essentially the same idea as used here can also give us the gamma function (indeed, the second derivative of the logarithm of the gamma function is $\zeta(2, q)$).
More generally, any time we have a condition "The average value of $F$ over the range from $q$ to $q + 1$ is $-f(q)$", there is at most one solution to this with the property that a sufficiently high derivative of it vanishes in the limit as its input is translated by addition of a large natural number (technically, what I mean by an $N$-th order derivative vanishing asymptotically is that its integral over any fixed size $N$-dimensional box vanishes asymptotically as the box is translated. I.e., technically what I am saying here is actually a statement about finite differences and not infinitesimal rates of change per se). This solution exists just in case applying sufficiently many differentiations (technically, finite differences) to $f$ makes the series $\sum_{n = 0}^{\infty} f^{(m)}(q + n)$ converge, in which case the solution for $F$ is given by taking this latter function of $q$ and re-integrating back up $m - 1$ times with suitable choice of additive constants at each stage.
A: One way to motivate the analytic (or meromorphic) continuation of the Riemann zeta function
$$ \zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}, \quad \mathrm{Re} s > 1$$
is to look at the continuous analogue
$$ \frac{1}{s-1} = \int_1^\infty \frac{1}{t^s}\ dt, \quad \mathrm{Re} s > 1$$
which clearly extends meromorphically to the whole complex plane.  So one now just has to understand the analyticity properties of the residual
$$ \int_1^\infty \frac{1}{t^s}\ dt - \sum_{n=1}^\infty \frac{1}{n^s}, \quad \mathrm{Re} s > 1.$$
For instance, using the Riemann sum type quadrature
$$ \int_n^{n+1} \frac{1}{t^s}\ dt = \frac{1}{n^s} + \int_n^{n+1} \frac{1}{t^s} - \frac{1}{n^s}\ dt$$
one can write this residual as
$$ \sum_{n=1}^\infty \int_n^{n+1} \frac{1}{t^s} - \frac{1}{n^s}\ dt;$$
since $\frac{1}{t^s} - \frac{1}{n^s} = O_s( \frac{1}{n^{\mathrm{Re} s+1}} )$, it is a routine application of the Fubini and Morera theorems to establish analytic continuation of the residual to the half-plane $\mathrm{Re} s > 0$.  Similarly, by using the trapezoidal rule type quadrature
$$ \int_n^{n+1} \frac{1}{t^s}\ dt = \frac{1}{2} \frac{1}{n^s} + \frac{1}{2} \frac{1}{(n+1)^s} + \int_n^{n+1} \frac{1}{t^s} - \frac{1}{n^s} - (t-n) (\frac{1}{(n+1)^s} - \frac{1}{n^s})\ dt$$
we can write the residual as
$$ -\frac{1}{2} + \sum_{n=1}^\infty \int_n^{n+1} \frac{1}{t^s} - \frac{1}{n^s} - (t-n) (\frac{1}{(n+1)^s} - \frac{1}{n^s})\ dt.$$
From Taylor's theorem with remainder the integrand here is $O_s( \frac{1}{n^{\mathrm{Re} s + 2}} )$, so now we obtain analytic continuation to the strip $\mathrm{Re} s > -1$.  One can keep going in this fashion using the Euler-Maclaurin formula, as mentioned in Ustinov's answer, to extend the range of meromorphic continuation to the rest of the complex plane.
Ultimately, meromorphic continuation in this case is a reflection of the natural numbers being so evenly spaced asymptotically that one can estimate sums over the natural numbers with reasonable accuracy in terms of sums over the half-line $[1,+\infty)$, where the error terms can be made as convergent as one wishes.  One can also use the Poisson summation formula to compare the sum and integral, which leads into the more traditional proof of meromorphic continuation based on theta functions etc..
