<< . .

. 49
( : 95)



. . >>

340 8. Polynomial Interpolation

De¬ne e = f ’ Π1 f . Since e(xj ) = 0 for all j = 0, . . . , K, Rolle™s theorem
h
infers the existence of ξj ∈ (xj , xj+1 ), for j = 0, . . . , K ’ 1 such that e (ξj ) = 0.
Since Π1 f is a linear function on each Ij , for x ∈ Ij we obtain
h
x x
e (x) = e (s)ds = f (s)ds,
ξj ξj

whence
xj+1
|e (x)| ¤ |f (s)|ds, for x ∈ [xj , xj+1 ]. (8.28)
xj

We recall the Cauchy-Schwarz inequality
1/2 1/2
β β β
u(x)v(x)dx ¤ 2 2
u (x)dx v (x)dx (8.29)
± ± ±

which holds if u, v ∈ L2 (±, β). If we apply this inequality to (8.28) we obtain
« 1/2 « 1/2
xj+1 xj+1
¬ · ¬ ·
|e (x)| ¤ |f (s)|2 ds
12 dx 
xj xj
« 1/2 (8.30)
xj+1
¬ ·
¤ h1/2  |f (s)|2 ds .
xj


To ¬nd a bound for |e(x)|, we notice that
x
e(x) = e (s)ds,
xj

so that, applying (8.30), we get
1/2
xj+1 xj+1
|e(x)| ¤ |e (s)|ds ¤ h |f (s)| ds
3/2 2
. (8.31)
xj xj

Then
xj+1 xj+1 xj+1 xj+1

|e (x)|2 dx ¤ h2 |f (s)|2 ds |e(x)|2 dx ¤ h4 |f (s)|2 ds,
and
xj xj xj xj

from which, summing over the index j from 0 to K ’ 1 and taking the square
root of both sides, we obtain
1/2 1/2
b b
|e (x)| dx ¤h |f (x)| dx
2 2
,
a a

and
1/2 1/2
b b
|e(x)| dx ¤h |f (x)| dx
2 2 2
,
a a
3
which is the desired estimate (8.27), with C1 = C2 = 1.
8.4 Hermite-Birko¬ Interpolation 341

8.4 Hermite-Birko¬ Interpolation
Lagrange polynomial interpolation can be generalized to the case in which
also the values of the derivatives of a function f are available at some (or
all) of the nodes xi .
Let us then suppose that (xi , f (k) (xi )) are given data, with i = 0, . . . , n,
n
k = 0, . . . , mi and mi ∈ N. Letting N = i=0 (mi + 1), it can be proved
(see [Dav63]) that, if the nodes {xi } are distinct, there exists a unique
polynomial HN ’1 ∈ PN ’1 , called the Hermite interpolation polynomial,
such that
(k) (k)
HN ’1 (xi ) = yi , i = 0, . . . , n k = 0, . . . , mi ,

of the form
n mi
(k)
HN ’1 (x) = yi Lik (x) (8.32)
i=0 k=0

(k)
where yi = f (k) (xi ), i = 0, . . . , n, k = 0, . . . , mi .
The functions Lik ∈ PN ’1 are called the Hermite characteristic polynomials
and are de¬ned through the relations

1 if i = j and k = p,
dp
(Lik )(xj ) =
dxp 0 otherwise.

De¬ning the polynomials
n mk +1
(x ’ xi )j x ’ xk
lij (x) = , i = 0, . . . , n, j = 0, . . . , mi ,
xi ’ xk
j! k=0
k=i



and letting Limi (x) = limi (x) for i = 0, . . . , n, we have the following recur-
sive formula for the polynomials Lij
mi
(k)
Lij (x) = lij (x) ’ j = mi ’ 1, mi ’ 2, . . . , 0.
lij (xi )Lik (x)
k=j+1

As for the interpolation error, the following estimate holds

f (N ) (ξ)
f (x) ’ HN ’1 (x) = „¦N (x) ∀x ∈ R
N!
where ξ ∈ I(x; x0 , . . . , xn ) and „¦N is the polynomial of degree N de¬ned
by

