Divided Difference

Polynomial Approximation: Practical Computations
Polynomial Approximation, Higher Order Matching

Beyond Hermite Interpolatory Polynomials
Numerical Analysis and Computing

Lecture Notes #5 Interpolation and Polynomial
Approximation
Divided Differences, and Hermite Interpolatory Polynomials
Joe Mahaffy,
hmahaffy@math.sdsu.edui
Department of Mathematics
Dynamical Systems Group
Computational Sciences Research Center
San Diego State University

San Diego, CA 92182-7720
http://www-rohan.sdsu.edu/jmahaffy
Spring 2010
Joe Mahaffy, hmahaffy@math.sdsu.edui
#5 Interpolation and Polynomial Approximation (1/41)

Outline

Representing Polynomials
Divided Differences
Different forms of Divided Difference Formulas

Osculating Polynomials
Hermite Interpolatory Polynomials
Computing Hermite Interpolatory Polynomials

Recap and Lookahead

Previously:
Nevilles Method to successively generate higher degree polynomial approximations at a specific point. If we need to compute the polynomial at many points, we have to re-run Nevilles method for each point.
O(n2 ) operations/point.
Algorithm: Nevilles Method
To evaluate the polynomial that interpolates the n + 1 points (xi , f (xi )), i = 0, . . . , n
at the point x:
1. Initialize Qi,0 = f (xi ).
2. FOR i = 1 : n
FOR j = 1 : i
(x xij )Qi,j1 (x xi )Qi1,j1
Qi,j =
xi xij
END
END
3. Output the Q-table.

Recap and Lookahead

Next:
Use divided differences to generate the polynomials themselves.
The coefficients of the polynomials. Once we have those, we can

quickly (remember Horners method?) compute the polynomial in
any desired points. O(n) operations/point.
Algorithm: Horners Method
Degree n; coefficients a0 , a1 , . . . , an ; x0
Input:
Output:
y = P(x0 ), z = P (x0 ).
1.
Set y = an , z = an
2.
For j = (n 1), (n 2), . . . , 1

Set y = x0 y + aj , z = x0 z + y
3.
Set y = x0 y + a0
4.
Output (y , z)
5.
End program

Divided Differences
If Pn (x) is the nth degree polynomial that agrees with f (x) at the
points {x0 , x1 , . . . , xn }, then we can (for the appropriate constants
{a0 , a1 , . . . , an }) write:
Pn (x) =
a0 + a1 (x x0 ) + a2 (x x0 )(x x1 ) +
+ an (x x0 )(x x1 ) (x xn1 )

Divided Differences
If Pn (x) is the nth degree polynomial that agrees with f (x) at the
points {x0 , x1 , . . . , xn }, then we can (for the appropriate constants
{a0 , a1 , . . . , an }) write:
Pn (x) =
a0 + a1 (x x0 ) + a2 (x x0 )(x x1 ) +
+ an (x x0 )(x x1 ) (x xn1 )
Note that we can evaluate this Horner-style, by writing

Pn (x) =
a0 + (x x0 ) (a1 + (x x1 ) (a2 +
+ (x xn2 ) (an1 + an (x xn1 )))) ,
so that each step in the Horner-evaluation consists of a

subtraction, a multiplication, and an addition.

Divided Differences
Finding the Constants {a0 , a1 , . . . , an }
Just Algebra
Given the relation

Pn (x) =
a0 + a1 (x x0 ) + a2 (x x0 )(x x1 ) +
+ an (x x0 )(x x1 ) (x xn1 )
at x0 : a0 = Pn (x0 ) = f (x0 ).
at x1 : f (x0 ) + a1 (x1 x0 ) = Pn (x1 ) = f (x1 )
f (x1 ) f (x0 )
.
x1 x 0
f (x1 ) f (x0 )
f (x2 ) f (x0 )
.
at x2 : a2 =
(x2 x0 )(x2 x1 ) (x2 x0 )(x1 x0 )
a1 =
This gets massively ugly fast! We need some nice clean

notation!

Divided Differences
Sir Isaac Newton to the Rescue: Divided Differences

Zeroth Divided Difference:
f [xi ] = f (xi ).
First Divided Difference:
f [xi , xi+1 ] =
f [xi+1 ] f [xi ]
.
xi+1 xi
Second Divided Difference:

f [xi , xi+1 , xi+2 ] =
f [xi+1 , xi+2 ] f [xi , xi+1 ]

.
xi+2 xi
k th Divided Difference:
f [xi , xi+1 , . . . , xi+k ] =
f [xi+1 , xi+2 , . . . , xi+k ] f [xi , xi+1 , . . . , xi+k1 ]

.
xi+k xi

Divided Differences
The Constants {a0 , a1 , . . . , an } Revisited

