Method of differences and Newton polynomials

Tags Math

I was reading about Babbage's Difference engine the other day, and stumbled upon a very interesting application of the forward differences method. It turns out that if we get a sequence generated by a polynomial, under certain conditions we can find the generating polynomial from just a few elements in the sequence.

For example, 0, 1, 5, 12, 22, 35, 51... is a sequence known as the pentagonal numbers, and we can use this technique to figure out that the polynomial \frac{3n^2}{2}-\frac{n}{2} generates it [1].

Notation

Let's start with some mathematical notation. We'll call the underlying function generating the sequence f(n). In our example, f(0)=0, f(1)=1, f(2)=5, f(3)=12 and so on.

  • The first difference is the sequence \Delta f(0)=f(1)-f(0), \Delta f(1)=f(2)-f(1), etc.
  • The second difference is the sequence \Delta^2 f(0)=\Delta f(1) - \Delta f(0), \Delta^2 f(1)=\Delta f(2) - \Delta f(1) etc.
  • In general, the k-th difference is: \Delta ^k f(n)=\Delta ^{k-1}f(n+1) - \Delta ^{k-1}f(n).
  • As a starting condition in the induction of differences, we can say that f(n) itself is the 0-th difference.

Difference table

We can construct the difference table for our sequence and observe its properties. In a difference table, the first column is n, which runs from 0 to whatever number of elements we have for the sequence. The second column is the values of the sequence at these n. Then come the first difference, the second difference and so on. For our sample sequence we get the table:

\[\begin{matrix} n & f(n) & \Delta f(n) & \Delta ^2 f(n) & \Delta ^3 f(n) \\ 0 & 0 & 1 & 3 & 0 \\ 1 & 1 & 4 & 3 & 0 \\ 2 & 5 & 7 & 3 & 0 \\ 3 & 12 & 10 & 3 & \\ 4 & 22 & 13 & & \\ 5 & 35 & & & \end{matrix}\]

Notice how at some point the differences become all-zero! We'll soon see why.

Obviously, we can construct such a table from only the column of f(n) - that's what we just did! A more interesting observation is that if we accept that all differences (columns) are 0 from a certain point, we can also construct this table from just the first row! For example, with f(0) and \Delta f(0) in hand, we know f(1); with \Delta f(0) and \Delta^2 f(0), we know \Delta f(1) etc. Try it!

Inferring the polynomial's degree from the table

Claim: if f(n) has degree k, then the k-th difference column in the table is constant.

Proof: this is a general k-th degree polynomial with coefficients a_k:

\[f(n)=a_k n^k + a_{k-1} n^{k-1} + \cdots + a_1 n + a_0\]

By definition of the first difference, if we expand the polynomial form and perform the subtraction per power of n:

\[\begin{align*} \Delta f(n) &= f(n+1)-f(n) = a_k (n+1)^{k} - a_k n^k + a_{k-1} (n+1)^{k-1} - a_{k-1} n^{k-1}+\cdots \\ &= \sum_{j=0}^{k} a_j(n+1)^j-a_j n^j \end{align*}\]

Using the binomial theorem we know that:

\[\begin{align*} (n+1)^j &= \sum_{i=0}^{j} \binom{j}{i}n^{j-i} \cdot 1^{i} \\ &= n^j + \binom{j}{1}n^{j-1}+\binom{j}{2}n^{j-2}+ \cdots \end{align*}\]

Therefore:

\[(n+1)^j - n^j = \binom{j}{1}n^{j-1}+\binom{j}{2}n^{j-2}+ \cdots\]

Now if we look at the sum we got for \Delta f(n) again, we'll notice that in each term, the j-th power of n gets canceled out. This means that the highest power of n in \Delta f(n) is going to be k-1.

We can similarly show that in \Delta^2 f(n), the highest power of n is going to be k-2. Therefore, the k-th difference \Delta^k f(n) will be constant \blacksquare.

