Interpolation and Approximation Theory

Finding a polynomial of at most degree n to pass through n + 1 points in the interval

[a, b] is referred to as interpolation.
Approximation theory deals with two types of problems.
Given a data set, one seeks a function best tted to this data set, for example,
given {(x
1
, y
1
), (x
2
, y
2
), , (x
n
, y
n
)}, one seeks a line y = mx + b which best ts
this data set.
Given an explicit function, one seeks a simpler function for representation, for
example, use 1 + x +
x
2
2!
+
x
3
3!
to represent e
x
.
Lagrange Polynomial Interpolation
Newtons Divided-Dierence Formula
Hermite Polynomial Interpolation
Cubic spline interpolation
Bezier curves
Cubic B-splines
Orthogonal functions
Trigonometric functions
Chebyshev polynomials
Legendre polynomials
Laguerre polynomials
Gamma functions
Beta functions
Bessel functions
Other Topics with Applications
Polynomial Approximation
Suppose that the function f(x) = e
x
is to be approximated by a polynomial of degree
2 over the interval [-1, 1]. The approximations by Taylor polynomial 1 + x + 0.5x
2
and
Chebyshev polynomial 1 + 1.17518x + 0.54309x
2
are given below.
1.5 1 0.5 0 0.5 1 1.5
0
0.5
1
1.5
2
2.5
3
Taylor Approximation for e
x
1.5 1 0.5 0 0.5 1 1.5
0
0.5
1
1.5
2
2.5
3
Chebyshev Approximation for e
x
Figure 1: Polynomial Approximations for e
x
over [-1, 1]
X=-1:0.1:1;
Y0=exp(X);
Y1=1.0000+1.0000*X+0.5000*X.^2; % Taylor Expansion
Y2=1.0000+1.17518*X+0.54309*X.^2; % Chebyshev Polynomial by Chaurchin
%Y2=1.0000+1.129772*X+0.532042*X.^2; % Chebyshev Polynomial from Textbook
V=[-1.5 1.5, 0 3];
subplot(2,1,1)
plot(X,Y0,b--,X,Y1,r-); axis(V); grid;
title(Taylor Approximation for e^x)
subplot(2,1,2)
plot(X,Y0,b--,X,Y2,r-); axis(V); grid;
title(Chebyshev Approximation for e^x)
Taylor Polynomial Approximation
Suppose that f C
n+1
[a, b] and x
0
[a, b] is a xed value. If x [a, b], then
f(x) = P
n
(x) + E
n
(x)
where P
n
(x) is a polynomial that can be used to approximate f(x) by
f(x) P
n
(x) =
n

k=0
f
(k)
(x
0
)
k!
(x x
0
)
k
having some c between x and x
0
such that
E
n
(x) =
f
(n+1)
(c)
(n + 1)!
(x x
0
)
n+1
|e P
15
(1)| = |e 2.718282818459| <
e
16!
<
3
16!
< 1.433844 10
13
|sin(x) P
9
(x)| <
1
10!
2.75574 10
7
for |x| 1, where
P
9
(x) = x
x
3
3!
+
x
5
5!

x
7
7!
+
x
9
9!
|cos(x) P
8
(x)| <
1
9!
2.75574 10
6
for |x| 1, where
P
8
(x) = 1
x
2
2!
+
x
4
4!

x
6
6!
+
x
8
8!
Polynomial Interpolation
We attempt to nd a polynomial of at most degree n to pass through n+1 points in the
interval [a, b].
[x
0
, y
0
]
t
, [x
1
, y
1
]
t
, , [x
n
, y
n
]
t
where
a x
0
< x
1
< < x
n
b
0 1 2 3 4 5 6
0
1
2
3
4
5
6
7
y=[5x
4
82x
3
+427x
2
806x+504]/24
Figure 2: Polynomial Passing Through Five Points
%
% Script File: func4.m
% A quadric function for interpolation: y=f(x)=[5x^4-82x^3+427x^2-806x+504]/24
%
X=0.6:0.1:5.2;
Y=(5*X.^4-82*X.^3+427*X.^2-806*X+504)/24.0;
V=[0 6, 0 7];
plot(X,Y,b-,[1 2 3 4 5],[2 1 5 6 1],ro); axis(V); grid
title(y=[5x^4-82x^3+427x^2-806x+504]/24)
Polynomials for Interpolation
Theorem: Suppose that the function y = f(x) is known at the n + 1 distinct points
[x
0
, y
0
]
t
, [x
1
, y
1
]
t
, , [x
n
, y
n
]
t
where
a x
0
< x
1
< < x
n
b
Then there is a unique polynomial P
n
(x) of degree at most n such that
P
n
(x
i
) = y
i
0 i n
If the error function E(x) = f(x) P
n
(x) is required, then we need to know f
(n+1)
(x)
whose bound of magnitude is
max{|f
(n+1)
(x)| : a x b}
A Lagrange polynomial of degree n
L
n,k
(x) =

