A formula for the power sums: $1^n+2^n+\dotsc +k^n=\,$? Is there explicit formula for the expression $1^n + 2^n + \dotsc + k^n\,$?
I know that for $n=1$ the explicit formula becomes $S=k(k+1)/2$ and for $n=3$ the formula becomes $S^2$. But what about general $n$?
I know there is a way using the Taylor expansion of $f(x)=1/(1-x)=1+x+x^2+\dotsc\;$, by differentiating it and then multiplying by $x$ and then differentiating again. Repeating this $n$ times, we get
$$\frac{d}{dx}(x\frac{d}{dx}(\dots x\frac{d}{dx}f(x))\dots )=1+2^nx^n+3^nx^n\dots.$$
Now do the same process but with the function $g(x)=x^{k+1}f(x)$. Then subtract them and we get $1+2^nx^n+\dots k^nx^n$. Because we have the explicit formulas $f(x)$ and $g(x)$ we can find the explicit formula by this process for arbitrary $n$. A big problem is that as $n$ grows, it is going take a lot of time finding the explicit formula. My question is therefore: are there other ways?
 A: A method which is more seldom used is that involving the Eulerian numbers. I found this solution myself by completely elementary means and "pattern-detection" only- so I liked it very much and I've made a small treatize about this. 
Unfortunately it is only in German, and since it is over 12 years old I don't want to translate it just now. However, the path to go and the development of the formulae go step by step, so it should be understandable/tractable and self-explanatory enough to see, how this can be done.      
In effect, one finds a solution of the form
$$ P(m,n) = 1+2^m+3^m+...+n^m = e_{m,1}\binom{n+1}{m+1}+ e_{m,2}\binom{n+2}{m+1}+ ... +e_{m,m} \binom{n+m}{m+1} $$
with the Eulerian triangle $\{e_{r,c}\}_{r,c=0..\infty} $
$$ \begin{array} {} 1 \\ .&1 \\ .&1&1 \\ .&1&4&1 \\ .&1&11&11&1\\ ... \end{array} $$
A: We know that $$\sum_{i=1}^ki^{n+1}-(i-1)^{n+1}=k^{n+1}.$$
But the expression  $i^{n+1}-(i-1)^{n+1}$ can be reduced to $$i^{n+1}-\left[\binom{n+1}{0}i^{n+1}-\binom{n+1}{1}i^{n}+...+(-1)^r\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\binom{n+1}{n+1}i^{0}\right]$$
which simplifies to $$\binom{n+1}{1}i^{n}+...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}+...+(-1)^{n}\binom{n+1}{n+1}i^{0}.$$
Thus $$\sum_{i=1}^k\left[\binom{n+1}{1}i^{n}+...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}+...+(-1)^{n}\binom{n+1}{n+1}i^{0}\right]=k^{n+1}.$$
Therefore $$(n+1)\sum_{i=1}^ki^{n}+\sum_{i=1}^k\left[-\binom{n+1}{2}i^{n-1}...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n}\right]=k^{n+1}.$$
We can then see that:$$(n+1)\sum_{i=1}^ki^{n}=k^{n+1}+\sum_{i=1}^k\left[\binom{n+1}{2}i^{n-1}...+(-1)^{r}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\right].$$
Finally for all 2 $\leq r \leq n+1$, $$\sum_{i=1}^ki^{n}=\frac{1}{n+1}\left(k^{n+1}+\sum_{i=1}^k\left[\binom{n+1}{2}i^{n-1}...+(-1)^{r}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\right]\right).$$
Note that this can be solved iteratively. E.g. once you find the value for n=1, you can use this for n=2, and n=3 and so on...
A: Some of the other answers say this in a more complicated way, but I've seen so much confusion around this question that I'd like to make it as clear as possible.
There is no good formula for $\sum_1^k i^n$ for the same reason that there is no good formula for $\int_0^x t(t-1)\ldots(t-n+1)dt$, namely that we are using the wrong basis for the space of polynomials.
For integrals, $t^n$ is a simpler choice (even an eigenbasis!): $\int_0^x t^n dt = \frac{1}{n+1} t^{n+1}$.  If we want to compute the integral above, we will probably expand out and then use this formula.
With summation, or "discrete calculus", it is the other way around.  $\sum_1^k i^n$ is messy, but $\sum_1^k i^{\underline{n}} = \sum_1^k i(i-1)\ldots(i-n+1)$ is very simple, 
just $\frac{1}{n+1} (k+1)^{\underline{n+1}} = \frac{1}{n+1} (k+1)k(k-1)\ldots (k-n+1)$.  (This doesn't give us an eigenbasis, but that's because we should really be summing up to $k-1$.)
I prefer to write this instead using binomial coefficients, ${i\choose n} = \frac{1}{n!}i^{\underline{n}}$.  Then the formula is $\sum_1^k {i\choose n} = {{k+1}\choose{n+1}}$.  This is easy to prove with a combinatorial argument or a number of other methods (see my post here).
But any formula for $\sum_1^k i^n$ will involve some complexity, and most likely be equivalent to expanding in terms of the "correct" basis.  Even the nice-looking Eulerian numbers formula has two layers of complexity, one in the outer sum and one in the combinatorics of the numbers themselves.
A: At the bottom, this is very close to the other approaches, but the details are somewhat different, and it adds a nice perspective to the question: 
Consider the function $$ B(t,x)=\frac{te^{tx}}{e^t-1}. $$ Expanding as a power series o=in $t$, we see that $$ B(t,x)=\sum_{m=0}^\infty B_m(x)\frac {t^m}{m!} $$ for some $B_m(x)$. (We can treat this formally, or check that it has a positive radius of convergence, and argue there.) 
Note that $B(t,x+1)=te^{xt}+B(t,x)$, so 
$$ \sum_{m=0}^\infty(B_m(x+1)-B_m(x))\frac{t^m}{m!}=te^{xt}=\sum_{m=0}^\infty\frac{mx^{m-1}t^m}{m!}, $$
or $B_m(x+1)-B_m(x)=mx^{m-1}$.
From this it follows easily that each $B_m(x)$ is a polynomial in $x$ of degree $m$, with rational coefficients; these are the Bernoulli polynomials, and we have $$ \sum_{i=0}^n (B_{m+1}(i+1)-B_{m+1}(i))=\sum_{i=0}^n(m+1)i^m, $$ so
 $$ \sum_{i=0}^n i^m = \frac{B_{m+1}(n+1)-B_{m+1}(0)}{m+1}. $$
