Compilation of proofs for the summation of natural squares and cubes

Question

I want to know different proofs for the following formulas,

$$ \sum_{i=1}^n{i^2} = \frac{(n)(n+1)(2n+1)}{6} $$

$$ \sum_{i=1}^n{i^3} = \frac{n^2(n+1)^2}{2^2} $$

Please do not mark this as duplicate, since what I specifically want is to be exposed to a variety of proofs using different techniques (I did not find such a compilation anywhere on the net)

I am only familiar with two proofs, one which uses expansion of $(x+1)^2 - x^2$ and $(x+1)^3 - x^3$ and the other which uses induction. I have provided link for the induction proof in a self-answer.

I am particularly interested in proof without words and proofs which use a unrelated mathematical concept (higher level math upto class 12 level is acceptable).

Also,

(1) Don't think I am being rude or anything, it is out of genuine interest that I am asking this question.

(2) Someone marked this as a duplicate of Methods to compute $\sum_{k=1}^nk^p$ without Faulhaber's formula

My question is different in three ways:

(i) I want to focus only on these two summation and not the general case,

(ii) Hence, it follows that the proofs which I am looking for a simpler than the ones provided in that link and are simpler (using images, pictures or high school algebra). What I want is to study new proofs. I believe it is a good practice when learning math to so this.

(iii) Since the proofs in the link are given for the general case, they are complicated and I am finding it hard to understand them. If someone is able to use the same method to the two cases in my question, then it would probably become much simpler and easier to digest.

Appendix

Feel free to make use of these topics in your answers,

Calculus

Basic Binomial Expansion

Coordinate Geometry

Algebra (upto what 18 year olds learn)

Taylor Series Expansions

Geometry (18 year old level)

basically.............math which 18 year old's learn on Earth.

If you want to err, then err on the higher math side:)

Answers Compilation List

1. By Newton series
2. By Sterling Numbers
3. By Induction
4. From the book Generatingfunctionology by Herbert Wilf

5. By generalizing the following pattern $$\begin{align} &\ \,4\cdot5\cdot6\cdot7\\=&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6+4\cdot5\cdot6\cdot7)\\-&(0\cdot1\cdot2\cdot3+1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6)\\ =&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot4+3\cdot4\cdot5\cdot4+4\cdot5\cdot6\cdot4)\\ \end{align}$$

6. By Lagrangian Interpolation

7. By Formal Differentiation
8. By the Euler-Maclaurin Summation Formula
9. By Assuming that the expression is a polynomial of degree $2$.
10. A Proof Without Words for the cube case
11. By integrating and assuming a error term.
12. SimplyBeatifulArt's Personal Approach

Answer 1

For the sum of cubes there is a nice illustration in the wikipedia article Squared Triangular Number Illustrating Nicomachus's theorem.

$$\sum_{i=1}^ni^3=\left(\sum_{i=1}^ni\right)^2$$

(Nicomachus Theorem by Cmglee)

Answer 2

By integrating and assuming an error term

Supposing known: $$\sum_{i=1}^{n}{1}=n \tag{*}$$ and $$\sum_{i=1}^{n}{i}=\frac{n(n+1)}{2}\tag{**}$$

Let $\sigma$ be: $$\sigma(n)=\sum_{i=1}^{n}{i^2}$$ This is similar to the integral: $$\Sigma(n)=\int_{0}^{n}{x^2\,dx}$$ up to an error $\epsilon(n)$. Therefore: $$\sigma(n) = \Sigma(n) +\epsilon(n) \tag{1}$$ This error may be obtained by the fact that the curve $x^2$ is always below of the actual value $i^2$, hence applying $(*)$ and $(**)$: $$\epsilon(n) = \sum_{i=1}^{n}{\left[i^2-\int_{i-1}^{i}{x^2\,dx}\right]}=\sum_{i=1}^{n}{\left[i-\frac{1}{3}\right]}=\frac{n(n+1)}{2}-\frac{n}{3}$$ Integrating $\Sigma(n)$ and adding them up in $(1)$ the desired result is reached: $$\sigma(n) = \frac{n^3}{3}+\frac{n(n+1)}{2}-\frac{n}{3}=\frac{n(2n+1)(n+1)}{6}$$

