Next post Previous post

July 19, 2026 at 9:07 PM by Dr. Drang

Yesterday I watched the most recent Numberphile video, in which Ben Sparks explains a few finite series and the equations that simplify the calculation of their sums. He derives the formulas graphically, and it’s all very cleverly done, especially the one for the sum of cubes.

The idea is to get a simple polynomial expression in n for

As I said, Ben’s graphical method is really clever and easy to understand, and it leads to this expression:

I wanted to see if I could get that same result using a nongraphical and less clever method. The method that came to mind was to generate a handful of values, put them in a table, and start taking differences.

Here are some values of n, N, and the first four differences:

The differences are calculated by looking in the preceding column and subtracting the value in the same row from the value in the following row. This sort of thing is fairly easy to do by hand but is really easy to do in a spreadsheet. Here’s a screenshot of the table in Numbers, where I’m displaying the formula for the first difference, Δ:

That formula can be filled into all the difference cells, which is why it’s so easy. (I included only the first six rows in the table above because I figured you’d trust me that all the values in the fourth difference column, Δ⁴, are the same.)

The constant values in the fourth difference column mean N is a fourth-degree polynomial in n:

Because the constant fourth difference is 6, we know that

This comes from the fact that differences are analogous to derivatives, and

Once we know the degree of the polynomial and have $a_4$, we can figure out the other coefficients by solving a set of four simultaneous equations for four different values of n and N. For example:

We could use any four of the six Ns we’ve calculated, but the first four are convenient. Let’s solve them in Python using NumPy and the solve function from the linalg submodule. We can do this interactively:

python:
from numpy import np

m = np.array([[1, 1, 1, 1], [1, 2, 4, 8], [1, 3, 9, 27], [1, 4, 16, 64]])
b = np.array([1 - 1**4/4, 9 - 2**4/4, 36 - 3**4/4, 100 - 4**4/4])
np.linalg.solve(m, b)

The last line returns

python:
array([ 1.77635684e-15, -3.99680289e-15,  2.50000000e-01,  5.00000000e-01])

Those first two elements of the solution array are at the lower limit of what a floating point number can be. So we can say

Therefore,

or, after collecting terms,

which is the formula Ben got in the video.

Doing basically the same thing in Mathematica is

m = {{1, 1, 1, 1}, {1, 2, 4, 8}, {1, 3, 9, 27}, {1, 4, 16, 64}};
b = {1 - 1^4/4, 9 - 2^4/4, 36 - 3^4/4, 100 - 4^4/4};
LinearSolve[m, b]

which returns

{0, 0, 1/4, 1/2}

This is the same as the Python answer but with no need to think about floating point precision.

A third way to solve the simultaneous equations is to do it in a spreadsheet. Here’s how that looks in Numbers:

where

• the block of yellow cells is the matrix m in the Python and Mathematica solutions;
• the block of magenta cells is the inverse of m;
• the block of cyan cells is the vector b in the Python and Mathematica solutions; and
• the block of gray cells is the solution for $a_0$ through $a_3$, determined by multiplying the magenta matrix by the cyan vector.

The spreadsheet uses a combination of MINVERSE and MMULT, functions that are also in Excel and Google Sheets. As with the Python solution, there are floating point artifacts here. They’re mostly hidden by my choice to show only four decimal places, but that -0.0000 for $a_1$ is a clue that these are not exact answers.

You might well ask why I’d bother using Python or Mathematica to solve the simultaneous equations if I already had a spreadsheet open to make the difference table. The answer is that I generally prefer working in an environment where my formulas and expressions are always visible. It makes things easier to debug, and I’m always debugging.

There are other ways to determine the coefficients. My favorite is a step-by-step procedure that simplifies the problem one polynomial degree at a time.

Given that we know $a_4=1/4$ from the difference table above, we can make a new table for

and its differences. By subtracting the fourth-degree term from N, we expect P to be a third-degree polynomial. Here are the first five rows of the difference table for P:

As expected, they become constant at the third difference. Using the same technique we used to get $a_4$, we can say

Moving on, we make a table for

and its differences. Here are the first four rows of that table:

These become constant at the second difference, and we can say

Finally, we make a table for

which is this:

Since all the R values are zero, the rest of the coefficients, $a_1$ and $a_0$, must be zero, and we’re done.

If it’s not obvious why $a_1=a_0=0$, consider this:

The only way for R to be zero for all values of n is if $a_1=a_0=0$.

Building these tables is fairly easy, but it’s not as fast as building and solving the simultaneous equations. Still, there’s some satisfaction in marching to the answer this way. It’s basically a recursive solution, where we’re simplifying the problem with each step.

Using difference tables to work out the formula for the sum of cubes isn’t as slick as the graphical method shown in the video, but it does have the advantage of not requiring any ingenuity. I’m not opposed to ingenuity, but sometimes you want to just set the problem up, turn the mathematical crank, and get the answer.

Next post Previous post