From a personal perspective, let me add that this version of the formula led to a nice result, discussed here: For any integer $k>2$ there are only finitely many primes of the form $p=kn+1$ such that $1^k,2^k,\dots,n^k$ are distinct modulo $p$. (This is false for $k=1,2$, and related to a nice open problem.)
A: There is a result which does not depend on the Bernouli numbers or the Eulerian numbers. I found it in my first year at university, so you can be sure it involves nothing too complicated.
The method is that of discrete calculus, first I define $ \Delta f(x) = f(x+1)-f(x) $ and $ \Sigma f(x) = \sum\limits_{k=0}^{x-1} f(k) $ as discrete equivalents of differentiation and integration respectively.
The result will follow from a discrete equivalent of integration by parts, which I first proove:
$ f(x)g(x) = \Sigma f(x+1)g(x+1) -\Sigma f(x)g(x) $
$ \hspace{30 pt} = \Sigma f(x+1)g(x+1) -\Sigma f(x)g(x) + \Sigma f(x)g(x+1) - \Sigma f(x)g(x+1) $
$ \hspace{30 pt} = \Sigma (f(x+1)-f(x))g(x+1) + \Sigma f(x)(g(x+1)-g(x)) $
$ \hspace{30 pt} = \Sigma (\Delta f(x))g(x+1) + \Sigma f(x) \Delta g(x) $
$ \therefore \space \Sigma f(x) \Delta g(x) = f(x)g(x) - \Sigma (\Delta f(x))g(x+1) $
The expression we want a formula for will be a polynomial one degree higher than the power n, so we choose f(x) and g(x) with consideration that $ \deg(f(x)g(x)) $ should be $ n+1 $.
For a dash of simplicity, we can get $ f(x) \Delta g(x) = x^{n} $ by choosing $ f(x) = x^{n-1} $ and $ g(x) = \frac{1}{2} x(x-1) $.
For the other terms in the equivalence we need:
$ f(x)g(x) = \frac{1}{2} (x-1) x^n $
$ (\Delta f(x))g(x+1) = ((x+1)^{n-1}-x^{n-1})(\frac{1}{2} x(x+1)) = \frac{1}{2} x(x+1) \sum\limits_{m = 0}^{n-2} \binom{n-1}{m} x^m $
$ \hspace{57 pt} = \frac{1}{2} x(x+1)((n-1) x^{n-2} + \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} x^m) $
$ \hspace{57 pt} = \frac{1}{2}(n-1) x^n + \frac{1}{2}(n-1) x^{n-1} + \frac{1}{2} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (x+1) x^{m+1} $
Finally, plugging these in and rearranging:
$ \sum\limits_{k=0}^{x-1} k^n = \frac{1}{2} (x-1) x^n - \sum\limits_{k=0}^{x-1} \left( \frac{1}{2}(n-1) k^n + \frac{1}{2}(n-1) k^{n-1} + \frac{1}{2} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} \right) $
$ 2\sum\limits_{k=0}^{x-1} k^n = (x-1) x^n - (n-1) \sum\limits_{k=0}^{x-1} k^n - (n-1) \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} $
$ (n+1) \sum\limits_{k=0}^{x-1} k^n = (x-1) x^n + (1-n) \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} $
$ \sum\limits_{k=0}^{x-1} k^n = \frac{x-1}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $
$ \therefore \sum\limits_{k=1}^{x} k^n = x^n + \sum\limits_{k=0}^{x-1} k^n = x^n + \frac{x-1}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $
$ \hspace{31 pt} = \frac{x+n}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=1}^{x-1} k^{n-1} - \sum\limits_{k=1}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $
Note that the power of k in the sums on the RHS does not exceed $ n-1 $.
Example with $ n = 3 $; the formula simplifies down to: 
$$ \sum\limits_{k=1}^{x} k^3 = \frac{1}{4} (x^4 + 3x^3 -3 \sum\limits_{k=1}^{x-1} k^2 - \sum\limits_{k=1}^{x-1} k) $$
Which further simplifies to the correct polynomial. Only the partial sums for $k^2$ and $k$ need be known, and this formula will provide that of any higher degree.
A: I will give an approach of complex analysis.
Write
$$f(k,n)=\sum_{i=1}^k i^n.$$
It is clear that $f_n(k):=f(k,n)$ is a polynomial of $k$ for fixed $n$.
Hence, we may pretend that $f_n$ is a function defined for all $x\in\mathbb{R}$.
Let $g_n(z)$ be the analyitic continuation of $f_n(k)$, we claim that
$$g_n(z+1)=g_n(z)+(z+1)^n.$$
Let
$$h_n(z)=g_n(z+1)-g_n(z)-(z+1)^n,$$
then $h_n(z)$ is a complex polynomial and has infinite zeroes since each
$z\in \{1,2,\cdots\}$ is. By the fundamental theroem of algebra, we know
$h_n$ vanishes identically for each fixed $n$, which is our claim. Hence,
for any $x\in \mathbb{R}$, we have
$$f(x+1,n)=f(x,n)+(x+1)^n.$$
Moreover, we have
$$f_n'(x+1)=f_n'(x)+n(x+1)^{n-1},$$
which implies
$$f_n'(k)=f_n'(0)+n\sum_{i=0}^{k-1}(i+1)^{n-1}=f_n'(0)+nf_{n-1}(k).$$
Similarly, we have
$$f_n'(x)=f_n'(0)+nf_{n-1}(x),\forall x\in\mathbb{R}.$$
By the fundamental theorem of integral, we have
$$f(k,n)=\int_0^{k}f_n'(x)dx+f(0,n)=n\int_0^{k}f(x,n-1)dx+kf_n'(0).$$
Take $k=1$, we have
$$f_n'(0)=1-n\int_0^1 f(x,n-1)dx. $$
We conclude that
$$f(k,n)=n\int_0^{k}f(x,n-1)dx+k\left(1-n\int_0^1 f(x,n-1)dx\right),$$
which gives us an induction formula.
For example, note that $f(k,1)=k(k+1)/2$, then
$$f(x,1)=\frac{x(x+1)}{2},$$
and then
$$f(k,2)=2\int_0^{k}\frac{x(x+1)}{2}dx+k\left(1-2\int_0^1 \frac{x(x+1)}{2}dx\right)=\frac{k(k+1)(2k+1)}{6}.$$
