Atkinson K. An Introduction to Numerical Analysis
Подождите немного. Документ загружается.
TRIGONOMETRIC INTERPOLATION 177
If
ian
I +
ibn
I
=/=
0,
then this function Pn(t)
is
called a trigonometric polynomial
of
degree n.
It
can
be shown by using trigonometric addition formulas that an
equivalent formulation is
n
Pn(t) = a
0
+ L
aJcos(t)J
1
+
/3
1
[sin(t)J
1
j=l
(3.8.3)
thus partially explaining
our
use
of
the word polynomial for such a function.
The
polynomial Pn(t) has period
2'1T
or
integral fraction thereof.
To
study
interpolation problems with Pn(t)
as
a solution,
we
must impose
2n
+ 1 interpolating conditions, since Pn(t) contains 2n + 1 coefficients a
1
,
b
1
.
Because
of
the
periodicity
of
the function
f(t)
and
the polynomial Pn(t), we also
require the interpolation nodes to lie in the interval
0
::;:
t <
2'1T
(or equivalently,
-
'1T
::;:
t <
'1T
or
0 < t
::;:
2'1T
).
Thus we assume
the
existence
of
the interpolation
nodes
(3.8.4)
and
we require Pn(t) to be chosen to satisfy
i = 0,
1,
...
,
2n
(3.8.5)
It
is shown
later
that this problem has a unique solution.
This interpolation problem has an explicit solution,
comparable
to the Lagrange
formula
(3.1.6)
for polynomial interpolation; this is dealt with in Problem 41.
Rather
than
proceeding with such a development, we first convert (3.8.4)-(3.8.5)
to
an
equivalent problem involving polynomials
and
functions
of
a complex
variable.
This
new formulation is the more
natural
mathematical setting for
trigonometric polynomial interpolation.
Using
Euler's formula
e;
8
= cos (
0)
+ i · sin ( 0)
i=Ff
we
obtain
e;e
+ e-;e
cos(O)
=
--
2
--
sin (0) =
--
2
-i-
Using
these
in
(3.8.2), we obtain
The
coefficients are related by
n
Pn(t) = L cieifr
j=
-n
1 :5:j::;: n
(3.8.6)
(3.8.7)
(3.8.8)
Given
{ cj}, the coefficients {a
1
,
b
1
}
are easily
obtained
by solving these latter
178
INTERPOLATION THEORY
equations. Letting z = e;
1
,
we
can rewrite (3.8.8) as the complex function
n
Pn(z) = L c
1
zj
(3.8.9)
j=
-n
The function znPn(z)
is
a polynomial of
degree~
2n.
To
reformulate the interpolation problem (3.8.4)-(3.8.5), let z
1
=
e;\
j = 0,
...
, 2n. With the restriction in (3.8.4), the numbers zj are distinct points on
the unit circle
JzJ
= 1
in
the complex plane. The interpolation problem
is
j = 0,
...
,2n
(3.8.10)
To see that this
is
always uniquely solvable, note that it
is
equivalent to
j = 0,
...
,2n
with
Q(z)
= znPn(z). This
is
a polynomial interpolation problem, with 2n + 1
distinct node points
z
0
,
•••
, z
2
n;
Theorem
3.1
shows there
is
a unique solution.
Also, the Lagrange formula (3.1.6) generalizes to
Q(z),
and thence to
P(z).
There are a number of reasons, both theoretical and practical, for converting
to the complex variable form of trigonometric interpolation. The most important
in
our
view
is
that interpolation and approximation by trigonometric polynomials
are intimately connected
to
the subject of differentiable functions of a complex
variable, and much of the theory
is
better understood from this pi!rspective. We
do
not
develop this theory, but a complete treatment
is
given
in
Henrici (1986,
chap. 13) and Zygmund (1959, chap. 10)
..
Evenly spaced interpolation The case of interpolation that
is
of most interest
in
applications
is
to use evenly spaced nodes
tj.
More precisely, define
2'TT
11
= j-2-n_+_1
j = 0,
±1,.±2,
...
(3.8.11)
The points t
0
,
•••
, t
2
n satisfy (3.8.4), and the points zj =
e;\
j =
0,
...
,
2n,
are
evenly spaced points on the unit circle
Jz
I =
1.
Note also that the points z
1
repeat
as
j increases by 2n +
1.
We now develop an alternative to the Lagrange form for
p,(t)
when the nodes
{
tj}
satisfy (3.8.11).
We
begin with the following lemma.
Lemma
4 For all integers k,
L
eikr
1
=
2n
+ 1
2n
{
j=O
0
(3.8.12)
The condition
e;
1
k
= 1
is
equivalent to k being an integer multiple of
2n
+
1.
TRIGONOMETRIC
INTERPOLATION 179
Proof
Let z =
e;
1
'.
Then using (3.8.11), eikr, =
eih,
and the sum in (3.8.12)
becomes
2n
S =
:E
zi
j=O
If
z =
1,
then this sums to
2n
+
1.
If
z
=t=
1,
then the geometric series
formula (1.1.8) implies
z2n+1-
1
S=
----
z-
1
Using (3.8.11), z
1
n+l
= e
1
"'ki =
1;
thus, S =
0.
The interpolation conditions (3.8.10) can be written
as
n
:E
ckeikrj =
/(
tJ
j =
0,
1,
...
,2n
k=
-n
•
(3.8.13)
To
find the coefficients
ck,
we
use Lemma
4.
Multiply equation j by
e-il\
then
sum over
j,
restricting I to satisfy - n
~
I
~
n. This yields
2n
n
2n
:E
:E
ckei(k-t)rj
=
:E
e-itr,f(tJ
(3.8.14)
j=O
k=
-n
j=O
Reverse the order of summation, and then use Lemma 4 to obtain
2n
'0
I
L
ei(k-t)rj
= ! k
=t=
J=O
t
2n
+ 1 k = I
Using this in (3.8.14),
we
obtain
1=-n,
...
,n
(3.8.15)
The coefficients { c_n,
...
,
en}
are called the finite Fourier transform of the data
{/(t
0
),
.••
,f(t
2
n)}. They yield an explicit formula for the trigonometric inter-
polating polynomial
Pn(t) of (3:8.8). The formula (3.8.15)
is
related to the
Fourier coefficients of
/(t):
-oo<l<oo
(3.8.16)
If
the trapezoidal numerical integration rule
[see
Section 5.1]
is
applied to these
integrals, using
2n
+ 1 subdivisions of
[0,
2'11"
],
then (3.8.15)
is
the result, provided
/(t)
is
periodic on
[0,2'11"].
We
next discuss the convergence of Pn(t)
to/(t).
180
INTERPOLATION THEORY
Theorem 3.6
Let
f(t)
be a continuous, periodic function, and let
27T
be an
integer multiple
of
its period. Define
Pn(f) = Infimum [
Max
1/(t)-
q(t)l]
deg(q).Sn
O.sr.s2or
(3.8.17)
with
q(t)
a trigonometric polynomial. Then the interpolating
function
Pn(t) from (3.8.8) and (3.8.15) satisfies
Max
1/(t)-
Pn(t)l::;;
c[ln(n
+ 2)]pn(/)
0.SI.S2or
n
~
0 (3.8.18)
The
constant c is independent
off
and
n.
Proof See
Zygmund
(1959, chap. 10, p. 19), since the proof is fairly
com-
plicated. Ill
The
quantity
PnU) is called a minimax error (see Chapter 4), and it can
be
estimated in a variety of ways.
The
most
important
bound
on PnU)
is
probably
that
of
D. Jackson. Assume that
f(t)
is k times continuously differentiable
on
[0,
27T
],
k
~
0,
and
further assume
J<k>(
t)
satisfies the condition
for some
0 < a
::;;
1.
(This
is
called a HOlder condition.) Then
n
~
1
(3.8.19)
with
ck(f)
independent
of
n.
For
a proof, see Meinardus (1967, p. 55).
An
alternative
error formula to that of (3.8.18)
is
given in Henrici (1986, cor.
13.6c), using
the
Fourier series coefficients (3.8.16) for
f(t).
Example
Consider
approximating
/(t)
=
esin(ll,
using the interpolating func-
tion
Pn(t).
The
maximum error
En=
Max
1/(t)-
Pn(t)l
O.st.s2.,.
for various values of n,
is
given in Table 3.15.
The
convergence
is
rapid.
Table 3.15
Error
in
trigonometric
polynomial interpolation
n
En
n
En
1
5.39E-
1
6
6.46E-
6
2
9.31E-
2
7
4.01E-
7
3
l.lOE-
2
8
2.22E-
8
4
1.11E-
3
9
l.lOE-
9
5
9.11E-
5
10
5.00E-
11
TRIGONOMETRIC
INTERPOLATION
181
The
fast Fourier transform The approximation
off(
t)
by
Pn(
t)
in the preceding
example was very accurate for small values of
n. In contrast, the calculation of
the finite Fourier transform
(3.8.15) in other applications will often require large
values of
n.
We
introduce a method that
is
very useful in reducing the cost of
calculating the coefficients {
c
1
}
when n
is
large.
Rather than using formula
(3.8.15),
we
consider the equivalent formula
1
m-1
"
"k
dk
=-
£....
w..:r
·
Jj
w =
e-2wijm
m
k = 0,
1,
...
, m - 1 (3.8.20)
m
j=O
with given data { f
0
,
•••
, fm-
1
}.
This
is
called a finite Fourier transform of order
m.
For
formula (3.8.15), let m =
2n
+
1.
We
can
allow k to be any integer,
noting that
-:oo<k<oo
(3.8.21)
Thus
it is sufficient to compute d
0
,
•••
, dm_
1
or
any other m consecutive
coefficients
d
k·
To
contrast the formula (3.8.20) with the alternative presented below,
we
calculate the cost of evaluating d
0
,
•••
, dm_
1
using (3.8.20). To evaluate dk, let
zk =
w;.
Then
1
m-1
dk
=-
L Jjz£
m
j=O
(3.8.22)
Using nested multiplication, this requires m - 1 multiplications and m - 1
additions. We ignore the division by
m, since often other factors are used. The
evaluation
of
zk
requires only 1 multiplication, since zk = wmzk_
1
,
k
;::_.;
2.
The
total cost of evaluating
d
0
,
•••
, dm_
1
is m
2
multiplications and
m(m
- 1)
additions.
To
introduce the main idea behind the fast Fourier transform, let m =
pq
with
p and q positive integers greater than
1.
Rewrite the definition (3.8.20) in the
equivalent form
1 p
-1
1
q-1
d = - " - " wk(l+pg)r
k
£....
£....
m
li+n
p 1=0 q
g=O
Use
w,t;
=
exp(
-2'1Ti/q)
=
wq.
Then
1p-1
[1q-l
l
d
--
kl
-
kg
k - L
Wm
L
Wq
fl+pg
. p 1=0 q g=O
k = 0, 1,
...
,
m-
1
Write
1
q-1
e<l)
= - " wkgr
k q
£....
q
1/+pq
g=O
O:::;/:::;p-1
(3.8.23)
1 p-1
d = - "
wkleU>
k
£....
m k
p 1=0
O:::;k:::;m-1
(3.8.24)
182 INTERPOLATION THEORY
Once {
ei'>}
is
known, each value of
dk
will
require p multiplications, using a
nested multiplication scheme
as
in
(3.8.22). The evaluation of (3.8.24)
will
require
mp
multiplications, assuming all
ei'>
have been computed previously. There
will
be a comparable number of additions.
We turn our attention
to
the computation of the quantities
ei'>.
The index k
ranges from 0 to m -
1,
but not all of these need
be
computed. Note that
1
q-1
eU>
= - "
w<k+q)gr
=
eU>
k+q
q
L...
q
J/+pq
k
g=O
because
wqq
=
1.
Thus
ei'>
needs to be calculated for only k =
0,
1,
...
, q -
1,
and then it repeats itself. For each
/,
{
e6'>,
... ,
e~
1
2r}
is
the finite Fourier
transform of the data { /
1
,
f
1
+p•
.••
, f
1
+p(q-l)
},
0
~
I
~
p -
1.
Thus the computa-
tion of {
eil)} amounts to the computation of p finite Fourier transforms of order
q (i.e., for data
of
length q ).
The fast Fourier transform amounts
to
a repeated use of this idea, to reduce
the computation to finite Fourier transforms of smaller and smaller order. To be
more specific, suppose
m =
2r
for some integer
r.
As
a first step, let p =
2,
q =
2r-l.
Then the computation of (3.8.24) requires 2m multiplications plus the
cost of evaluating two finite Fourier transforms of order
q =
2r-l.
For each of
these, repeat the process recursively. There
will
be r levels in this process,
resulting eventually in the evaluation of finite Fourier transforms of order one.
The total number of multiplications
is
given
by
which sums to
2rm =
2m
·
log
2
m
(3.8.25)
Thus the number of operations
is
proportional
to
m · log
2
m, in contrast to the
m
2
operations of the nested multiplication algorithm of (3.8.22). When m
is
large, say 2
10
or
.larger, this results in an enormous savings in calculation time.
For the particular case of
m =
2r,
a more careful accounting will show that only
m · log
2
m multiplications are actually needed, and there
is
a variant procedure
that requires only half
of this.
For
other values of m, there are modifications of the preceding procedure that
will still yield
an
operations count proportional
to
m · log m, just as previously
shown,
but
the case of m =
2r
leads
to
the greatest savings.
For
a further
discussion
of
this topic,
we
refer the reader
to
Henrici (1986, chap. 13). Henrici
also contains a discussion of the stability of the algorithm when rounding errors
are taken into account. The use of the fast Fourier transform has revolutionized
many subjects, making computationally feasible many calculations that previ-
ously were
not
practical.
DISCUSSION
OF
THE LITERATURE 183
Discussion of the Literature
As noted in the introduction, interpolation theory
is
a foundation for the
development of methods in numerical integration and differentiation, approxima-
tion theory, and the numerical solution of differential equations. Each of
these·
topics
is
developed in the following chapters, and the associated literature
is
discussed at that point. Additional results on interpolation theory are given in
de Boor (1978), Davis (1963), Henrici
(1982,
chaps. 5 and
7),
and Hildebrand
(1956).
For
an historical account of many of the topics of this Ghapter, see
Goldstine (1977).
The introduction of digital computers produced a revolution
in
numerical
analysis, including interpolation theory. Before the
use
of digital computers, hand
computation
was
necessary, which meant that numerical methods were used that
minimized the need for computation.
Such methods were often more complicated
than the methods now
used
on
computers, taking special advantage of the unique
mathematical characteristics of each problem. These methods also made exten-
sive use of tables,
to
avoid repeating calculations done by others; interpolation
formulas based on finite differences were used extensively. A large subject
was
created, called the finite difference calculus, and it
was
used in solving problems
in several areas of numerical analysis and applied mathematics. For a general
introduction to this approach to numerical analysis, see Hildebrand (1956) and
the references contained therein.
The use of digital computers has changed the needs of other areas for
interpolation theory, vastly reducing the need for finite difference based interpo-
lation formulas. But there
is
still an important place for both hand computation
and the use of mathematical tables, especially for the more complicated functions
of
mathematical physics. Everyone doing numerical work should possess an
elementary book of tables such
as
the well-known CRC tables. The National
Bureau of Standards tables of Abramowitz and
Stegun (1964) are an excellent
reference for nonelementary functions. The availability of sophisticated hand
calculators and microcomputers makes possible a new
level
of hand (or personal)
calculation.
Piecewise polynomial approximation theory has been very popular since the
early
1960s, and it
is
finding
use
in
a variety of fields. For example, see Strang
and Fix (1973, chap.
1)
for applications
to
the solution of boundary value
problems for ordinary differential equations, and see Pavlidis (1982, chaps.
10-12)
for applications in computer graphics. Most of the interest
in
piecewise
polynomial functions has centered on spline functions. The beginning of the
theory of spline functions
is
generally credited to Schoenberg
in
his
1946
papers,
and he has been prominent
in
helping
to
develop the subject [e.g., see Schoenberg
(1973)]. There
is
now
an extensive literature on spline functions, involving many
individuals and groups. For general surveys, see Ahlberg et
al.
(1967), de Boor
(1978), and Schumaker (1981).
Some of the most widely used computer software
for using spline functions
is
based on the programs in de Boor (1978). Versions of
these are available in the
IMSL and NAG numerical analysis libraries.
Finite Fourier transforms, trigonometric interpolation, and associated topics
are quite old topics; for example,
see
Goldstine (1977,
p.
238) for a discussion of
184
INTER
PO
LA
TION
THEORY
Gauss'
work on trigonometric interpolation. Since
Fourier
series and Fourier
transforms are
important
tools in much of applied mathematics, it
is
not
surprising
that
there
is
a great deal of interest in their discrete approximations.
Following
the
famous paper
by
Cooley and Tukey (1965) on the fast Fourier
transform, there has been a large increase in the use
of
finite Fourier transforms
and
associated topics.
For
example, this has led to very fast methods for solving
Laplace's partial differential equation
on
rectangular regions, which
we
discuss
further in
Chapter
8.
For
a classical account of trigonometric interpolation, see
Zygmund (1959, chap.
10), and for a more recent survey
of
the entire area
of
finite Fourier analysis, see Henrici (1986, chap. 13).
Multivariate polynomial interpolation theory is a rapidly developing area,
which for reasons
of
space has been omitted in this text.
The
finite element
method
for solving partial differential equations makes extensive use of multi-
variate interpolation theory, and some
of
the better presentations
of
this theory
are contained in books on the finite element method.
For
example, see Jain
(1984), Lapidus
and
Pinder (1982), Mitchell and Wait (1977),
and
Strang and Fix
(1973). More recently, work in computer graphics has led to new developments
[see Barnhill (1977) and
Pavlidis (1982, chap.
13)}.
Bibliography
Abramowitz, M., and I. Steguri, eds. (1964). Handbook
of
Mathematical Func-
tions.
National
Bureau of Standards, Washington, D.C. (Available now
from Dover, New York.)
Ahlberg,
J.,
E.
Nilson,
and
J.
Walsh (1967). The Theory
of
Splines and Their
Applications. Academic Press, New York.
Barnhill, R. (1977). Representation
and
approximation
of
surfaces. In Mathe-
matical Software
III,
J.
Rice, ed., pp. 69-120. Academic Press, New York.
De
Boor,
C.
(1978). A Practical Guide
to
Splines. Springer-Verlag, New York.
De
Boor,
C.
(1985). Convergence of cubic spline interpolation with the not-a-knot
condition. Math. Res. Ctr. Tech. Rep. 1876, Madison, Wis.
CRC
Standard Mathematical Tables. Chern.
Rubber
Publ. Co., Cleveland. (Pub-
lished yearly).
Cooley, J.,
and
J. Tukey (1965).
An
algorithm for the machine calculation of
complex
Fourier
series. Math. Comput.,
19,
297-301.
Davis,
P.
(1963). Interpolation and Approximation.
Ginn
(Blaisdell), Boston.
Goldstine,
H.
(1977). A History
of
Numerical Analysis from the 16th through the
19th Century. Springer-Verlag, New York.
Henrici, P. (1982). Essentials
of
Numerical Analysis. Wiley, New York.
Henrici, P. (1986). Applied and Computational Complex Analysis, Vol.
3.
Wiley,
New
York.
PROBLEMS
185
Hildebrand, F. (1956). Introduction
to
Numerical Analysis. McGraw-Hill,
New
York.
Isaacson,
E.,
and H. Keller
(1966).
Analysis
of
Numerical Methods. Wiley, New
York.
Jain,
M.
(1984). Numerical Solution
of
Differential Equations, 2nd ed.,
Wiley,
New
Delhi.
Lapidus,
L., and G. Pinder
(1982).
Numerical Solution
of
Partial Differential
Equations
in
Science and Engineering. Wiley, New York.
Lorentz, G.,
K.
Jetter, and
S.
Riemenschneider (1983).
Birkhoff
Interpolation. In
Encylopedia
of
Mathematics and Its Applications, Vol.
19.
Addison-Wesley,
Reading, Mass.
Meinardus, G. (1967).
Approximation
of
Functions: Theory and Numerical Meth-
ods
(trans!.
L.
Schumaker). Springer-Verlag, New York.
Mitchell,
A.,
and
R.
Wait
(1977).
The Finite Element Method
in
Partial Differen-
tial Equations.
Wiley, London.
Pavlidis, T. (1982). Algorithms for Graphics and Image Processing. Computer
Science
Press, Rockville, Md.
Powell,
M.
(1981). Approximation Theory and Methods. Cambridge Univ. Press,
Cambridge, England.
Schoenberg,
I.
(1946). Contributions to the approximation of equidistant data
by
analytic functions. Quart. Appl. Math, 4(Part A), 45-99; (Part
B),
112-141.
Schoenberg,
I.
(1973). Cardinal Spline Interpolation." Society for Industrial and
Applied Mathematics,
Philadelphia.
Schumaker,
L.
(1981). Spline Functions: Basic Theory: Wiley, New York.
Stoer, J., and
R.
Bulirsch
(1980).
Introduction
to
Numerical Analysis. Springer-
Verlag, New
York.
Strang, G., and G. Fix
(1973).
An
Analysis
of
the Finite Element Method.
Prentice-Hall, Englewood
Cliffs,
N.J.
Zygmund,
A. (1959). Trigonometric Series,
Vols.
1 and
2.
Cambridge Univ. Press,
Cambridge.
Problems
I. Recall the Vandermonde matrix X given in (3.1.3),
and
define
1
Xo
x2
0
xn
0
Vn(x) = det
1
2
x:-r
xn-1
Xn-1
1 X
x2
xn
186 INTERPOLATION THEORY
(a) Show that Vn(x)
is
a polynomial of degree n, and that its roots are
Xo,
.•.
, xn-1· Obtain the formula
Hint:
Expand the last
row
of Vn(x) by minors
to
show that Vn(x)
is
a
polynomial of degree
n and
to
find
the coefficient of the term xn.
(b) Show
det
(X)
=
V,(xn)
=
TI
(x;-
x)
0:5j<i:5:.n
2.
For
the basis functions
lj,n(x)
given
in
(3.1.5), prove that for any n
~
1,
n
L
lj,n(x)
= 1
j=O
for all x
3. Recall the Lagrange functions 1
0
(x),
...
, ln(x), defined in (3.1.5) and then
rewritten in a slightly different form
in
(3.2.4), using
Let
wj
=
['lt;(xj)]-
1
•
Show that the polynomial Pn(x) interpolating
f(x)
can be written
as
n
L [
wJ(xj)
]!(x-
xj)
Pn
(X) =
_j_-_O---:n------
L
w/(x-
x)
j-0
provided x
is
not a node point. This
is
called the barycentric representation
of
Pn(x), giving it
as
a weighted sum of the values
{/(x
0
),
•••
,
f(xn)}.
For
a discussion of the use of this representation, see Henrici
(1982,
p.
237).
4. Consider linear interpolation in a table of values of
eX,
0
.::;;
x
.::;;
2,
with
h = .01. Let the table values be given
to
five
significant digits,
as
in the
CRC
tables. Bound the error of linear interpolation, including that part due
to the rounding errors
in
the table entries.
5. Consider linear interpolation in a table of cos (
x)
with x given in degrees,
0
.::;;
x
.::;;
90°, with a stepsize h = 1' =
'io
degree. Assuming that the table
entries are given to
five
significant digits, bound the total error of interpola-
tion.
6. Suppose you are to
I11ake
a table of values of sin ( x )
1
0
.::;;
x
.::;;
'1T
/2,
with a
stepsize
of
h.
Assume linear interpolation
is
to
be used with the table, and
suppose the total error, including the effects due
to
rounding in table
entries, is
to
be at most
10-
6
•
What should h equal (choose it in a