We had
at x0 : a0 = Pn (x0 ) = f (x0 ).
at x1 : f (x0 ) + a1 (x1 x0 ) = Pn (x1 ) = f (x1 )
f (x1 ) f (x0 )
.
x1 x 0
f (x1 ) f (x0 )
f (x2 ) f (x0 )
.
at x2 : a2 =
(x2 x0 )(x2 x1 ) (x2 x0 )(x1 x0 )
a1 =
Clearly:
a0 = f [x0 ],
a1 = f [x0 , x1 ].
We may suspect that a2 = f [x0 , x1 , x2 ], that is indeed so (a little

bit of careful algebra will show it), and in general
ak = f[x0 , x1 , . . . , xk ].

Divided Differences
Algebra: Chasing down a2 = f [x0 , x1 , x2 ]

a2
f (x2 ) f (x0 )
f (x1 ) f (x0 )
(x2 x0 )(x2 x1 )
(x2 x1 )(x1 x0 )
(f (x2 ) f (x0 ))(x1 x0 ) (f (x1 ) f (x0 ))(x2 x0 )

(x2 x0 )(x2 x1 )(x1 x0 )
(x1 x0 )f (x2 ) (x2 x0 )f (x1 ) + (x2 x0 x1 + x0 )f (x0 )

(x2 x0 )(x2 x1 )(x1 x0 )
(x1 x0 )f (x2 ) (x1 x0 + x2 x1 )f (x1 ) + (x2 x1 )f (x0 )

(x2 x0 )(x2 x1 )(x1 x0 )
(x1 x0 )(f (x2 ) f (x1 )) (x2 x1 )(f (x1 ) f (x0 ))

(x2 x0 )(x2 x1 )(x1 x0 )
(f (x2 ) f (x1 ))
(f (x1 ) f (x0 ))
(x2 x0 )(x2 x1 )
(x2 x0 )(x1 x0 )
f [x0 , x1 ]
f [x1 , x2 ]
= f[x0 , x1 , x2 ]
x2 x 0
x2 x 0
(!!!)

Divided Differences
Newtons Interpolatory Divided Difference Formula

Hence, we can write
Pn (x) = f [x0 ] +
n
X
k=1
"
f [x0 , . . . , xk ]
k1
Y
m=0
(x xm ) .
Pn (x) = f [x0 ] +
f [x0 , x1 ](x x0 ) +
f [x0 , x1 , x2 ](x x0 )(x x1 ) +
f [x0 , x1 , x2 , x3 ](x x0 )(x x1 )(x x2 ) +
This expression is known as Newtons Interpolatory Divided
Difference Formula.

Divided Differences
Computing the Divided Differences (by table)

x
x0
x1
x2
x3
x4
x5
f(x)
f [x0 ]
1st Div. Diff.

f [x0 , x1 ] =
f [x1 ]f [x0 ]
x1 x0
f [x1 , x2 ] =
f [x2 ]f [x1 ]
x2 x1
f [x2 , x3 ] =
f [x3 ]f [x2 ]
x3 x2
f [x3 , x4 ] =
f [x4 ]f [x3 ]
x4 x3
f [x4 , x5 ] =
f [x5 ]f [x4 ]
x5 x4
f [x1 ]
f [x2 ]
f [x3 ]
f [x4 ]
f [x5 ]
2nd Div. Diff.
f [x0 , x1 , x2 ] =
f [x1 ,x2 ]f [x0 ,x1 ]

x2 x0
f [x1 , x2 , x3 ] =
f [x2 ,x3 ]f [x1 ,x2 ]

x3 x1
f [x2 , x3 , x4 ] =
f [x3 ,x4 ]f [x2 ,x3 ]

x4 x2
f [x3 , x4 , x5 ] =
f [x4 ,x5 ]f [x3 ,x4 ]

x5 x3
Note: The table can be extended with three 3rd divided

differences, two 4th divided differences, and one 5th divided
difference.

Divided Differences
Algorithm: Computing the Divided Differences

Algorithm: Newtons Divided Differences
Given the points (xi , f (xi )), i = 0, . . . , n.
Step 1: Initialize Fi,0 = f (xi ), i = 0, . . . , n
Step 2:
FOR i = 1 : n
FOR j = 1 : i
Fi,j1 Fi1,j1
Fi,j =
xi xij
END
END
Result: The diagonal, Fi,i now contains f [x0 , . . . , xi ].

Divided Differences
A Theoretical Result: Generalization of the Mean Value Theorem

Theorem (Generalized Mean Value Theorem)
Suppose that f C n [a, b] and {x0 , . . . , xn } are distinct number in
[a, b]. Then (a, b) :
f [x0 , . . . , xn ] =
f (n) ()
.
n!
For n = 1 this is exactly the Mean Value Theorem...

So we have extended to MVT to higher order derivatives!
What is the theorem telling us?
Newtons nth divided difference is in some sense an approximation to the nth derivative of f .

Divided Differences
Newton vs. Taylor...

