How to sum up

\(1^2+2^2+\ldots+n^2 ?\)

Seems this question is going round in the net. These sums are attributed to various authors, the last of them is Bernoulli. I show you two tricks to do them.

Assume, we already know the lower powers, here

\(1+2+\ldots+n = \dfrac{n(n+1)}{2}.\)

Then, surprisingly, the following is useful

\(\sum\limits_{i=1}^n (i+1)^3 = \sum\limits_{i=1}^n (i^3+3i^2+3i+1).\)

If we expand the right hand side and cancel the same sums on the left and right hand side, we get

\((n+1)^3 = 1 + 3 \sum\limits_{i=1}^n i^2 + \dfrac{3n(n+1)}{2} + n.\)

From that, we get

\(\sum\limits_{i=1}^n i^2 = \dfrac{n(n+1)(2n+1)}{6}.\)

I found the factorization with Maxima.

>factor(((n+1)^3-1-3*n*(n+1)/2-n)/3)

                          n (n + 1) (2 n + 1)
                          -------------------
                                   6

Another trick is to observe that

\(p_k(n) = \sum\limits_{i=1}^n i^k\)

must be a polynomial of degree k+1. Why? Because, more generally, for any polynomial q of degree k there is polynomial p of degree k+1 such that for all x

\(p(x+1) – p(x)  = q(x+1).\)

We make sure that this is true for n from 1 to k by setting p(0)=q(0) and

\(p(1)=q(0)+q(1), \, \ldots, \, p(k+1)=q(0)+\ldots+q(k+1).\)

This is an ordinary interpolation problem in k+2 points and has a unique solution p of degree k+1. Now we note that

\(p(x+1)-p(x)\)

is a polynomial of degree k only, since the highest coefficient vanishes. And this polynomial agrees with the polynomial q(x+1) of degree k in the k+1 points x=0 to x=k. So they must be equal.

We apply this to q(x)=x^k. By Lagrange interpolation, we get

\(p_k(n) = \sum\limits_{m=0}^{k+1} s_m L_m(n)\)

with

\(s_m = 1+2^k+\ldots+m^k\)

and the Lagrange polynomials in 0,…,k+1

\(L_m(x) = \dfrac{\omega(x)}{\omega'(m)(x-m)}\)

with

\(\omega(x) = x(x-1)\cdots(x-(k+1)).\)

I don’t know if that counts as a formula for the sum. But it can be computed with Maxima in EMT easily.

>k &:= 2
 2
>function omega(x) &= prod(x-m,m,0,k+1)
 
                       (x - 3) (x - 2) (x - 1) x
 
>function domega(x) &= factor(diff(omega(x),x))
 
                                     2
                       2 (2 x - 3) (x  - 3 x + 1)
 
>function p(n) &= factor(sum(sum(i^k,i,1,m) ...
     *omega(n)/(domega(m)*(n-m)),m,0,k+1))
 
                          n (n + 1) (2 n + 1)
                          -------------------
                                   6

Here is the result for k=3.

>k &:= 3
 3
>function omega(x) &= prod(x-m,m,0,k+1)
 
                   (x - 4) (x - 3) (x - 2) (x - 1) x
 
>function domega(x) &= factor(diff(omega(x),x))
 
                      4       3        2
                   5 x  - 40 x  + 105 x  - 100 x + 24
 
>function p(n) &= factor(sum(sum(i^k,i,1,m) ...
     *omega(n)/(domega(m)*(n-m)),m,0,k+1))
 
                               2        2
                              n  (n + 1)
                              -----------
                                   4

Check

>p(1:10)
 [1,  9,  36,  100,  225,  441,  784,  1296,  2025,  3025]
>cumsum((1:10)^3)
 [1,  9,  36,  100,  225,  441,  784,  1296,  2025,  3025]