The one for your second required proof involves the previous result, and it is done in a similar fashion. Try it!

Answer 3

My personal approach:

Let $f(x,p)$ be a polynomial of $x$ with degree $p+1$ such that for $n\in\mathbb N$, we have

$$f(n,p)=\sum_{k=1}^nk^p$$

This function satisfies the recursive relation:

$$f(x,p)=a_px+\int_0^xf(t,p-1)~\mathrm dt$$

where

$$a_p=1-p\int_0^1f(t,p-1)~\mathrm dt$$

For example, if know that

$$f(x,1)=\frac{x(x+1)}2$$

Then it follows that $a_1=\frac12$ and

$$f(x,2)=\frac12x+\int_0^x\frac{t(t+1)}2~\mathrm dt=\frac{x(x+1)(2x+1)}6$$

which is the formula for $\sum k^2$.

---

As described in this answer, we have

$$\sum_{k=1}^nk^p=\sum_{k=1}^p{p\brace k}\frac{(n+1)^{\underline{k+1}}}{k+1}$$

Where ${p\brace k}$ is a stirling number of the second kind and $a^{\underline b}=a(a-1)(a-2)\dots(a-b+1)$ is a falling factorial. For your cases,

$$\sum_{k=1}^nk^2={2\brace1}\frac{(n+1)(n)}2+{2\brace2}\frac{(n+1)(n)(n-1)}3\\=\frac{n(n+1)}2+\frac{n(n+1)(n-1)}3\\=\frac{n(n+1)(2n+1)}6$$

$$\sum_{k=1}^nk^3={3\brace1}\frac{(n+1)(n)}2+{3\brace2}\frac{(n+1)(n)(n-1)}3+{3\brace3}\frac{(n+1)(n)(n-1)(n-2)}4\\=\frac{n(n+1)}2+n(n+1)(n-1)+\frac{n(n+1)(n-1)(n-2)}4\\=\left[\frac{n(n+1)}2\right]^2$$

Of course, for normal purposes, just use the table of values to evaluate basic stirling numbers.

Answer 4

I'm surprised that no one has pointed out the elementary Newton series, which I'll describe below.

---

The forward difference operator $Δ$ defined as: $ \def\nn{\mathbb{N}} \def\zz{\mathbb{Z}} \def\lfrac#1#2{{\large\frac{#1}{#2}}} \def\lbinom#1#2{{\large\binom{#1}{#2}}} $

$Δ = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) - f(n) ) )$

Namely, for any function $f$ on $\zz$ and $n \in \zz$, we have $Δ(f)(n) = f(n+1) - f(n)$.

If you think of the functions as sequences (infinite in both directions), then taking the forward difference means replacing each term with the value of the next term minus itself. For example if you repeatedly take the forward difference of the sequence of cubes:

...,-27,-8,-1, 0, 1, 8,27,...
..., 19, 7, 1, 1, 7,19,37,...
...,-12,-6, 0, 6,12,18,24,...
...,  6, 6, 6, 6, 6, 6, 6,...
...,  0, 0, 0, 0, 0, 0, 0,...
...,  0, 0, 0, 0, 0, 0, 0,...

This powerful abstraction makes it easy to get a lot of things. For instance, the numbers obtained here can be easily used to obtain the general formula for sum of cubes, as you desire.

General method for indefinite summation

The key is that:

$Δ\left( \text{int $n$} \mapsto \lbinom{n}{k+1} \right) = \left( \text{int $n$} \mapsto \lbinom{n}{k} \right)$ for any $k \in \zz$.

This is to be expected because it follows directly from extending Pascal's triangle naturally, namely if we define $\lbinom{n}{k}$ by the recurrence:

$\lbinom{n}{0} = 1$ for any $n \in \zz$.

$\lbinom{0,k+1}{0} = 0$ for any $k \in \nn$.