n
j=k
(x x
j
)

n
j=k
(x
k
x
j
)
Error Formula for Lagrange Polynomial
f(x) =
n

k=0
f(x
k
)L
n,k
(x) +
f
(n+1)
(
x
)
(n + 1)!
n

k=0
(x x
k
)
for some unknown number
x
that lies in the smallest interval that contains x
0
, x
1
, , x
n
,
and x.
Polynomials in Newton Form
P
n
(x) = P
n1
(x) + a
n
n1

j=0
(x x
j
)
Polynomials in Chebyshev Form
P
n
(x) =
0
+
1
T
1
(x) +
2
T
2
(x) + +
n
T
n
(x)
where
T
n
(x) = cos(ncos
1
x), T
0
(x) 1, T
1
(x) = x, T
2
(x) = 2x
2
1, T
3
(x) = 4x
3
3x.
Hermite Polynomials H
n
(x)
An Example for Polynomial Interpolation
We look for polynomials of degree at most 3 to interpolate the following four points.
x 5 -7 -6 0
y 1 -23 -54 -954
Table 1: P
3
(x) = 4x
3
+ 35x
2
84x 954
.
Solution in Lagrange form
P
3
(x) = 1
(x+7)(x+6)(x0)
(5+7)(5+6)(50)
+ (23)
(x5)(x+6)(x0)
(75)(7+6)(70)
+ (54)
(x5)(x+7)(x0)
(65)(6+7)(60)
+ (954)
(x5)(x+7)(x6)
(05)(0+7)(0+6)
Solution in Newton form
P
3
(x) = 1 + 2(x 5) + 3(x 5)(x + 7) + 4(x 5)(x + 7)(x + 6)
Solution in Chebyshev form
P
3
(x) = 936.5 81T
1
(x) + 17.5T
2
(x) + T
3
(x)
where
T
n
(x) = cos(ncos
1
x), T
0
(x) 1, T
1
(x) = x, T
2
(x) = 2x
2
1, T
3
(x) = 4x
3
3x.
Divided Dierences
Suppose that the function y = f(x) is known at the n + 1 points
[x
0
, f(x
0
)]
t
, [x
1
, f(x
1
)]
t
, , [x
n
, f(x
n
)]
t
, where a x
0
< x
1
< < x
n
b
The n + 1 zeroth divided dierences of f are dened as
f[x
i
] = f(x
i
) 0 i n
The rst divided dierences of f are dened as
f[x
i
, x
i+1
] =
f[x
i+1
] f[x
i
]
x
i+1
x
i
0 i n 1
The kth divided dierences can be inductively dened by
f[x
i
, x
i+1
, , x
i+k1
, x
i+k
] =
f[x
i+1
, x
i+2
, , x
i+k
] f[x
i
, x
i+1
, , x
i+k1
]
x
i+k
x
i
0 i nk
The nth divided dierence is
f[x
0
, x
1
, , x
n
] =
f[x
1
, x
2
, , x
n
] f[x
0
, x
1
, , x
n1
]
x
n
x
0
It can be shown that the nth Lagrange interpolation polynomial w.r.t. x
0
< x
1
< < x
n
can be expressed as Newton (interpolatory) divided-dierence formula
P
n
(x) = f[x
0
] + f[x
0
, x
1
](x x
0
) + f[x
0
, x
1
, , x
n
](x x
0
)(x x
1
) (x x
n1
)
= f[x
0
] +

n
k=1
f[x
0
, x
1
, , x
k
](x x
0
)(x x
1
) (x x
k1
)
(1)
Newton (interpolatory) divided-dierence formula has simpler form when x
j
x
j1
=
h 1 j n. Let x = x
0
+ sh, then x x
i
= (s i)h, then the formula ?? becomes
P
n
(x) = P
n
(x
0
+ sh) = f[x
0
] +

