MATH/APPM 4650
Class notes, February 11

Divided differences

Recall that the Lagrange polynomial, given three points {x₀, x₁, x₂}, is P₂(x) = f(x₀) [(x-x₁)(x-x₂)] / [(x₀-x₁)(x₀-x₂)] + f(x₁) [(x-x₀)(x-x₂)] / [(x₁-x₀)(x₁-x₂)] + f(x₂) [(x-x₀)(x-x₁)] / [(x₂-x₀)(x₂-x₁)].
Although it's easy to write an algorithm to find this, there are a few problems with this form:

1. Having the previous Lagrange polynomial doesn't much help us get the next one.
2. The polynomial is not in Horner (nested) form. So it will take a long time to evaluate it at any particular x.
3. It's not obvious how many terms one would need to get a good approximation.
4. As all the points approach x₀, the Lagrange polynomial should approach the Taylor polynomial about x₀. That's not obvious from this form.

An alternative formula for the same thing was discovered by Newton, using divided differences.

P₂(x) = f[x₀] + f[x₀, x₁] (x-x₀) + f[x₀, x₁, x₂] (x-x₀) (x-x₁).
The divided differences are computed recursively using f[x₀] = f(x₀)
f[x₀, x₁] = [ f[x₁] - f[x₀] ] / [ x₁ - x₀ ]
f[x₀, x₁, x₂] = [ f[x₁, x₂] - f[x₀,x₁] ] / [ x₂ - x₀ ]
etc.

We can make a table out of these:

dd[0,0]
dd[1,0]	dd[1,1]
dd[2,0]	dd[2,1]	dd[2,2]
dd[3,0]	dd[3,1]	dd[3,2]	dd[3,3]

Here the general recursion formula is dd[j,i] = ( dd[j,i-1] - dd[j-1,i-1] ) / ( x[j] - x[j-i] ). So each term depends both on the term to the left, and the term above that one.
Hence each term along the diagonal (the one we want) depends on the entire triangle formed to the left and above it.


Here is a pseudocode algorithm to compute the divided differences, given an array of values "xvals."

dd[0,0] = f(xvals[0])
for j from 1 to Max do
   dd[j,0] = f(xvals[j])
   for i from 1 to j do
      dd[j,i] = (dd[j,i-1]-dd[j-1,i-1]) / (xvals[j]-xvals[j-i])
   end
end


The polynomial is then P_n(x) = dd[0,0] + dd[1,1] (x - x₀) + dd[2,2] (x - x₀) (x - x₁) + ... + dd[n,n] (x - x₀) (x - x₁) ... (x - x_n-1).
Obviously to compute this for any particular value of x, we should use Horner's method.
This algorithm computes the j^th interpolating polynomial at x, given the array dd[0,0] through dd[j,j] and the array xvals[0] through xvals[j-1].

poly[j](x) = dd[j,j]
for k from 1 to j do
   poly[j](x) = (x-xvals[j-k])*poly[j](x) + dd[j-k,j-k]
end do



Practically, these programs rely on being able to store the coefficients as doubly-indexed arrays.

In Maple, you could first declare the arrays you need as follows. You list 0 (the first index), two periods, and Max (the last index).
dd := array(0..Max, 0..Max);
xvals := array(0..Max);
and then access the terms using
dd[j,i]
xvals[i]


In MATLAB, things are a bit different. Although you can work with arrays without declaring them first, in the long run it's better to declare them.
You declare an array by filling up an array of the desired size with zeros. Also, you cannot start with index zero (which is sometimes awkward).
dd = zeros(Max+1, Max+1)
xvals = zeros(Max+1)
You access the terms using
dd(j,i)
xvals(i)


In Mathematica, you declare an array as follows.
Array[dd, {Max,Max}, {0,0}]
Array[xvals, Max, 0]
The terms are accessed using
dd[j,i]
xvals[i]

---

We then come to a few important questions.

1. How should we choose the points x₀, x₁, ..., x_n so as to minimize the error?
2. How many terms will we need to take to get a good approximation?