„¦N (x) = (x ’ x0 )m0 +1 (x ’ x1 )m1 +1 . . . (x ’ xn )mn +1 . (8.33)
342 8. Polynomial Interpolation

Example 8.6 (osculatory interpolation) Let us set mi = 1 for i = 0, . . . , n.
In this case N = 2n + 2 and the interpolating Hermite polynomial is called the
osculating polynomial, and it is given by
n
(1)
HN ’1 (x) = yi Ai (x) + yi Bi (x)
i=0


where Ai (x) = (1 ’ 2(x ’ xi )li (xi ))li (x)2 and Bi (x) = (x ’ xi )li (x)2 , for i =
0, . . . , n, with
n
1
li (xi ) = , i = 0, . . . , n.
xi ’ xk
k=0,k=i

As a comparison, we use Programs 65 and 67 to compute the Lagrange and
Hermite interpolating polynomials of the function f (x) = sin(4πx) on the interval
[0, 1] taking four equally spaced nodes (n = 3). Figure 8.4 shows the superposed
graphs of the function f (dashed line) and of the two polynomials Πn f (dotted

line) and HN ’1 (solid line).


1.5



1



0.5



0



’0.5



’1



’1.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1



FIGURE 8.4. Lagrange and Hermite interpolation for the function
f (x) = sin(4πx) on the interval [0, 1]



Program 67 computes the values of the osculating polynomial at the ab-
scissae contained in the vector z. The input vectors x, y and dy contain the
interpolation nodes and the corresponding function evaluations of f and
f , respectively.

Program 67 - hermpol : Osculating polynomial
function [herm] = hermite(x,y,dy,z)
n = max(size(x)); m = max(size(z)); herm = [];
for j = 1:m
xx = z(j); hxv = 0;
for i = 1:n,
den = 1; num = 1; xn = x(i); derLi = 0;
for k = 1:n,
8.5 Extension to the Two-Dimensional Case 343

if k ˜= i, num = num*(xx-x(k)); arg = xn-x(k);
den = den*arg; derLi = derLi+1/arg;
end
end
Lix2 = (num/den)ˆ2; p = (1-2*(xx-xn)*derLi)*Lix2;
q = (xx-xn)*Lix2; hxv = hxv+(y(i)*p+dy(i)*q);
end
herm = [herm, hxv];
end




8.5 Extension to the Two-Dimensional Case
In this section we brie¬‚y address the extension of the previous concepts to
the two-dimensional case, referring to [SL89], [CHQZ88], [QV94] for more
details. We denote by „¦ a bounded domain in R2 and by x = (x, y) the
coordinate vector of a point in „¦.


8.5.1 Polynomial Interpolation
A particularly simple situation occurs when „¦ = [a, b] — [c, d], i.e., the
interpolation domain „¦ is the tensor product of two intervals. In such a
case, introducing the nodes a = x0 < x1 < . . . < xn = b and c = y0 <
y1 < . . . < ym = d, the interpolating polynomial Πn,m f can be written as
n m
Πn,m f (x, y) = i=0 j=0 ±ij li (x)lj (y), where li ∈ Pn , i = 0, . . . , n, and
lj ∈ Pm , j = 0, . . . , m, are the characteristic one-dimensional Lagrange
polynomials with respect to the x and y variables respectively, and where
±ij = f (xi , yj ).
The drawbacks of one-dimensional Lagrange interpolation are inherited
by the two-dimensional case, as con¬rmed by the example in Figure 8.5.

Remark 8.1 (The general case) If „¦ is not a rectangular domain or if
the interpolation nodes are not uniformly distributed over a Cartesian grid,
the interpolation problem is di¬cult to solve, and, generally speaking, it is
preferable to resort to a least-squares solution (see Section 10.7). We also
point out that in d dimensions (with d ≥ 2) the problem of ¬nding an
interpolating polynomial of degree n with respect to each space variable on
n + 1 distinct nodes might be ill-posed.
Consider, for example, a polynomial of degree 1 with respect to x and y
of the form p(x, y) = a3 xy +a2 x+a1 y +a0 to interpolate a function f at the
nodes (’1, 0), (0, ’1), (1, 0) and (0, 1). Although the nodes are distinct, the
problem (which is nonlinear) does not in general admit a unique solution;
actually, imposing the interpolation constraints, we end up with a system
that is satis¬ed by any value of the coe¬cient a3 .
344 8. Polynomial Interpolation