n
k=1
s(s 1) (s k + 1)h
k
f[x
0
, x 1, , x
k
]
= f[x
0
] +

n
k=1
_
s
k
_
k!h
k
f[x
0
, x
1
, , x
k
]
= f[x
0
] +

n
k=1
_
s
k
_

k
f(x
0
)
= f[x
n
] +

n
k=1
(1)
k
_
s
k
_

k
f(x
n
)
Hermite Interpolation and Polynomial
If f C
1
[a, b] and a x
0
< x
1
< < x
n
b, the unique polynomial of least degree
which agrees with f and f

at x
0
, x
1
, , x
n
is the polynomial of degree at most 2n + 1
given by
H
2n+1
(x) =
n

j=0
f(x
j
)H
n,j
+
n

j=0
f

(x
j
)

H
n,j
(x)
where
H
n,j
(x) = [1 2(x x
j
)L

n,j
(x
j
)]L
2
n,j
(x)

H
n,j
(x) = (x x
j
)L
2
n,j
(x
j
)
L
n,j
(x) =
(x x
0
)(x x
1
) (x x
j1
)(x x
j+1
) (x x
n
)
(x
j
x
0
)(x
j
x
1
) (x
j
x
j1
)(x
j
x
j+1
) (x
j
x
n
)
Show that H
2n+1
(x
k
) = f(x
k
) and H

2n+1
(x
k
) = f

(x
k
) k = 0, 1, , n.
Error Formula
If f C
2n+2
[a, b], then
f(x) = H
2n+1
(x) +
f
(2n+2)
(
x
)
(2n + 2)!
(x x
0
)
2
(x x
1
)
2
(x x
n
)
2
for some
x
(a, b).
Cubic Spline Interpolation
Given a function f dened on [a, b] and a set of n+1 nodes a = x
0
< x
1
< < x
n
=
b, a cubic spline interpolant, S, for f is a function that satises the following conditions:
(1) For each j = 0, 1, , n 1, S(x) is a cubic polynomial, denoted by S
j
(x), on the
subinterval [x
j
, x
j+1
).
(2) S(x
j
) = f(x
j
) for each j = 0, 1, , n.
(3) S
j+1
(x
j+1
) = S
j
(x
j+1
) for each j = 0, 1, , n 2.
(4) S

j+1
(x
j+1
) = S

j
(x
j+1
) for each j = 0, 1, , n 2.
(5) S

j+1
(x
j+1
) = S

j
(x
j+1
) for each j = 0, 1, , n 2.
(6) One of the following sets of boundary conditions is satised:
(a) S

(x
0
) = S

(x
n
) = 0 (natural or free boundary);
(b) S

(x
0
) = f

(x
0
) and S

(x
n
) = f

(x
n
) (clamped boundary).
x 0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0
f(x) 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25
x 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3
f(x) 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25
Table 2: A ruddy duck in ight
Finding A Cubic Spline Interpolant
Let S
j
(x) = a
j
+b
j
(xx
j
)+c
j
(xx
j
)
2
+d
j
(xx
j
)
3
, h
j
= x
j+1
x
j
, for 0 j n1,
From (2), a
j
= S
j
(x
j
) = f(x
j
), 0 j n 1, and denote a
n
= f(x
n
).
From (3), a
j+1
= a
j
+ b
j
(x
j+1
x
j
) + c
j
(x
j+1
x
j
)
2
+ d
j
(x
j+1
x
j
)
3
, 0 j n 2.
(A) a
j+1
= a
j
+ b
j
h
j
+ c
j
h
2
j
+ d
j
h
3
j
, 0 j n 1, where a
n
= f(x
n
).
Similarly, S

j
(x) = b
j
+ 2c
j
(x x
j
) + 3d
j
(x x
j
)
2
, 0 j n 1.
(B) b
j+1
= b
j
+ 2c
j
h
j
+ 3d
j
h
2
j
, 0 j n 1 by (4).
Dene c
n
=
1
2
S

(x
n
), and by using (5), we have
(C) c
j+1
= c
j
+3d
j
h
j
, 0 j n1, and c
n1
+3d
n1
h
n1
= c
n
= 0 by using (6)(a).
(C

) d
j
=
1
3h
j
(c
j+1
c
j
), 0 j n 1, substitute (C

) into (A) and (B), we have