Using Newtons Divided Differences...
PnN (x) = f [x0 ] + f [x0 , x1 ](x x0 ) +
f [x0 , x1 , x2 ](x x0 )(x x1 ) +
f [x0 , x1 , x2 , x3 ](x x0 )(x x1 )(x x2 ) +
Using Taylor expansion
PnT (x) = f (x0 ) + f (x0 )(x x0 ) +
1
2
2! f (x0 )(x x0 ) +
1
3
3! f (x0 )(x x0 ) +
It makes sense that the divided differences are approximating the
derivatives in some sense!

Divided Differences
Simplification: Equally Spaced Points

When the points {x0 , . . . , xn } are equally spaced, i.e.
h = xi+1 xi , i = 0, . . . , n 1,
we can write x = x0 + sh, x xk = (s k)h so that
Pn (x) = Pn (x0 + sh) =
n
X
s(s 1) (s k + 1)hk f [x0 , . . . , xk ].
k=0
Using the binomial coefficients,

s(s 1) (s k + 1)
s
=
k
k!
Pn (x0 + sh) = f[x0 ] +
n
X
s
k=1
k! hk f[x0 , . . . , xk ].
This is Newtons Forward Divided Difference Formula.


Divided Differences
Notation, Notation, Notation...

Another form, Newtons Forward Difference Formula is constructed
by using the forward difference operator :
f (xn ) = f (xn+1 ) f (xn )
using this notation:
1
f (x1 ) f (x0 )
= f (x0 ).
x1 x0
h
1
1 f (x1 ) f (x0 )
= 2 2 f (x0 ).
f [x0 , x1 , x2 ] =
2h
h
2h
f [x0 , x1 ] =
1
k f (x0 ).
k! hk
Thus we can write Newtons Forward Difference Formula
n
X
s
k f(x0 ).
Pn (x0 + sh) = f[x0 ] +
k
f [x0 , . . . , xk ] =
k=1

Divided Differences
Notation, Notation, Notation... Backward Formulas

If we reorder {x0 , x1 , . . . , xn } {xn , . . . , x1 , x0 }, and define the backward
difference operator :
f (xn ) = f (xn ) f (xn1 ),
we can define the backward divided differences:
f [xn , . . . , xnk ] =
1
k f (xn ).
k! hk
We write down Newtons Backward Difference Formula

n
X
k s
k f(xn ),
(1)
Pn (x) = f[xn ] +
k
k=1
where
s
k
= (1)k
s(s + 1) (s + k 1)
.
k!

Divided Differences
Forward? Backward? Im Confused!!!

x
x0
x1
x2
x3
f (x)
f [x0 ]
1st Div. Diff.

f[x0 , x1 ] =
f[x1 ]f[x0 ]
x1 x0
f [x1 , x2 ] =
f [x2 ]f [x1 ]
x2 x1
f [x2 , x3 ] =
f [x3 ]f [x2 ]
x3 x2
f [x3 , x4 ] =
f [x4 ]f [x3 ]
x4 x3
f [x1 ]
f [x2 ]
f [x3 ]
x4
f [x4 ]
x5
f [x5 ]
f[x4 , x5 ] =
2nd Div. Diff.
f[x0 , x1 , x2 ] =
f[x1 ,x2 ]f[x0 ,x1 ]

x2 x0
f [x1 , x2 , x3 ] =
f [x2 ,x3 ]f [x1 ,x2 ]

x3 x1
f [x2 , x3 , x4 ] =
f [x3 ,x4 ]f [x2 ,x3 ]

x4 x2
f[x3 , x4 , x5 ] =
f[x4 ,x5 ]f[x3 ,x4 ]

x5 x3
f[x5 ]f[x4 ]
x5 x4
Forward: The fwd div. diff. are the top entries in the table.
Backward: The bwd div. diff. are the bottom entries in the table.

Divided Differences
Forward? Backward? Straight Down the Center!

The Newton formulas works best for points close to the edge of
the table; if we want to approximate f (x) close to the center, we
have to work some more...
x
x2
f (x)
f [x2 ]
x1
f [x1 ]
x0
f[x0 ]
x1
f [x1 ]
x2
f [x2 ]
x3
f [x3 ]
1st Div. Diff.

f [x2 , x1 ]
f[x1 , x0 ]
f[x0 , x1 ]
f [x1 , x2 ]
f [x2 , x3 ]
2nd Div. Diff.

f [x2 , x1 , x0 ]
f[x1 , x0 , x1 ]
f [x0 , x1 , x2 ]
f [x1 , x2 , x3 ]
3rd Div. Diff.
f[x2 , x1 , x0 , x1 ]
f[x1 , x0 , x1 , x2 ]
f [x0 , x1 , x2 , x3 ]
4th Div. Diff.
f[x2 , x1 , x0 , x1 , x2 ]
f [x1 , x0 , x1 , x2 , x3 ]
We are going to construct Stirlings Formula a scheme using

centered differences. In particular we are going to use the blue
(centered at x0 ) entries, and averages of the red (straddling the x0
point) entries.