8
0.5

0.4 6

0.3
4
0.2
2
0.1
0
0

’2
’0.1
5
5
5
5
0
0
0
0

’5
’5 ’5
’5




FIGURE 8.5. Runge™s counterexample extended to the two-dimensional case:
interpolating polynomial on a 6 — 6 nodes grid (left) and on a 11 — 11 nodes grid
(right). Notice the change in the vertical scale between the two plots

8.5.2 Piecewise Polynomial Interpolation
In the multidimensional case, the higher ¬‚exibility of piecewise interpola-
tion allows for easy handling of domains of complex shape. Let us suppose
that „¦ is a polygon in R2 . Then, „¦ can be partitioned into K nonover-
lapping triangles (or elements) T , which de¬ne the so called triangulation
of the domain which will be denoted by Th . Clearly, „¦ = T . Suppose
T ∈Th
that the maximum length of the edges of the triangles is less than a positive
number h. As shown in Figure 8.6 (left), not any arbitrary triangulation is
allowed. Precisely, the admissible ones are those for which any pair of non
disjoint triangles may have a vertex or an edge in common.


T1
T1
T2
aT
T2
y y 3

FT
1
aT T
1
T1
aT
T 2
T1
T2 x x
T2
0 1
FIGURE 8.6. The left side picture shows admissible (above) and non admissible
(below) triangulations while the right side picture shows the a¬ne map from the
ˆ
reference triangle T to the generic element T ∈ Th

Any element T ∈ Th , of area equal to |T |, is the image through the a¬ne
ˆ
map x = FT (ˆ ) = BT x + bT of the reference triangle T , of vertices (0,0),
x
8.5 Extension to the Two-Dimensional Case 345

ˆ
(1,0) and (0,1) in the x = (ˆ, y ) plane (see Figure 8.6, right), where the

invertible matrix BT and the right-hand side bT are given respectively by

x2 ’ x1 x3 ’ x1
, bT = (x1 , y1 )T ,
BT = (8.34)
y2 ’ y1 y3 ’ y1

(l)
while the coordinates of the vertices of T are denoted by aT = (xl , yl )T
for l = 1, 2, 3.

li (x,y)
li (x,y)
1
1

zi
zi



li (x) li (x)
1
1

zi
zi

FIGURE 8.7. Characteristic piecewise Lagrange polynomial, in one and two space
dimensions. Left, k = 0; right, k = 1

The a¬ne map (8.34) is of remarkable importance in practical computa-
tions, since, once a basis has been generated for representing the piecewise
ˆ
polynomial interpolant on T , it is possible, applying the change of coor-
dinates x = FT (ˆ ), to reconstruct the polynomial on each element T of
x
Th . We are thus interested in devising local basis functions, which can be
fully described over each triangle without needing any information from
adjacent triangles.
For this purpose, let us introduce on Th the set Z of the piecewise interpo-
lation nodes zi = (xi , yi )T , for i = 1, . . . , N , and denote by Pk („¦), k ≥ 0,
the space of algebraic polynomials of degree ¤ k in the space variables x, y
± 
 
 
k
Pk („¦) = p(x, y) = aij x y , x, y ∈ „¦ .
ij
(8.35)
 
 
i,j=0
i+j¤k


Finally, for k ≥ 0, let Pc („¦) be the space of piecewise polynomials of degree
k
¤ k, such that, for any p ∈ Pc („¦), p|T ∈ Pk (T ) for any T ∈ Th . An ele-
k
mentary basis for Pc („¦) consists of the Lagrange characteristic polynomials
k
li = li (x, y), such that li ∈ Pc („¦) and
k

li (zj ) = δij , i, j = 1, . . . , N, (8.36)
346 8. Polynomial Interpolation

where δij is the Kronecker symbol. We show in Figure 8.7 the functions li for
k = 0, 1, together with their corresponding one-dimensional counterparts.

<< . .

. 49
( : 95)



. . >>