(D) a
j+1
= a
j
+ b
j
h
j
+
h
2
j
3
(2c
j
+ c
j+1
), 0 j n 1
(E) b
j+1
= b
j
+ h
j
(c
j
+ c
j+1
), 0 j n 1, or
(E

) b
j
= b
j1
+ h
j1
(c
j1
+ c
j
), 1 j n
From (D), we have
(F) b
j
=
1
h
j
(a
j+1
a
j
)
h
j
3
(2c
j
+ c
j+1
), 0 j n 1, or
(F

) b
j1
=
1
h
j1
(a
j
a
j1
)
h
j1
3
(2c
j1
+ c
j
), 1 j n.
Substitute (F) and (F

) into (E

), we have
(G) h
j1
c
j1
+2(h
j1
+h
j
)c
j
+h
j
c
j+1
=
3
h
j
(a
j+1
a
j
)
3
h
j1
(a
j
a
j1
), for 1 j n1.
Thus the problem is reduced to solving Ac = h with (n1) equations and (n1) unknown
variables c = [c
1
, c
2
, , c
n1
]
t
by using the boundary conditions c
0
=
1
2
S

(x
0
) = 0 and
c
n
=
1
2
S

(x
n
) = 0.
Once {c
j
, 0 j n 1} are solved, {d
j
, 0 j n1} and {b
j
, 0 j n1} could
be easily solved by using (C

) and (F

), respectively.
_

_
2(h
0
+ h
1
) h
1
0 0 0
h
1
2(h
1
+ h
2
) h
2
0 0
0 h
2
.
.
. h
3

.
.
.
.
.
. 0
.
.
.
.
.
. 0
.
.
.
.
.
. 0 h
n3
.
.
. h
n2
0 0 h
n2
2(h
n2
+ h
n1
)
_

_
c
1
c
2
c
3
.
.
.
.
.
.
c
n1
_

_
= h
where
h =
_

_
3
h
1
(a
2
a
1
)
3
h
0
(a
1
a
0
)
3
h
2
(a
3
a
2
)
3
h
1
(a
2
a
1
)
.
.
.
.
.
.
.
.
.
3
h
n1
(a
n
a
n1
)
3
h
n2
(a
n1
a
n2
)
_

_
Cubic Spline Interpolant for A Ruddy Duck
%
% Script File: cspline.m
% Cubic Spline Interpolation for a rubby duck of 21 points
%
n=21;
fin=fopen(duck.txt);
fgetL(fin);
X=fscanf(fin,%f,n);
Y=fscanf(fin,%f,n);
X0=0.9:0.4:13.3;
Y0=spline(X,Y,X0);
plot(X,Y,b--o,X0,Y0,r-); axis([0.5 13.5, -1, 5]); grid
legend(Sample Points of A Duck,Cubic Spline Interpolant);
title(Cubic spline interpolant for a ruddy duck)
2 4 6 8 10 12
1
0
1
2
3
4
5
Cubic spline interpolant for a ruddy duck
Sample Points of A Duck
Cubic Spline Interpolant

Figure 3: Cubic Spline Interpolant for A Ruddy Duck

Bezier Curves and B-splines
Bezier curves and B-splines are widely used in computer graphics and computer-aided
design. These curves have good geometric property in that in changing one of the points
we change only one portion of the tted curve, a local eect. For cubic splines, changing
only one point might have a global eect.
Bezier curves are named after the French engineer, Pierre Bezier of the Renault
Automobile Company. He developed them in the early 1960s to ll a need for curves
whose shape can be practically controlled by changing a few parameters.
The nth degree Bezier polynomial determined by n + 1 points is given by
P(u) =
n