$\lbinom{n+1}{k+1} = \lbinom{n}{k+1} + \lbinom{n}{k}$ for any $k \in \nn$ and $n \in \zz$.

Now consider any function $f$ on $\zz$ such that $f(n) = \sum_{k=0}^d a_k \lbinom{n}{k}$ for any $n \in \zz$. Then we have for any $m \in \nn_{\le d}$:

$Δ^m(f)(n) = \sum_{k=0}^{d-m} a_{k+m} \lbinom{n}{k}$ for any $n \in \zz$.

And hence:

$Δ^m(f)(0) = a_m$.

Which immediately gives Newton's series:

$f(n) = \sum_{k=0}^d Δ^k(f)(0) \lbinom{n}{k}$ for any $n \in \zz$.

From a high-level perspective, this is the discrete version of the Taylor series.

This works for any polynomial function $f$ on $\zz$, since $D^k(f)$ is the zero function once $k$ is larger than the degree of $f$, so we can use it to immediately find the series for $(\text{int n} \mapsto n^3)$, and then just take the anti-difference by shifting the coefficients of the series the other way. The undetermined constant that appears will drop out once we perform a definite sum like if we want the sum of the first $m$ cubes.

---

Sum of $p$ powers

For example if we want $\sum_{k=1}^{n-1} k^3$ we first find the iterated forward differences of the sequence of cubes $( n^3 )_{n \in \zz}$:

..., 0, 1, 8,27,64,...
..., 1, 7,19,37,...
..., 6,12,18,...
..., 6, 6,...
..., 0,...

So we immediately get $n^3 = 0 \binom{n}{0} + 1 \binom{n}{1} + 6 \binom{n}{2} + 6 \binom{n}{3}$ and hence $\sum_{k=0}^{n-1} k^3 = 0 \binom{n}{1} + 1 \binom{n}{2} + 6 \binom{n}{3} + 6 \binom{n}{4} = \lfrac{n(n-1)}{2} \Big( 1 + \lfrac{6(n-2)}{3} \big( 1 + \lfrac{n-3}{4} \big) \Big) = \Big( \lfrac{n(n-1)}{2} \Big)^2$.

Computation efficiency

This is far more efficient than the method you mentioned in your question (namely by taking summation on both sides of $(n+1)^3-n^3 = 3n^2+3n+1$ and telescoping) because the Newton series is easy to manipulate and easy to compute. For sum of $p$-powers we only need $O(p^2)$ arithmetic operations to find the forward-differences and then $O(p^2)$ more to simplify the series form into a standard polynomial form. In contrast, the other method requires $O(p^3)$ arithmetic operations.

Answer 5

By formal differentiation:

We express the sum as a polynomial in $x$ and evaluate it at $x=1$ (taking the limit, by L'Hospital).

$$S_0(n)=\sum_{k=1}^nx^k=\frac{x^{n+1}-x}{x-1}\to n.$$

$$S_1(n)=\sum_{k=1}^nkx^k=xS_0'(n)=\frac{nx^{n+2}-(n+1)x^{n+1}+x}{(x-1)^2}\to\frac{n^2+n}2.$$

$$S_2(n)=\sum_{k=1}^nk^2x^k=xS_1'(n)=\frac{n^2x^{n+2}-(2n^2+2n-1)x^{n+1}+(n^2+2n)x^n-x-1}{(x-1)^3}\to\frac{2n^3+3n^2+n}6.$$

$$S_3(n)=\sum_{k=1}^nk^3x^k=xS_2'(n)=\frac{(n^3-n^2)x^{n+3}-(3n^3-5n+2)x^{n+1}+(3n^3+3n^2-5n-1)x^n-(n^3+2n^2)x^{n-1}+2n+4}{(x-1)^4}\to\frac{n^4+2n^3+n^2}4.$$

---

The computation can be somewhat eased by the change of variable $x=t+1$, such that you keep working with polynomials (and the Binomial theorem).

$$S_0(n)=\sum_{k=0}^{n-1}x^k=\frac{x^n-1}{x-1}=\frac{(t+1)^n-1}t=\sum_{i=1}^n\binom nit^{i-1}.$$

$$S_1(n)=xS_0'(n)=\sum_{k=0}^{n-1}kx^k=(t+1)\sum_{i=2}^n\binom ni(i-1)t^{i-2}=\sum_{i=2}^n\binom ni(i-1)(t^{i-1}+t^{i-2})\to\binom n1.$$

$$S_2(n)=xS_1'(n)=\sum_{k=0}^{n-1}kx^k=(t+1)\sum_{i=2}^n\binom ni(i-1)((i-1)t^{i-2}+(i-2)t^{i-3})\\ =\sum_{i=2}^n\binom ni(i-1)((i-1)t^{i-1}+(2i-3)t^{i-2}+(i-2)t^{i-3})\to\binom n2+2\binom n3.$$

Answer 6

There is a neat proof for the sum of squares using generating functions in the book Generatingfunctionology by Herbert Wilf, pages 35 - 36. First take the generating function $\sum_{n=0}^N x^n = \frac{x^{N+1}-1}{x-1}$. Then, we convert it into the generating function $\sum_{n=0}^N n^2x^n$ by applying $(xD)^2$ to both sides (this is a standard technique that you can learn about in the book): \begin{align*} (xD)^2\sum_{n=0}^N x^n &= (xD)\sum_{n=0}^N x\cdot nx^{n-1} = (xD)\sum_{n=0}^N nx^n = \sum_{n=0}^N x\cdot n^2x^{n-1} = \sum_{n=0}^N n^2x^n = (xD)^2\frac{x^{N+1}-1}{x-1} \end{align*} Now evaluate both sides at $x=1$ (take the limit as $x\to1$) and massage the resulting expression until you get what you want! Bashing out the derivative and limit on the right hand side takes forever though, I didn't have the guts to do it out.

Of course, you could do the same for $n^3$ but that would be even more bashing than for $n^2$.

Answer 7

$$\begin{align} &\ \,4\cdot5\cdot6\cdot7\\=&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6+4\cdot5\cdot6\cdot7)\\-&(0\cdot1\cdot2\cdot3+1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot5+3\cdot4\cdot5\cdot6)\\ =&(1\cdot2\cdot3\cdot4+2\cdot3\cdot4\cdot4+3\cdot4\cdot5\cdot4+4\cdot5\cdot6\cdot4)\\ \end{align}$$

By generalizing the pattern, this is enough to establish

$$n^{(m)}=\sum_{k=1}^n k^{(m)}-\sum_{k=0}^{n-1} k^{(m)}=(m+1)\sum_{k=1}^nk^{(m-1)},$$ where $a^{(b)}$ denotes the raising factorial.

From this,

$$n=\sum_{k=1}^n 1=S_0(n)$$

$$n(n+1)=2\sum_{k=1}^n k=S_1(n)$$ $$n(n+1)(n+2)=3\sum_{k=1}^n k(k+1)=S_2(n)$$ $$n(n+1)(n+2)(n+3)=4\sum_{k=1}^n k(k+1)(k+2)=S_3(n).$$

They are the "linear", triangular, tetrahedral, "hypertetrahedral" numbers.

Then

$$\sum_{k=1}^n 1=S_0(n)=n$$

$$\sum_{k=1}^n k=\frac12S_1(n)=\frac{n^2+n}2$$

$$\sum_{k=1}^n k^2=\frac13S_2(n)-\sum_{k=1}^n k=\frac{2n^3+3n^2+n}6$$

$$\sum_{k=1}^n k^3=\frac14S_3(n)-3\sum_{k=1}^n k^2-2\sum_{k=1}^n k=\frac{n^4+2n^3+n^2}4$$

One can continue in a systematic way.

Answer 8

While it might be considered overkill, use of the Euler-Maclaurin Summation Formula facilitates evaluation of the sum $\sum_{i=1}^n i^m$ for any integer $m\ge 0$.