Two observations for extra credit:

  1. Note that the claim goes one way - if f(n) is a k-th degree polynomial, the k-th difference is constant. What we observe in the table is the k-th difference is constant, so can we infer that the function is a k-th degree polynomial? Not in the general case! The sequence could be generated by some higher-degree polynomial, or by a completely different kind of function. That said, since we assume f(n) is polynomial and are seeking to find the simplest (lowest degree) one, this inference is valid.
  2. Did you notice the equivalence to derivatives of polynomials? The k-th difference of a polynomial of degree k is constant... but the same is true for the k-th derivative! This is not by chance - since we have a discrete domain, differences play largely the same role as derivatives for continuous functions. This isn't a rigorous proof - but think about the definition of derivatives (the one with the limit) - what do you get when you take \Delta x (also sometimes called h) and set it to 1?

Finding the coefficients with the Newton polynomial

Polynomial interpolation can fit any N points (with distinct x values) with a N-1 degree polynomial. One way of finding such a polynomial was discovered by Isaac Newton and is called the Newton polynomial.

Our problem of finding a polynomial that generates a given set of points can be reduced to this interpolation problem, since we've just figured out the degree of the generating polynomial! Looking at the difference table, we've found when the difference becomes constant, and that gives us the polynomial's degree k. So all we need is the first k+1 points.

Here's how to develop the Newton polynomial from scratch; we'll start with the first few coefficients and will then generalize for any k.

The Newton polynomial for our set of forward differences can be expressed as follows:

\[f(n) = b_0 + b_1 n + b_2 n (n-1) + \cdots + b_k n(n-1)(n-2)\cdots (n-k+1)\]

This polynomial is constructed in a clever way; notice that for any p, when we calculate f(p) all the elements starting with b_{p+1} will be multiplied by zero and vanish. This helps us determine this polynomial's coefficients in a gradual manner.

We'll start with the point (0, f(0)), substituting it into the Newton polynomial:

\[f(0) = b_0\]

This gives us the first coefficient b_0. Next, let's look at (1, f(1)):

\[f(1) = b_0 + b_1 \cdot 1 = b_0 + b_1\]

Since we know that b_0=f(0), we can infer that b_1=f(1)-f(0). Another way to express that is b_1=\Delta f(0).

Continuing to (2, f(2)):

\[\begin{align*} f(2) &= b_0 + b_1 \cdot 2 + b_2 \cdot 2 \cdot 1 \\ &= b_0 + 2 b_1 + 2! b_2 \\ &= f(0) + 2 \Delta f(0) + 2! b_2 \end{align*}\]

We've substituted the values of b_0 and b_1 that we've found earlier; let's solve for b_2 now:

\[\begin{align*} b_2 &= \frac{f(2) - f(0) - 2 \Delta f(0)}{2!} \\ &= \frac{f(2) - f(1) + f(1) - f(0) - 2 \Delta f(1)}{2!} \\ &= \frac{\Delta f(1) + \Delta f(0) - 2 \Delta f(0)}{2!} \\ &= \frac{\Delta f(1) - \Delta f(0)}{2!} \end{align*}\]

The last line's numerator is - by definition - \Delta^2 f(0):

\[b_2 = \frac{\Delta^2 f(0)}{2!}\]

We can keep going with this (feel free to do (3, f(3)) as an exercise), but the emerging generalization is that:

\[b_i = \frac{\Delta^i f(0)}{i!}\]

And a concise way to write Newton's polynomial for forward differences is:

\[f(n) = f(0) + \sum_{i=1}^{k} \frac{\Delta^i f(0)}{i!}g_i(n)\]

Where g_i(n) is:

\[\prod_{j=0}^{i-1} (n-j)\]

Note that we only use the differences for f(0), meaning that we need just the first row of the difference table! Let's try it for the pentagonal numbers example.

First, we've determined that \Delta^2 f(n) is a constant, so the degree of the polynomial is 2. We only need to calculate until i=2 in the sum:

\[f(n)=f(0)+\frac{\Delta f(0)}{1!} n + \frac{\Delta^2 f(0)}{2!} n(n-1)\]

Substituting the values from the difference table we have for n=0, we get:

\[\begin{align*} f(n)&=0+n+\frac{3}{2}n(n-1) \\ &=n + \frac{3n^2-3n}{2} \\ &=\frac{3 n^2}{2} -\frac{n}{2} \end{align*}\]

Which is exactly what we expected!

Another example

Let's work through another example, taking the sequence -8, -12, -6, 16, 60... We'll start by constructing the difference table:

\[\begin{matrix} n & f(n) & \Delta f(n) & \Delta ^2 f(n) & \Delta ^3 f(n) \\ 0 & -8 & -4 & 10 & 6 \\ 1 & -12 & 6 & 16 & 6 \\ 2 & -6 & 22 & 22 & \\ 3 & 16 & 44 & & \\ 4 & 60 & & & \\ \end{matrix}\]

The difference \Delta^3 f(n) appears to be constant, so we can generate this sequence with a degree 3 polynomial. Let's use the first line of the table to construct the Newton polynomial for it.

\[\begin{align*} f(n)&=f(0)+\frac{\Delta f(0)}{1!} n + \frac{\Delta^2 f(0)}{2!} n(n-1) + \frac{\Delta^3 f(0)}{3!} n(n-1)(n-2)\\ &= -8 + \frac{-4}{1}n + \frac{10}{2}n(n-1) + \frac{6}{6}n(n-1)(n-2) \\ &= -8-4n+5n^2-5n+n^3-3n^2+2n \\ &= n^3+2n^2-7n-8 \end{align*}\]

Verifying that this polynomial generates our sequence as its first 5 elements is an easy exercise.

Note that given 5 elements, we can always find a 4th-degree polynomial fitting it. Here we found a 3rd-degree one, though, leveraging the technique of differences. This becomes more acute if we have more elements in the sequence - using differences it's often possible to find significantly simpler polynomials.

Recap

For an integer sequence, if this sequence is generated by a polynomial we can figure out which polynomial it is - given enough elements. We start by constructing a difference table and noticing if a column becomes constant from some point on. This tells us the degree of the generating polynomial. With that in hand, we can use the Newton polynomial to discover a polynomial that generates the sequence.

Resources

Appendix: Another way to find the polynomial's coefficients

Here's another way to discover the coefficients of the polynomial; the following discusses the coefficients of the highest power, but can be generalized to other powers as well. Using Newton's polynomial is simpler, though, so this is just an appendix for some extra practice in manipulating such polynomials.

Claim: if f(n) has degree k, its coefficient a_k is equal to \frac{\Delta^k f(n)}{k!}. Since we've proven that \Delta^k f(n) is a constant, we can find the precise value of a_k this way.

Proof: Let's go back to the sum formulation of \Delta f(n) from the previous proof.

\[\sum_{j=0}^{k} a_j(n+1)^j-a_j n^j\]

Expanding the binomial, we get:

\[\sum_{j=0}^{k} a_j \left [ \binom{j}{1}n^{j-1}+\binom{j}{2}n^{j-2}+ \cdots \right ]\]

The highest power of n here is k-1; its coefficient comes only from the first term of the binomial expansion for j=k, and is equal to k a_k. In order not to deal with long sums, let's just focus on the highest-degree term in \Delta f(n):

\[\Delta f(n) = k a_k n^{k-1} + \gamma\]

Where \gamma represents other elements with lower powers of k, so we don't care about them for this discussion [2]. Let's move on to the next difference:

\[\begin{align*} \Delta^2 f(n) &= \Delta f(n+1) - \Delta f(n) \\ &= k a_k (n+1)^{k-1} + \gamma_1 - k a_k n^{k-1} + \gamma_2 \\ &= k a_k \left [ n^{k-1} + \binom{k-1}{1} n^{k-2} + \cdots \right ] - k a_k n^{k-1} + \gamma \\ &= k a_k \left [ \binom{k-1}{1} n^{k-2} + \cdots \right ] + \gamma \\ &= k(k-1)a_k n^{k-2} + \gamma \end{align*}\]

We've just found the coefficient of the highest power of n in \Delta^2 f(n). It's clear that if we continue doing this, \Delta^k f(n) will be:

\[\Delta^k f(n)= a_k k(k-1)(k-2)\cdots 1 = a_k k!\]

In other words, a_k = \frac{\Delta^k f(n)}{k!}, meaning that we can know the coefficient of the highest power of n in f(n) from the k-th difference \blacksquare.

[1]Naturally, there's an infinitude of potential functions that generate this sequence; it's more precise to say we're looking for the simplest (lowest degree) polynomial that would do this.
[2]In the following calculation, we're playing loose with the \gammas to represent "anything with lower powers of n that we don't care about".

For comments, please send me an email.