Polynomial interpolation and Approximation
Dr. John Travis

Weierstrass Approximation Theorem:  If feC[a,b], then there exists a polynomial P(x) such that P(x) is arbitrarily close to f(x) on [a,b].

Pf:  Wolog, assume [a,b] = [0,1].
Indeed, if not we can always apply the theorem to the function
g(x) = f(a+x(b-a))
to get g(0)=f(a) and g(1)=f(b).
Further, wolog assume that f(0) = f(1) = 0.
Indeed, if not, we can apply the theorem to the function
h(x)  = f(x) - f(a) - [f(b) - f(a)]*(x-a)/(b-a)
to get h(a) = h(b) = 0.
In either case, we can then convert the resulting polynomial P(x) back to the original settings.

Define Qn(x) = cn(1 - x2)n, for n=1, 2, 3, ... where cn is chosen so that

1 =
 ó
õ		 1

-1

 Qn(x) dx 
Lemma:  (1 - x2)n> (1 - nx2), for 0< x < 1.
Pf of Lemma:  Consider h(x) = (1 - x2)n - (1 - nx2).
Notice h(0) = 0 and h'(x) = 2nx(1 - (1-x2)n-1) > 0 for 0 < x < 1.
        => h(x) is increasing with a minimum of zero
        => h(x) is non-negative for 0 < x < 1.
Note,
  ó
õ		 1

-1

 (1-x2)n dx =
2  ó
õ		 1

0

 (1-x2)n dx >
2
 ó
õ		 1

1/n1/2

 (1-x2)n dx >
2
 ó
õ		 1

1/n1/2

 (1 - nx2) dx =
                                    = 4/ (3 Ön )
                                    > 1 / Ön.
Therefore, cn < Ön.
Now, for d < |x| < 1,
Qn(x) = cn(1 - x2)n < Ön (1 - d2)n
and so Qn(x) -> 0 uniformly for d< |x| < 1.
Consider
Pn(x) = ó
õ		 1-x

-x

 f(x+t)Qn(t) dt  =
ó
õ		 1

0

 f(t)Qn(t-x) dt
using a simple change of variable.
Notice, the last integral is clearly a polynomial in x since the variable t which involves the function f(t) disappears with the integration and the only place the x appears is in the Qn(t-x) term, where Qn is a polynomial.
Thus {Pn} is a sequence of polynomials.
Since f is a continuous function, given e > 0, choose d > 0 such that
|(x+t) - x| < d => |f(x+t) - f(x)| < e/2.
Let M = sup |f(x)|.
Then, for 0 < x < 1, noting that we have chosen Qn so that it integrates to 1,
| Pn(x) - f(x) | =
| ó
õ		 1

-1

 f(x+t)Qn(t)dt  - f(x) ó
õ		  1

-1

 Qn(t) dt |
< ó
õ		 1

-1

 |f(x+t)-f(x)|Qn(t)dt
, by pulling the absolute value on the inside
< 2M ó
õ		 -d

-1

 Qn(t)dt  + 
e/2 ó
õ		  d

-d

 Qn(t)dt +
2M  ó
õ		 1

d

 Qn(t)dt
                      < 4M Ön (1 - d2)n + e/2, noting that Qn(t) is even on the first integral.
                        < e, for all n large enough...say N
For this N, PN(x) is a finite degree polynomial within e distance from f(x), as desired.
Remark: This is an existence theorem. In practice, finding the polynomial might be very difficult and of very high degree.

Taylor Polynomials:  Given feCn+1[a,b] and some value c with a<c<b,