Invoking the EMSF for $f(i)=i^2$ yields

$$\begin{align} \sum_{i=1}^n i^2&=\int_0^n x^2\,dx +\frac12 (n^2)+\frac1{12}\left(2n\right)\\\\ &=\frac{n^3}{3}+\frac12 n^2 +\frac16 n\\\\ &=\frac n6 (2n+1)(n+1) \end{align}$$

and invoking the EMSF for $f(i)=i^3$ yields

$$\begin{align} \sum_{i=1}^n i^3&=\int_0^n x^3\,dx +\frac12 (n^3)+\frac1{12}\left(3n^2\right)\\\\ &=\frac{n^4}{4}+\frac{n^3}{2}+\frac{n^2}{4}\\\\ &=\frac{n^2}4 (n+1)^2 \end{align}$$

as expected!

Answer 9

By Lagrangian interpolation.

The expression of the sum of the cubes must be a quartic polynomial, because its first order difference is cubic. Furthermore, it has no constant term because the sum of no numbers is zero.

Hence, the average of the cubes is a cubic polynomial, by the four points

$$(1,1),\left(2,\frac{1+2^3}2\right),\left(3,\frac{1+2^3+3^3}3\right),\left(4,\frac{1+2^3+3^3+4^4}4\right).$$

The resquested polynomial is the Lagrangian interpolant.

---

Average of ones (constant, by $(1,1)$):

$$\overline S_0(n)=1$$

Average of naturals (linear, by $(1,1)$ and $(2,\frac32)$):

$$\overline S_1(n)=\frac{\frac32-1}{2-1}(n-1)+1=\frac{n+1}2$$

Average of squares (quadratic, by $(1,1)$, $(2,\frac52)$ and $(3,\frac{14}3)$):

$$\overline S_2(n)=\cdots=\frac{2n^2+3n+1}2$$

Average of cubes (cubic, by $(1,1)$, $(2,\frac92)$, $(3,\frac{26}3)$ and $(4,\frac{45}2)$):

$$\overline S_3(n)=\cdots=\frac{n^3+2n^2+n}4$$

---

An explicit formula can be obtained for any degree, by means of the Lagrangian interpolation formula for equidistant points.

Newtonian interpolation could be easier.

Answer 10

One technique is the Newton series, which permits termwise summation just as the Taylor series permits termwise integration. There is also the more powerful but less efficient method of indefinite summation, which is unsurprisingly the discrete version of integration by parts. $ \def\zz{\mathbb{Z}} \def\lbinom#1#2{{\large\binom{#1}{#2}}} $

---

Let $R$ be the right-shift operator defined as:

$R = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) ) )$

Namely, for any function $f$ on $\zz$ and $n \in Z$, we have $R(f)(n) = f(n+1)$.

For convenience we also define the summation operator:

$Σ = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto \sum_{k=0}^{n-1} f(k) ) )$

Then we have the important property that $ΔΣ(f) = f$ for any function $f$ on $\zz$, analogous to the fundamental theorem of calculus.

Now we can derive summation by parts from the product rule in the same manner as for the continuous version. Consider any functions $f,g$ on $\zz$:

$Δ(f·g)(n) = f(n+1) g(n+1) - f(n) g(n) = f(n+1) Δ(g)(n) - Δ(f)(n) g(n)$ for any $n \in \zz$.

And hence the discrete product rule (with the usual pointwise sum and product of functions):

$Δ(f·g) = R(f)·Δ(g) + Δ(f)·g$.

Now by substituting $f$ with $Σ(f)$ and taking summation on both sides we get summation by parts:

$Δ(Σ(f)·g) = R(Σ(f))·Δ(g) + f·g$.

$Σ(f·g) = Σ(f)·g - Σ(R(Σ(f))·Δ(g)) + c$ for some constant function $c$ on $\zz$.

---

For convenience we let "$Δ_n E(n)$" denote "$E(n+1)-E(n)$" and "$Σ_n E(n)$" denote "$\sum_{k=0}^{n-1} E(k)$".