Divided Differences
Stirlings Formula Approximating at Interior Points

Assume we are trying to approximate f (x) close to the interior
point x0 :
Pn (x)
P2m+1 (x) = f [x0 ] + sh

+ s 2 h2 f [x1 , x0 , x1 ]
f [x1 , x0 ] + f [x0 , x1 ]
2
f [x2 , x1 , x0 , x1 ] + f [x1 , x0 , x1 , x2 ]
2
+ s 2 (s 2 1)h4 f [x2 , x1 , x0 , x1 , x2 ]
+ s(s 2 1)h3
+ ...
+ s 2 (s 2 1) (s 2 (m 1)2 )h2m f [xm , . . . , xm ]
+ s(s 2 1) (s 2 m2 )h2m+1
f [xm1 , . . . , xm ] + f [xm , . . . , xm+1 ]

2
If n = 2m + 1 is odd, otherwise delete the last two lines.


Divided Differences
Summary: Divided Difference Formulas

Newtons Interpolatory Divided Difference Formula
Pn (x)
f [x0 ] + f [x0 , x1 ](x x0 ) + f [x0 , x1 , x2 ](x x0 )(x x1 ) +

f [x0 , x1 , x2 , x3 ](x x0 )(x x1 )(x x2 ) +
Newtons Forward Divided Difference Formula

Pn (x0 + sh) = f [x0 ] +
n
X
s
k! hk f [x0 , . . . , xk ]
k
k=1
Newtons Backward Difference Formula

Pn (x) = f [xn ] +
Reference: Binomial Coefficients

s
s(s 1) (s k + 1)
,
=
k
k!
n
s
X
(1)k
k f (xn )
k
k=1
s
s(s + 1) (s + k 1)
= (1)k
k
k!

Combining Taylor and Lagrange Polynomials

A Taylor polynomial of degree n matches the function and its
first n derivatives at one point.


A Lagrange polynomial of degree n matches the function values
at n + 1 points.


at n + 1 points.
Question: Can we combine the ideas of Taylor and Lagrange to
get an interpolating polynomial that matches both the
function values and some number of derivatives at multiple points?


at n + 1 points.
Answer: To our euphoric joy, such polynomials exist! They are
called Osculating Polynomials.


at n + 1 points.
Answer: To our euphoric joy, such polynomials exist! They are
called Osculating Polynomials.
The Concise Oxford Dictionary:
Osculate 1. (arch. or joc.) kiss. 2. (Biol., of species, etc.) be related through
intermediate species etc., have common characteristics with another or with each
other. 3. (Math., of curve or surface) have contact of higher than first order with,
meet at three or more coincident points.

In Painful Generality
Given (n + 1) distinct points {x0 , x1 , . . . , xn } [a, b], and non-negative

integers {m0 , m1 , . . . , mn }.
Notation: Let m = max{m0 , m1 , . . . , mn }.
The osculating polynomial approximation of a function f C m [a, b]
at xi , i = 0, 1, . . . , n is the polynomial (of lowest possible order) that
agrees with
{f (xi ), f (xi ), . . . , f (mi ) (xi )} at xi [a, b], i.
The degree of the osculating polynomial is at most
M =n+
n
X
mi .
i=0
In the case where mi = 1, i the polynomial is called a Hermite

Interpolatory Polynomial.

The Existence Statement
If f C 1 [a, b] and {x0 , x1 , . . . , xn } [a, b] are distinct, the unique

polynomial of least degree ( 2n + 1) agreeing with f (x) and f (x) at
{x0 , x1 , . . . , xn } is
H2n+1 (x) =
n
X
f(xj )Hn,j (x) +
j=0
where
n
X
n,j (x),
f (xj )H
j=0
h
i
Hn,j (x) = 1 2(x xj )Ln,j (xj ) L2n,j (x)
n,j (x) = (x xj )L2 (x),
H
n,j
and Ln,j (x) are our old friends, the Lagrange coefficients:
n
Y
Ln,j (x) =
i=0, i6=j
x xi
.
xj x i
Further, if f C 2n+2 [a, b], then for some (x) [a, b]