f(x) = f(c) + f'(c)(x-c) + f''(c) (x-c)2/2! + f'''(c) (x-c)3/3! + ... f(n)(c) (x-c)n/n! + Rn(x)
where Rn(x) = f(n+1)(z) (x-c)n+1/(n+1)!, where z is some point between c and x.
Pf:  Write the polynomial part above as Pn(x).  So we want to prove f(x) = Pn(x) + Rn(x).
Certainly for x=c, the statement holds.
Assume x is not equal to c and define Rn(x) = f(x) - Pn(x).
Define
g(t) = f(x) - f(t) - f'(t)(x-t) - f''(t) (x-t)2/2! - f'''(t) (x-t)3/3! - ... f(n)(t) (x-t)n/n! - Rn(x) (x-t)n+1/(x-c)n+1.
Notice,
g'(t) = - f(n+1)(t) (x-t)n/n! + (n+1)Rn(x) (x-t)n/(x-c)n+1,
for all t between c and x.
Moreover, for the fixed x, g(c) = 0 and g(x) = 0.
Thus, by the Mean Value Theorem (actually, the special case known as Rolle's Theorem), there is a number z between c and x such that g'(z) = 0.
Substituting this into the expression above yields
0 = - f(n+1)(z) (x-z)n/n! + (n+1)Rn(x) (x-z)n/(x-c)n+1,
and by solving yields
Rn(x) = f(n+1)(z) (x-c)n+1/(n+1)!
Finally, since g(c)=0, plugging this back into the definition of g(c) yields
0 = f(x) - f(c) - f'(t)(x-c) - f''(c) (x-c)2/2! - f'''(c) (x-c)3/3! - ... f(n)(c) (x-c)n/n! - Rn(x) (x-c)n+1/(x-c)n+1,
from which the result easily follows.
These are created by using Taylor's theorem and throwing away the remainder term. Bad Ex: Approximate f(x)=1/x about c=1 using higher and higher Taylor's polynomials. The approximation gets worse and worse as the degree of the polynomial gets larger.

LaGrange Polynomials: Given n+1 values of the function (x0,f(x0)) , (x1,f(x1)) , ... , (xn,f(xn)), find the polynomial of degree at most n which passes through all these points.

Derivation:

  1. Utilize a general nth degree polynomial and plug in the data values to get a system of equations in the coefficients.
  2. Use Lagrange cardinal formulas. Writing these using a Horner-type arrangement gives a faster implementation.
  3. Newton's Divided Differences...generally used when actually computing the polynomial.
Ex: Approximate f(x)=1/x using x0=1, x1=2, etc. for a few values. This gives a much better approximation than Taylors did.

Lagrange error bound Theorem:  Let x0,x1, , ... , xn be distinct numbers in the interval [a,b] and let feCn+1[a,b].  If P(x) is the Lagrange Interpolating Polynomial, then for each x in [a,b], there exists a number z in (a,b) such that

f(x) = P(x) +  f(n+1)(z) (x-x0) (x-x1) ... (x-xn)/(n+1)!.
Pf:  Easily, if x = xk, for one of the given x values, then f(x) = P(x) and the remainder is zero.
So, suppose x is a fixed value different that all of the given data values {xk}.  Consider
g(t) = f(t) - P(t) - [f(x) - P(x)] (t-x0) (t-x1) ... (t-xn) / [(x-x0) (x-x1) ... (x-xn)]
Certainly g has n+1 continuous derivatives.
Now, for t = xk,  g(xk) = 0.
Moreover, g(x) = 0.
Thus, g(t) has n+2 zeros in the interval [a,b].
So, by the generalized Rolle's Theorem, there exists z in (a,b) such that
g(n+1)(z) = 0, or
0 = f(n+1)(z) - P(n+1)(z) - [f(x) - P(x)] D(n+1)[ (t-x0) (t-x1) ... (t-xn) ] / [(x-x0) (x-x1) ... (x-xn)].
Note, P(t) is a polynomial of degree n so its (n+1)st derivative will be zero.
Further, the term (t-x0) (t-x1) ... (t-xn) is a polynomial of degree n+1 so its (n+1)st derivative will be a constant.  Indeed, the leading term in the product is tn+1 whose (n+1)st derivative is (n+1)!
Therefore, the last equation becomes
0 = f(n+1)(z) - [f(x) - P(x)] (n+1)! / [(x-x0) (x-x1) ... (x-xn)].
Solving for f(x) yields the desired result.
Lagrange Uniqueness Theorem: Any two polynomials of the same degree n which agree at n+1 points are equal.
Pf: Suppose f(x), g(x) are of degree n and agree at n+1 points.
Then, r(x)=f(x)-g(x) is of degree at most n with n+1 roots gives r(x)=0.
Ex: Construct a parabola passing through (1,4), (-1,0) and (2,9). Then p(x)= ax2 + bx + c and by plugging in the points,
4 = a + b + c
0 = a - b + c
9 = 4a + 2b + c
which is a 3x3 system of equations to solve. If the polynomial p(x) were of higher degree, then hard to solve.

But if we write p(x)=c + b(x-1) + a(x-1)(x+1), then plugging in the points yields (solving as we go along)

4 = c
0 = c - 2b, or b=2
9 = c + b + 3a, or a=1,
which is very easy to solve successively. This second form is called Newton's form.

Newton's Divided Differences: Successively build up the polynomial by finding the coefficients as above.

(x0,f0) gives p0(x) = f0, passes through the one point.
(x1,f1) gives p1(x) = p0(x) + (x-x0)b1, passes through the first two points if b1=(f1-f0)/(x1-x0).
(x2,f2) gives p2(x) = p1(x) + (x-x0)(x-x1)b2, passes through the first three points if b2=?.
Continuing in this manner shows the coefficients bi are given by the divided difference table, page 63, etc.
f[xj,...,xk] = (f[xj+1,...,xk] - f[xj,...,xk-1]) / (xk - xj), for j<k.
If x0 < x1 < ... < xn, this is called Newton's Forward Divided Difference Polynomial.
If x0 > x1 > ... > xn, this is called Newton's Backward Divided Difference Polynomial.

If Dx is constant, the Divided Difference Formulas can be written in a closed form involving the binomial coefficents.
 

Program Subroutine To Determine the Newton's Divided Differences:

cj=f(xj) for j=0,1,...,n
for k= 1:n
    for j=n:-1:k
        cj:=(cj-cj-1)/(xj-xj-k)
    end
end
Then, the Newton's Interpolation polynomial is given by (with the cj from above):
    p(x)    = c0 + c1(x-x0) + c2(x-x0)(x-x1) + c3(x-x0)(x-x1)(x-x2) + ... + cn(x-x0)(x-x1)...(x-xn-1).
                = c0 + (x-x0){c1 + (x-x1) {c2 + (x-x2) {c3 + ... + (x-xn-2) {cn-1 + cn(x-xn-1) }...} } }.

Homework:  Go to some busy local store and collect some data regarding the number of customers that enter the store during a given amount of time.  Carefully count the cummulative number that enter the store within 2, 5, 7 and 10 minutes (or other similar, non-uniform time steps but no more than five points).  Using this data:

Piecewise Interpolation: Instead of using one polynomial (possibly of very high degree) for the entire interval of interest, use several polynomials each defined on small intervals. Then, to obtain the y-value for a given x-value, one must first determine which formula to use. Then, evaluate that formula for the given x-value. If x0 < x1 < ... < xn, and x0 < x < xn, then one can determine which interval x belongs to by considering the iteration:
k = 0
Repeat
    k = k+1
    (Use the kth formula...)
Until x < xk
Return k
Notice, if we desire to evaluate this piecewise approximation for a fixed number of equally spaced points in each sub-interval, we can evaluate the basis polynomials once and store the results in vectors. These then can be reused for each interval.

Linear Splines:  Piece linear segments Lk(x) = ak + bk(x-xk-1) together to interpolate y-values only at the ends:

Cubic Splines: Piece cubics Sk(x) = ak + bk(x-xk-1) + ck(x-xk-1)2 + dk(x-xk-1)3, for xk-1 < x < xk, together in such a way that the resulting approximation is continuous and has 2 continuous derivatives (no visible kinks). So, we impose the conditions: The derivation for these equations...

Hence, we have 4n unknowns (the vectors of coefficients a, b, c and d) and 4n-2 equations. By adding two additional equations, we can uniquely solve for the unknown vectors.

End Conditions:

  ó
õ		 1

-1

 (1-x2)n dx =