We can therefore compute $\sum_{k=0}^{n-1} k^3$ as follows:

$\sum_{k=0}^{n-1} k^3 = Σ_n\Big(\lbinom{n}{0}·n^3\Big)$   [Recall this definition of binomial coefficients.]

$\ = \lbinom{n}{1}·n^3 - Σ_n\Big(\lbinom{n+1}{1}·(3n^2+3n+1)\Big) + c$, where $c$ is some constant (independent of $n$),

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + Σ_n\Big(\lbinom{n+2}{2}·(6n+6)\Big) + c'$,

    where $c'$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - Σ_n\Big(\lbinom{n+3}{3}·6\Big) + c''$,

    where $c''$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - \lbinom{n+3}{4}·6 + c''$.

By substituting $n = 0$ we readily get $c'' = 0$ and hence simplifying gives the desired formula.

---

Using this we can also easily compute things like $\sum_{k=0}^{n-1} k^3 3^k$ by applying it three times, each time reducing the degree of the polynomial part. Incidentally, there is another purely discrete technique to compute this sum, namely subtracting it from $3 \sum_{k=1}^n k^3 3^k = n^3 3^{n+1} + \sum_{k=1}^n (k-1)^3 3^k$, to reduce the polynomial degree, and repeating three times. But interestingly this technique does not work for $\sum_{k=1}^n k^3$.

Answer 11

By indeterminate coefficients.

The sum of the cubes must be a quartic polynomial with no independent term, so

$$S_3(n)=an^4+bn^3+cn^2+dn.$$

With $n=1,2,3,4$,

$$\begin{cases}a+b+c+d&=1,\\16a+8b+4c+2d&=9,\\ 81a+27b+9c+3d&=36,\\ 256a+64b+16c+4d&=100.\\ \end{cases}$$

$$\begin{cases}14a+6b+2c&=7,\\ 78a+24b+6c&=33,\\ 252a+60b+12c&=96.\\ \end{cases}$$

$$\begin{cases} 36a+6b&=12,\\ 168a+24b&=54.\\ \end{cases}$$

$$\begin{cases} 24a&=6.\\ \end{cases}$$

Then $b=\dfrac12,c=\dfrac14,d=0$.

---

There is a faster way, using

$$S_3(n)-S_3(n-1)=n^3\\ =a(n^4-(n-1)^4)+b(n^3-(n-1)^3)+c(n^2-(n-1)2)+d(n-(n-1))\\ =4an^3+(-6a+3b)n^2+(4a-3b-2c)n+(-a+b-c+d).$$

Then by identifying the coefficients, you get an easy triangular system:

$$\begin{cases}-a+b-c+d&=0,\\4a-3b+2c&=0,\\-6a+3b&=0,\\4a&=1.\end{cases}$$

Notice that the system is easy to write, as the columns are the rows of Pascal's triangle with alternating signs and the last coefficient missing.

Answer 12

This is not the most elegant way to do it, but there is a direct visualization of a related sum which directly implies that first formula.

Starting one step lower...

Consider the following visual proof for the sum of counting numbers: they arrange into a triangle which can be overlaid upon itself to form a square, duplicating the diagonal: The implication is that the sum of counting numbers, being visually represented as one of these triangles, must be $(n^2 + n)/2$ which can afterwards be simplified to $n(n+1)/2,$ and the key focus that I want to highlight is this double-counting of stuff on the diagonal.

This is the basic mechanism that we're going to use to directly visualize the sum of the first $m$ odd squares.

Going from odds to general sums.

There is also a slightly more direct visualization also that the sum of the first $n$ odd numbers is the $n^\text{th}$ square number: to see this just color the top left corner of an $n\times n$ square green, then color the three cells beside it blue, then color the 5 cells beside those cells red, and so on... one immediately sees that each color highlights a subsequent odd number. They are odd, because they are made up of two equal sides plus the one square on the diagonal.

