There are some convenient formulas for the summation of the first $n$ integers which are the $a$th powers of other integers, e.g.
$$ \sum_{i=0}^n i = \frac {n(n+1)}2$$ $$ \sum_{i=0}^n i^2 = \frac {n(n + 1)(2n + 1)}6 $$
Similarly, there are formulae for integration of monomials ($ x^n $), such as
$$ \int x^2\,dx = \frac {x^3}3 , \int x\,dx = \frac {x^2}2 $$
On replacing $n + 1$ by $n$, $2n + 1$ by $2n$, and $a_n + b $ by $a_n$ for n'th term in the summation formulas, we get the integration formulas, e.g
Replace $n + 1$ by $n$, we get $ n^2/2 $ for $n(n+1)/2$ Replace $n+1$, $2n+1$, by $n, 2n$, then for $n(2n+1)(n+1)$, we get $$ \frac {n(n)(2n)}6 = \frac {2n^3}6 = \frac {n^3}3 $$
This also works for $x^3$, and $x^4$. Why does this happen? Does it work for higher powers? (I've checked that it works up to $n=8$ so far.)
NOTE:I'm learning some discrete mathematics, and I already know a little bit of calculus.
Integration is a continuous version of the sum. So your sum is actually a discrete approximation of the integral. It will have the same leading term, but the next terms differ (the sum approximates a smooth function with rectangles, as in the Riemann sum, so there are some discrepancies).
The discrete version for arbitrary power is given by the Faulhaber's formula [ http://en.wikipedia.org/wiki/Faulhaber%27s_formula ]. With a bit of thought you can see this can also help you derive the connection between a sum of a series ($\sum_n a_n$) and an integral ($\int a(x)dx$). With careful manipulation of the Taylor series of $a(x)$ under the integral, using the Faulhaber's formula in the process, you will arive at the Euler-Maclaurin formula [ http://en.wikipedia.org/wiki/Euler_maclaurin_formula ] which is very useful for calculating the series (frequently it's easier to integrate). You will notice how the Bernoulli numbers appear in both the Faulhaber's and Euler-Maclaurin's series.
It's quite fun to try and do this on paper -- you learn a lot.
EDIT:
A smooth function can be approximated as a polynomial around some point... the coefficients are derivatives of the function: $$f(x+h)=f(x)+f'(x)h+\frac{1}{2!}f''(x)h^2+\cdots$$ where $h$ is a displacement from the chosen origin $x$. Now you can approximate the values $f(1)$, $f(2)$, $f(n)$ and so on by setting $x=0$, $h=n$: $$f(n)=f(0)+f'(0)n+\frac{1}{2!}f''(0)n^2+\cdots$$ Now, take the sum you want to evaluate: $$\sum_{n=0}^N f(n)=f(0)\sum_n^N 1 +f'(0)\sum_n^N n+\frac12 f''(0)\sum_n^N n^2+\cdots$$ $$=f(0) N +f'(0)\frac{N(N+1)}{2}+\frac{1}{2!}f''(0)\frac{N(N+1)(2N+1)}{6}+\cdots$$ So, if you have the expression for $n$-th term of some series that you want to sum ($\sum a_n=\sum f(n)$) then you can express it as a sum of derivatives times the Faulhaber's polynomials that you listed in your original post.
Now, if you instead of summing just try to integrate the function (again use the Taylor series): $$\int_0^N f(x)dx=f(0)N+f'(0)\frac{N^2}{2}+\frac{1}{2!}f''(0)\frac{N^3}{3}+\cdots$$
You can now see the result is very similar... it only differs in a few higher-order terms. If you evaluate $\sum_{n=0}^N f(n)-\int_0^N f(x)dx$ you get the residual in the Euler-Maclaurin formula, which is just the extra terms that you noticed in comparison of integral to the Faulhaber's polinomial. Congratulations, you almost derived a very famous and mysterious formula from classical calculus.
In this form, it only matches the lower limit of the conventional Euler-Maclaurin, if you want the full version with derivatives on the upper limit, you have to write $\sum_{n=0}^N=\sum_{n=0}^\infty-\sum_{n={N+1}}^\infty$ and do a shifted version of above formula for the second sum (starting at $N+1$ instead of $0$.
Of course... the convergence of this is questionable (it's an asymptotic series) -- the Taylor series usually doesn't converge on all reals.
EDIT2:
And how to get the Faulhaber's polynomials? You get them by setting the polynomial with undetermined coefficients and comparing terms:
$$\sum_{n=0}^N n^k = a_0 +a_1 N + a_2 N^2+\cdots+ a_{k+1}N^{k+1}$$ $$\sum_{n=0}^N n^k=\sum_{n=0}^{N-1}n^k+N^k$$ Use the first one in the second: $$a_0 +a_1 N + a_2 N^2+\cdots+ a_{k+1}N^{k+1}=a_0 +a_1 (N-1) + a_2 (N-1)^2+\cdots+ a_{k+1}(N-1)^{k+1}+N^k$$ From $N=0$ you get $a_0=0$. Then collect terms with $N$, terms with $N^2$ and so on. You get a set of linear equations for coefficients $a_n$. What you get are the Faulhaber's polynomials which (check the wiki) can alternatively also be expressed with Bernoulli numbers and binomial coefficents.
For $k=3$, for instance, you get $$a_1 N + a_2 N^2+a_3 N^3 + a_4 N^4=a_1 (N-1) + a_2 (N-1)^2+a_3 (N-1)^3+a_4 (N-1)^4+N^3$$ which after expansion (note that the left side cancels out with the leading terms), you get: $$0=-a_1+a_2(1-2N)+a_3(-1+3N-3N^2)+a_4(1-4N+6N^2-4N^3)+N^3$$ split by powers: $$a_1-a_2+a_3-a_4=0$$ $$2a_2-3a_3+4a_4=0$$ $$3a_3-6a_4=0$$ $$4a_4=1$$ The system is always upper triangular (actually, if you write down the Matrix, it's a part of the Pascal triangle with alternating signs, there's no need to actually expand the polynomials, just write the matrix and fix the right column to (0,0,0,0,...1)), so the solution is a simple back-substitution: $$a_4=1/4$$ $$a_3=1/2$$ $$a_2=1/4$$ $$a_1=0$$ Leading to $$\sum_{n=0}^N n^3=\frac{N^2+2N^3+N^4}{4}$$
So everything can actually be done by hand. Have fun with $\sum n^4$ :)