i=0
C
n
i
(1 u)
ni
u
i
P
i
Bezier cubics are commonly used. For 0 u 1, denote
x(u) = (1 u)
3
x
0
+ 3(1 u)
2
ux
1
+ 3(1 u)u
2
x
2
+ u
3
x
3
y(u) = (1 u)
3
y
0
+ 3(1 u)
2
uy
1
+ 3(1 u)u
2
y
2
+ u
3
y
3
Then
dx
du
= 3(x
1
x
0
),
dy
du
= 3(y
1
y
0
) at u = 0.
dy
dx
=
y
1
y
0
x
1
x
0
at P
0
,
dy
dx
=
y
2
y
3
x
2
x
3
at P
3
An Algorithm for drawing a Bezier curve
for i = 0, 3n 1, 3
for u = 0, 1, u
x(u) = (1 u)
3
x
i
+ 3(1 u)
2
ux
i+1
+ 3(1 u)u
2
x
i+2
+ u
3
x
i+3
y(u) = (1 u)
3
y
i
+ 3(1 u)
2
uy
i+1
+ 3(1 u)u
2
y
i+2
+ u
3
y
i+3
plot(x(u), y(u))
endfor
endfor
B-splines
The B-splines (basis of splines) are like Bezier curves in that they do not ordinarily
pass through the given data points. They can be of any degree, but cubic B-splines are
commonly used.
Given the points P
i
(x
i
, y
i
), i = 0, 1, , n, a portion of a cubic B-spline for the
interval (P
i
, P
i+1
), i = 1, 2, , n 1, is computed by
B
i
(u) =
2

k=1
b
k
P
i+k
where
b
1
=
(1 u)
3
6
, b
0
=
u
3
2
u
2
+
2
3
, b
1
=
u
3
2
+
u
2
2
+
u
2
+
1
6
, b
2
=
u
3
6
u-cubics act as weighting factors on the coordinates of the four successive points to
generate the curve, for example, at u = 0, the weights are [
1
6
,
2
3
,
1
6
, 0]; at u = 1, the
weights are [0,
1
6
,
2
3
,
1
6
].
An Algorithm for drawing a cubic B-spline
for i = 1, n 2
for u = 0, 1, u
x = x
i
(u)
y = y
i
(u)
plot(x,y)
endfor
endfor
where
x
i
(u) =
(1 u)
3
6
x
i1
+ [
u
3
2
u
2
+
2
3
]x
i
+ [
u
3
2
+
u
2
2
+
u
2
+
1
6
]x
i+1
+
u
3
6
x
i+2
y
i
(u) =
(1 u)
3
6
y
i1
+ [
u
3
2
u
2
+
2
3
]y
i
+ [
u
3
2
+
u
2
2
+
u
2
+
1
6
]y
i+1
+
u
3
6
y
i+2
Note that a B-spline does not necessarily pass through any point of P

i
s.
Approximation Theory
Approximation theory deals with two types of problems.
Given a data set, one seeks a function best tted to this data set, for example,
given {(x
1
, y
1
), (x
2
, y
2
), , (x
n
, y
n
)}, one seeks a line y = mx + b which best ts
this data set.
Given an explicit function, one seeks a simpler function for representation, for
example, use 1 + x +
x
2
2!
+
x
3
3!
to represent e
x
.
Orthogonal Functions
The set of functions {
0
,
1
, ,
n
} is said to be orthogonal for the interval [a, b]
with respect to the weight function w if
_
b
a

i
(x)
k
(x)w(x)dx =
_

k
> 0 if i = k
0 if i = k
(2)
{
0
,
1
, ,
n
} is said to be orthonormal if, in addition,
k
= 1 for 0 k n.
{1, cos x, sin x, , cos kx, sin kx, } with respect to w(x) 1 is orthogonal for the
interval [0, 2].
{
1

2
,
1

cos x,
1

sin x, ,
1

cos kx,
1

sin kx, } with respect to w(x) 1 is or-

thonormal for the interval [0, 2].
The set of Chebyshev polynomials {cos(ncos
1
x)}

n=0
is orthogonal with respect to
w(x) =
1

1x
2
for the interval [1, 1].
The set of Chebyshev polynomials {
1

[cos(ncos
1
x)]

n=1
} is orthonormal with
respect to w(x) =
1

1x
2
for the interval [1, 1].
The set of Legendre polynomials {P
n
(x) =
1
2
n
n!

d
n
(x
2
1)
n
dx
n
} is orthogonal with respect
to w(x) 1 for the interval [1, 1]. Note that
_
1
1
P
m
(x)P
n
(x)dx =
_

_
2
2n+1
for m = n
0 for m = n
(3)
Any high-order Legendre polynomial may be derived using the recursion formula
P
n
(x) =
2n 1
n
xP
n1
(x) +
n 1
n
P
n2
(x) (4)
Note that
P
0
(x) = 1, P
1
(x) = x, P
2
(x) =
1
2
(3x
2
1), P
3
(x) =
1
2
(5x
3
3x)