One can then use this sum of odd numbers to prove the same formula. Let me only consider the case where we end on an even number: sum up the counting numbers from $1$ to $2m$; this consists of the first $m$ odd numbers plus the first $m$ even numbers (if the "first" even number is regarded as 2, not 0). Subtracting 1 from each even number gives the corresponding odd number, so the sum is $2m^2 + m = n^2/2 + n/2.$ One can make a similar argument for odds or just add $n+1$ to this form, which is very similar to a proof by induction but lacks the inductive step.

Extending this directly

The obvious 3D generalization of that first overlay diagram concerns the sum of odd squares, as we build a staggered pyramid with a $1\times1$ square for its topmost layer supported by a $3\times3$ square for its second layer supported by a $5\times5$ square for its third layer, and so on, up to a $(2n-1)\times(2n-1)$ square at the bottom. I don't have a great 3D visualization package for this sort of thing so the quickest I could mock up looks like this:

We can use 6 of these staggered pyramids to complete a supercube with sides of length $2n-1$, if we can count the overlaid squares. We start by counting every unit subcube of the overall supercube once, $(2n-1)^3$, and this gets most of what we need. Looking at just one of the pyramids above we see that there is the "topmost" square (in blue) which needs to be counted 6 times total, beside it each of the edges of length $n-1$ (in gold) needs to be counted 3 times total, and flanked by those edges there is an "inner triangle" which needs to be counted 2 times total. Since we've already counted each of these once we need 5 more for the teal cube, 2 more of each gold cube, and one more of these grassy cubes forming the inner triangle on the outside. Of these, the inner triangles are the only ones that look daunting -- until you realize that they are a sum of the first $n-1$ odd numbers and therefore are just $(n-1)^2,$ so subtracting the one time we've already counted each of these squares we need to add $5 + 2\cdot(n-1)\cdot(\text{# of edges}) + 1\cdot(n-1)^2(\text{# of inner triangles})$ yielding the expression: $$\begin{align}6\sum_{k=1}^{n} (2k-1)^2 =& (2n-1)^3 + 12 (n - 1)^2 + 16 (n-1) + 5\\ =&~(2n-1)~2n~(2n+1).\end{align}$$Now if we have a sum of the squares from $1$ to $4m^2$ we find that this is a sum of the first $m$ odd squares plus the first $m$ even squares, but each even square is a subsequent odd number above an odd square: $3 + 7 + 11 + \dots + (4m-1).$ Subtracting $1$ from all $m$ of these and dividing by $2$ gives this as $2(1 + 3 + 5 + \dots + (2m-1)) + m = 2m^2 + m,$ so our final result is $$ \sum_{k=1}^{2m} k^2 = \frac{(2m-1)~2m~(2m+1)}{3} + 2m^2 + m = \frac13 m (4m+1)(2m + 1),$$ from which the result follows.

Answer 13

I've came across this post and wanted to add another proof using generating functions. First of all we will denote $$g(n)=\sum_{k=1}^nk^3,\;\;f(n)=\sum_{k=1}^nk$$ and their generating functions by

$$G(x)=\sum_{n=0}^{\infty}g(n)x^{n},\;\;F(x)=\sum_{n=0}^{\infty}f(n)x^{n},S(x)=\sum_{n=0}^{\infty}f^{2}(n)x^{n}$$ Notice that $$ \sum_{n=1}^{\infty}\left(f(n)-f(n-1)\right)x^{n}=\sum_{n=1}^{\infty}nx^{n}=x\sum_{n=1}^{\infty}nx^{n-1}=x\frac{d}{dx}\sum_{n=0}^{\infty}x^{n}=x\frac{d}{dx}\frac{1}{1-x}=\frac{x}{(1-x)^{2}} $$ On the other hand $$ \sum_{n=1}^{\infty}\left(f(n)-f(n-1)\right)x^{n}=\sum_{n=1}^{\infty}f(n)x^{n}-\sum_{n=1}^{\infty}f(n-1)x^{n}\stackrel{f(0)=0}{=}F(x)-x\sum_{n=0}^{\infty}f(n)x^{n}=(1-x)F(x) $$ So we will conclude that $F(x)=\frac{x}{(1-x)^{3}}$. Now, observe that

$$ \sum_{n=1}^{\infty}\left(g(n)-g(n-1)\right)x^{n}=\sum_{n=1}^{\infty}g(n)x^{n}-\sum_{n=1}^{\infty}g(n-1)x^{n}\stackrel{g(0)=0}{=}G(x)-x\sum_{n=0}^{\infty}g(n)x^{n}=(1-x)G(x) $$ By using the identity $n^{3}=n(n-1)(n-2)+3n(n-1)+n$ we will get \begin{align*} \sum_{n=1}^{\infty}\left(g(n)-g(n-1)\right)x^{n} & =\sum_{n=1}^{\infty}n^{3}x^{n}=\\ \\ =\sum_{n=1}^{\infty}n(n-1)(n-2)x^{n} & +\sum_{n=1}^{\infty}3n(n-1)x^{n}+\sum_{n=1}^{\infty}nx^{n}=\\ \\ =x^{3}\sum_{n=1}^{\infty}n(n-1)(n-2)x^{n-3} & +3x^{2}\sum_{n=1}^{\infty}n(n-1)x^{n-2}+x\sum_{n=1}^{\infty}nx^{n-1}=\\ \\ =x^{3}\frac{d^{3}}{dx^{3}}\frac{1}{1-x}+3x^{2}\frac{d^{2}}{dx^{2}} & \frac{1}{1-x}+x\frac{d}{dx}\frac{1}{1-x}=\\ \\ =\frac{6x^{3}}{(1-x)^{4}}+\frac{6x^{2}}{(1-x)^{3}}+ & \frac{x}{(1-x)^{2}}=\frac{x^{3}+4x^{2}+x}{(1-x)^{4}} \end{align*} So $G(x)=\frac{x^{3}+4x^{2}+x}{(1-x)^{5}}$. Notice that $$ \sum_{n=1}^{\infty}\left(f^{2}(n)-f^{2}(n-1)\right)x^{n}=\sum_{n=1}^{\infty}f^{2}(n)x^{n}-\sum_{n=1}^{\infty}f^{2}(n-1)x^{n}\stackrel{f^{2}(0)=0}{=}S(x)-x\sum_{n=0}^{\infty}f^{2}(n)x^{n}=(1-x)S(x) $$ On the other hand \begin{align*} \sum_{n=1}^{\infty}\left(f^{2}(n)-f^{2}(n-1)\right)x^{n} & =\sum_{n=1}^{\infty}\underset{n}{\underbrace{(f(n)-f(n-1))}}\underset{2f(n-1)+n}{\underbrace{(f(n)+f(n-1))}}x^{n}=2\sum_{n=1}^{\infty}nf(n-1)x^{n}+\sum_{n=1}^{\infty}n^{2}x^{n}=\\ \\ =2x^{2}\sum_{n=0}^{\infty}nf(n)x^{n-1}+2x & \sum_{n=0}^{\infty}f(n)x^{n}+x^{2}\sum_{n=1}^{\infty}n(n-1)x^{n-2}+x\sum_{n=1}^{\infty}nx^{n-1}=\\ \\ =2x^{2}\frac{d}{dx}F(x)+2xF(x) & +x^{2}\frac{d^{2}}{dx^{2}}\frac{1}{1-x}+x\frac{d}{dx}\frac{1}{1-x}=\\ \\ =\frac{2x^{2}(1-x)+6x^{3}}{(1-x)^{4}}+ & \frac{2x^{2}}{(1-x)^{3}}+\frac{2x^{2}}{(1-x)^{3}}+\frac{x}{(1-x)^{2}}=\frac{x^{3}+4x^{2}+x}{(1-x)^{4}} \end{align*} Thus $S(x)=\frac{x^{3}+4x^{2}+x}{(1-x)^{5}}=G(x)$. Because the generating functions of the sequences $f^2(n),g(n)$ are equal, the sequences $g(n),f^{2}(n)$ themselves are equal.