f (x) = H2n+1 (x) +
Qn
xi )2 (2n+2)
f
((x)).
(2n + 2)!
i=0 (x

Thats Hardly Obvious Proof Needed!
Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
(i,j is Kroneckers delta).

Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
n,j (xi ) = 0.
If follows that when i 6= j: Hn,j (xi ) = H

Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
n,j (xi ) = 0.
i
h
(
Hn,j (xj ) = 1 2(xj xj )Ln,j (xj ) 1 = 1
When i = j:
n,j (xj ) = (xj xj )L2 (xj ) = 0.
H
n,j

Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
n,j (xi ) = 0.
i
h
(
Hn,j (xj ) = 1 2(xj xj )Ln,j (xj ) 1 = 1
When i = j:
n,j (xj ) = (xj xj )L2 (xj ) = 0.
H
n,j
Thus, H2n+1 (xj ) = f(xj ).

Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
n,j (xi ) = 0.
i
h
(
Hn,j (xj ) = 1 2(xj xj )Ln,j (xj ) 1 = 1
When i = j:
n,j (xj ) = (xj xj )L2 (xj ) = 0.
H
n,j

(x)
Hn,j
=
=
[2Ln,j (xj )]L2n,j (x) + [1 2(x xj )Ln,j (xj )] 2Ln,j (x)Ln,j (x)
h
i
Ln,j (x) 2Ln,j (xj )Ln,j (x) + [1 2(x xj )Ln,j (xj )] 2(x)Ln,j

Recall:
Ln,j (xi ) = i,j =
0,
1
if i 6= j
if i = j
1 of 2
n,j (xi ) = 0.
i
h
(
Hn,j (xj ) = 1 2(xj xj )Ln,j (xj ) 1 = 1
When i = j:
n,j (xj ) = (xj xj )L2 (xj ) = 0.
H
n,j

(x)
Hn,j
=
=
[2Ln,j (xj )]L2n,j (x) + [1 2(x xj )Ln,j (xj )] 2Ln,j (x)Ln,j (x)
h
i
Ln,j (x) 2Ln,j (xj )Ln,j (x) + [1 2(x xj )Ln,j (xj )] 2(x)Ln,j
(x): H (x ) = 0 when i 6= j.
Since Ln,j (x) is a factor in Hn,j
n,j i

Proof, continued...
(x )
Hn,j
j
[2Ln,j (xj )] L2n,j (xj )

| {z }
2Ln,j (xj ) + 1 2 Ln,j (xj ) = 0
[1 2 (xj xj ) Ln,j (xj )] 2 Ln,j (xj ) Ln,j (xj )

| {z }
| {z }
0

Proof, continued...
(x )
Hn,j
j

| {z }
2Ln,j (xj ) + 1 2 Ln,j (xj ) = 0

| {z }
| {z }
0
i.e. Hn,j (xi ) = 0, i.

Proof, continued...
(x )
Hn,j
j

| {z }
2Ln,j (xj ) + 1 2 Ln,j (xj ) = 0

| {z }
| {z }
0
i.e. Hn,j (xi ) = 0, i.

(x)
H
n,j
=
=
L2n,j (x) + 2(x xj )Ln,j (x)Ln,j (x)

h
i
Ln,j (x) Ln,j (x) + 2(x xj )Ln,j (x)

Proof, continued...
(x )
Hn,j
j

| {z }
2Ln,j (xj ) + 1 2 Ln,j (xj ) = 0

| {z }
| {z }
0
i.e. Hn,j (xi ) = 0, i.

(x)
H
n,j
=
=

h
i
(xi ) = 0, since Ln,j (xi ) = i,j .

If i 6= j: H
n,j
i
h
If i = j: H (xj ) = 1 1 + 2(xj xj )L (xj ) = 1.

n,j
n,j

Proof, continued...
(x )
Hn,j
j

| {z }
2Ln,j (xj ) + 1 2 Ln,j (xj ) = 0

| {z }
| {z }
0
i.e. Hn,j (xi ) = 0, i.

(x)
H
n,j
=
=

h
i
(xi ) = 0, since Ln,j (xi ) = i,j .

If i 6= j: H
n,j
i
h
If i = j: H (xj ) = 1 1 + 2(xj xj )L (xj ) = 1.

n,j
n,j
Hence, H2n+1 (xi ) = f (xi ), i.


Uniqueness Proof
Assume there is a second polynomial G (x) (of degree
2n + 1) interpolating the same data.

Uniqueness Proof
Define R(x) = H2n+1 (x) G (x).

Uniqueness Proof
Then by construction R(xi ) = R (xi ) = 0, i.e. all the xi s are
zeros of multiplicity at least 2.

Uniqueness Proof
Q
This can only be true if R(x) = q(x) ni=0 (x xi )2 , for some
q(x).

Uniqueness Proof
Q
q(x).
If q(x) 6 0 then the degree of R(x) is 2n + 2, which is a
contradiction.

Uniqueness Proof
Q
q(x).
If q(x) 6 0 then the degree of R(x) is 2n + 2, which is a
contradiction.
Hence q(x) 0 R(x) 0 H2n+1 (x) is unique.

Main Use of Hermite Interpolatory Polynomials

One of the primary applications of Hermite Interpolatory Polynomials is
the development of Gaussian quadrature for numerical integration. (To
be revisited later this semester.)
The most commonly seen Hermite interpolatory polynomial is the cubic
one, which satisfies
H3 (x0 ) = f (x0 ),
H3 (x1 ) = f (x1 ),
H3 (x0 ) = f (x0 )
H3 (x1 ) = f (x1 ).
it can be written explicitly as

H3 (x)
=
+
0
1 + 2 xxx
1 x0
x
1 + 2 xx11x
0
ih
ih
x1 x
x1 x0
xx0
x1 x0
i2
i2
f (x0 ) + (x x0 )
f (x1 ) + (x x1 )
h
h
x1 x
x1 x0
xx0
x1 x0
i2
i2
f (x0 )
f (x1 ).
It appears in some optimization algorithms (see Math 693a, linesearch

algorithms.)

Computing from the Definition is Tedious!

However, there is good news: we can re-use the algorithm for Newtons
Interpolatory Divided Difference Formula with some modifications in
the initialization.
We double the number of points, i.e. let
{y0 , y1 , . . . , y2n+1 } = {x0 , x0 + , x1 , x1 + , . . . , xn , xn + }
Set up the divided difference table (up to the first divided differences),
and let 0 (formally), and identify:
f (xi ) = lim
0
f [xi + ] f [xi ]
,
to get the table [next slide]...

Hermite Interpolatory Polynomial using Modified Newton Divided Differences
y
y0 = x0
f(x)
f [y0 ]
y1 = x0
f [y1 ]
y2 = x1
f [y2 ]
y3 = x1
f [y3 ]
y4 = x2
f [y4 ]
y5 = x2
f [y5 ]
y6 = x3
f [y6 ]
y7 = x3
f [y7 ]
y8 = x4
f [y8 ]
y9 = x4
f [y9 ]
1st Div. Diff.

f [y0 , y1 ] =
f (y
2nd Div. Diff.

0)
f [y1 , y2 ]
f [y2 , y3 ] =
f (y
2)
f [y3 , y4 ]
f [y4 , y5 ] = f (y4 )
f [y5 , y6 ]
f [y6 , y7 ] =
f (y
6)
f [y7 , y8 ]
f [y8 , y9 ] =
f (y
8)
f [y0 , y1 , y2 ]
f [y1 , y2 , y3 ]
f [y2 , y3 , y4 ]
f [y3 , y4 , y5 ]
f [y4 , y5 , y6 ]
f [y5 , y6 , y7 ]
f [y6 , y7 , y8 ]
f [y7 , y8 , y9 ]
3rd Div. Diff.
f [y0 , y1 , y2 , y3 ]
f [y1 , y2 , y3 , y4 ]
f [y2 , y3 , y4 , y5 ]
f [y3 , y4 , y5 , y6 ]
f [y4 , y5 , y6 , y7 ]
f [y5 , y6 , y7 , y8 ]
f [y6 , y7 , y8 , y9 ]

H3 (x) revisited...
Old notation
H3 (x)
=
+
ih
i2
ih
i2
h
x1 x
x1 x
xx0
0
1 + 2 xxx
f
(x
)
+
1
+
2
f (x1 )
0
x
x
x
x
x
x
x
1
0
1
0
1
0
1
0
i2
i2
h
h
x
0
f (x0 ) + (x x1 ) xxx
f (x1 ).
(x x0 ) xx11x
0
1 x0
Divided difference notation

H3 (x)
f (x0 ) + f (x0 )(x x0 ) + f [x0 , x0 , x1 ](x x0 )2

+ f [x0 , x0 , x1 , x1 ](x x0 )2 (x x1 ).
Or with the y s...

H3 (x)
f (y0 ) + f (y0 )(x y0 ) + f [y0 , y1 , y2 ](x y0 )(x y1 )

+ f [y0 , y1 , y2 , y3 ](x y1 )(x y2 )(x y3 ).

H3 (x) Example
4
0.2
0.4
x0 = 0,
f (x0 ) = 0,
f (x0 ) = 4,
0.6
0.8
x1 = 1
f (x1 ) = 3,
f (x1 ) = 1
H3 (x) = 4x x 2 3x 2 (x 1)

H3 (x) Example
Not Very Pretty Computations
Example
x0 = 0; x1 = 1;
fv0 = 0; fpv0 = 4;
fv1 = 3; fpv1 = -1;
% This is the data
y0
y1
y2
y3
% Initializing the table
=
=
=
=
x0;
x0;
x1;
x1;
f0=fv0;
f1=fv0;
f2=fv1;
f3=fv1;
f01 = fpv0;
f12 = (f2-f1)/(y2-y1);
f23 = fpv1;
% First divided differences
f012 = (f12-f01)/(y2-y0);
f123 = (f23-f12)/(y3-y1);
% Second divided differences
f0123 = (f123-f012)/(y3-y0);
% Third divided difference
x=(0:0.01:1);
H3 = f0 + f01*(x-y0) + f012*(x-y0).*(x-y1) + ...
f0123*(x-y0).*(x-y1).* (x-y2);

Algorithm: Hermite Interpolation
Algorithm: Hermite Interpolation, Part #1

Given the data points (xi , f (xi ), f (xi )), i = 0, . . . , n.
Step 1: FOR i=0:n
y2i = xi , Q2i,0 = f (xi ), y2i+1 = xi , Q2i+1,0 = f (xi )
Q2i+1,1 = f (xi )
IF i > 0
Qi,0 Qi1,0
Q2i,1 =
y2i y2i1
END
END

Algorithm: Hermite Interpolation
Algorithm: Hermite Interpolation, Part #2

Step 2: FOR i = 2 : (2n + 1)
FOR j = 2 : i
Qi,ji Qi1,j1
Qi,j =
.
yi yij
END
END
Result: qi = Qi,i , i = 0, . . . , 2n + 1 now contains the coefficients for
2n+1
k1
X
Y
qi
H2n+1 (x) = q0 +
(x yj ) .
k=1
j=0

Higher Order Osculating Polynomials
1 of 3
So far we have seen the osculating polynomials of order 0 the

Lagrange polynomial, and of order 1 the Hermite interpolatory
polynomial.
It turns out that generating osculating polynomials of higher order
is fairly straight-forward; and we use Newtons divided
differences to generate those as well.

1 of 3
So far we have seen the osculating polynomials of order 0 the

Lagrange polynomial, and of order 1 the Hermite interpolatory
polynomial.
It turns out that generating osculating polynomials of higher order
is fairly straight-forward; and we use Newtons divided
differences to generate those as well.
n,k
Given a set of points {xk }nk=0 , and {f () (xk )}k=0,=0
; i.e. the
function values, as well as the first k derivatives of f in xk . (Note
that we can specify a different number of derivatives in each point.)
Set up the Newton-divided-difference table, and put in (k + 1)

duplicate entries of each point xk , as well as its function value
f (xk ).

2 of 3
Run the computation of Newtons divided differences as usual;

with the following exception:
Whenever a zero-denominator is encountered i.e. the divided
difference for that entry cannot be computed due to duplication of a point use a derivative instead. For mth divided
1 (m)
differences, use m!
f (xk ).
On the next slide we see the setup for two point in which two
derivatives are prescribed.


y
y0 = x0
f(x)
f [y0 ]
y1 = x0
f [y1 ]
y2 = x0
f [y2 ]
y3 = x1
f [y3 ]
y4 = x1
f [y4 ]
y5 = x1
f [y5 ]
1st Div. Diff.
3 of 3
2nd Div. Diff.
f [y0 , y1 ] =
f (x
0)
f [y1 , y2 ] =
f (x
0)
f [y2 , y3 ]
f [y3 , y4 ] =
f (x
1)
f [y4 , y5 ] =
f (x
1)
f [y0 , y1 , y2 ] = 12 f (x0 )
f [y1 , y2 , y3 ]
f [y2 , y3 , y4 ]
f [y3 , y4 , y5 ] = 12 f (x1 )
3rd Div. Diff.
f [y0 , y1 , y2 , y3 ]
f [y1 , y2 , y3 , y4 ]
f [y2 , y3 , y4 , y5 ]
3rd and higher order divided differences are computed as usual

in this case.
On the next slide we see four examples of 2nd order osculating
polynomials.

Examples...
1
0.8
0.6
f(0) = +0.00
f(1) = 1.00
df(0) = +1.00
df(1) = 1.00
ddf(0) = +0.00
ddf(1) = +0.00
0.4
0.2
0
0.2
0.4
f(1) = +1.00
0.6
df(1) = +1.00
ddf(1) = +0.00
0.8
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8

Examples...
1
0.8
0.6
f(0) = +0.00
f(1) = 1.00
df(0) = 1.00
df(1) = 1.00
ddf(0) = +0.00
ddf(1) = +0.00
0.4
0.2
0
0.2
0.4
f(1) = +1.00
0.6
df(1) = +1.00
ddf(1) = +0.00
0.8
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8

Examples...
0.5
f(0) = +0.00
f(1) = 1.00
df(0) = +0.00
df(1) = +1.00
ddf(0) = +10.00
ddf(1) = +3.00
f(1) = +1.00
0.5
df(1) = +2.00
ddf(1) = 1.00
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8

Examples...
1
0.8
0.6
f(0) = +0.00
f(1) = 1.00
df(0) = 1.00
df(1) = 1.00
ddf(0) = 10.00
ddf(1) = +10.00
0.4
0.2
0
0.2
0.4
f(1) = +1.00
0.6
df(1) = 1.00
ddf(1) = 10.00
0.8
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8

Examples...
1
0.8
0.6
f(0) = +0.00
f(1) = 1.00
df(0) = +1.00
df(1) = 1.00
ddf(0) = +0.00
ddf(1) = +0.00
0.8
0.6
0.4
0.4
0.2
0.2
0.2
0.2
0.4
f(0) = +0.00
f(1) = 1.00
df(0) = 1.00
df(1) = 1.00
ddf(0) = +0.00
ddf(1) = +0.00
0.4
f(1) = +1.00
0.6
f(1) = +1.00
0.6
df(1) = +1.00
ddf(1) = +0.00
0.8
df(1) = +1.00
ddf(1) = +0.00
0.8
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8
1
0.8
1
0.6
0.5
f(0) = +0.00
f(1) = 1.00
df(0) = +0.00
df(1) = +1.00
ddf(0) = +10.00
ddf(1) = +3.00
f(0) = +0.00
f(1) = 1.00
df(0) = 1.00
df(1) = 1.00
ddf(0) = 10.00
ddf(1) = +10.00
0.4
0.2
0
0.2
0.4
f(1) = +1.00
f(1) = +1.00
0.5
0.6
df(1) = 1.00
df(1) = +2.00
ddf(1) = 10.00
0.8
ddf(1) = 1.00
1
1
1
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8

We have encountered methods by these fellows
Sir Isaac Newton, 4 Jan 1643 (Woolsthorpe, Lincolnshire, England) 31

March 1727.
Joseph-Louis Lagrange, 25 Jan 1736 (Turin, Sardinia-Piedmont (now
Italy)) 10 April 1813.
Johann Carl Friedrich Gauss, 30 April 1777 (Brunswick, Duchy of
Brunswick (now Germany)) 23 Feb 1855.
Charles Hermite, 24 Dec 1822 (Dieuze, Lorraine, France) 14 Jan 1901.
The class website contains links to short bios for these (and other)
mathematicians, click on Mathematics Personae, who have contributed
to the material covered in this class. It makes for interesting reading and
puts mathematics into a historical and political context.

Divided Difference

Uploaded by

Copyright:

Available Formats

Divided Difference

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Divided Difference

Uploaded by

Copyright:

Available Formats

Polynomial Approximation: Practical Computations

Polynomial Approximation, Higher Order Matching

Numerical Analysis and Computing

San Diego State University

#5 Interpolation and Polynomial Approximation (1/41)

Polynomial Approximation: Practical Computations

Polynomial Approximation: Practical Computations

Polynomial Approximation, Higher Order Matching

Beyond Hermite Interpolatory Polynomials

Joe Mahaffy, hmahaffy@math.sdsu.edui

#5 Interpolation and Polynomial Approximation (2/41)

Polynomial Approximation: Practical Computations

Recap and Lookahead

#5 Interpolation and Polynomial Approximation (3/41)

Polynomial Approximation: Practical Computations

Recap and Lookahead

The coefficients of the polynomials. Once we have those, we can

For j = (n 1), (n 2), . . . , 1

#5 Interpolation and Polynomial Approximation (4/41)

Polynomial Approximation: Practical Computations

Joe Mahaffy, hmahaffy@math.sdsu.edui

#5 Interpolation and Polynomial Approximation (5/41)

Polynomial Approximation: Practical Computations

Note that we can evaluate this Horner-style, by writing

so that each step in the Horner-evaluation consists of a

#5 Interpolation and Polynomial Approximation (5/41)

Polynomial Approximation: Practical Computations

Finding the Constants {a0 , a1 , . . . , an }

Given the relation

This gets massively ugly fast! We need some nice clean

#5 Interpolation and Polynomial Approximation (6/41)

Polynomial Approximation: Practical Computations

Sir Isaac Newton to the Rescue: Divided Differences

Second Divided Difference:

f [xi+1 , xi+2 ] f [xi , xi+1 ]

f [xi+1 , xi+2 , . . . , xi+k ] f [xi , xi+1 , . . . , xi+k1 ]

Joe Mahaffy, hmahaffy@math.sdsu.edui

#5 Interpolation and Polynomial Approximation (7/41)

Polynomial Approximation: Practical Computations

The Constants {a0 , a1 , . . . , an } Revisited

We may suspect that a2 = f [x0 , x1 , x2 ], that is indeed so (a little

#5 Interpolation and Polynomial Approximation (8/41)

Polynomial Approximation: Practical Computations

Algebra: Chasing down a2 = f [x0 , x1 , x2 ]

(f (x2 ) f (x0 ))(x1 x0 ) (f (x1 ) f (x0 ))(x2 x0 )

(x1 x0 )f (x2 ) (x2 x0 )f (x1 ) + (x2 x0 x1 + x0 )f (x0 )

(x1 x0 )f (x2 ) (x1 x0 + x2 x1 )f (x1 ) + (x2 x1 )f (x0 )

(x1 x0 )(f (x2 ) f (x1 )) (x2 x1 )(f (x1 ) f (x0 ))

Joe Mahaffy, hmahaffy@math.sdsu.edui

#5 Interpolation and Polynomial Approximation (9/41)

Polynomial Approximation: Practical Computations

Newtons Interpolatory Divided Difference Formula

#5 Interpolation and Polynomial Approximation (10/41)

Polynomial Approximation: Practical Computations

Computing the Divided Differences (by table)

1st Div. Diff.

2nd Div. Diff.

f [x1 ,x2 ]f [x0 ,x1 ]

f [x2 ,x3 ]f [x1 ,x2 ]