Shape-preserving C(2) functional interpolation via parametric cubics

Numerical Algorithms 0 (2005) ?–?

1

Shape-preserving C 2 functional interpolation via parametric cubics P. Lamberti a

a

C. Manni

a

Department of Mathematics, University of Torino - Italy E-mail: {lamberti, manni}@dm.unito.it

The paper proposes a method for the construction of a shape preserving C 2 function interpolating a given set of data. The constructed interpolant is a parametric cubic curve. The shape of the curve can be easily controlled via tension parameters which have an immediate geometric interpretation. The approximation order is investigated and numerical examples are presented. Keywords: Interpolation, Shape preserving, Parametric curves AMS Subject classification: 41A05, 41A25, 41A29, 65D05

1.

Introduction

This paper deals with the problem of generating smooth functions interpolating a prescribed set of data and preserving certain shape properties of the data such as positivity, monotonicity or convexity. In recent years, this problem has received considerable attention because applied mathematicians and CAD users need to avoid inflections and oscillations extraneous to the expected behavior of the data. Numerous methods with such features are available. Among these, the so called tension methods (see for example [3], [4], [13], [17], [19], [21] and references quoted therein) are probably the most well known. Basically they consist of constructing a smooth function as a collection of several patches between the data points that depend, in a local or in a global way, on a set of tension parameters. Selecting these parameters in a suitable way, it is possible to stretch the patches and thus control the shape of the interpolant in order to reproduce the shape of the data. In each scheme the values of tension parameters play a crucial role and have a significant impact on both visual and analytic properties of the constructed

2

P. Lamberti, C. Manni / Shape-preserving functional interpolation

function. The strategy to select the appropriate parameters is very important as well, and a geometric meaning could be useful. However, classical tension methods (rationals, exponentials, polynomials of variable degrees, ...) present some common general problems for the various techniques. Roughly speaking, we can say that rationals and exponentials can produce accurate shape preserving interpolants ([5], [9], [18], [21] and references quoted therein). However, their tension parameters do not have an immediate geometric interpretation. On the other hand, the tension effect in variable degree polynomials can be easily described geometrically by the B´ezier control points (see [2], [3], [4]), but the tension parameters (i.e. the degrees) can assume only integer values and the constructed shape preserving interpolant is, generally, only second order accurate. Recently ([14], [15], [16]) new tension methods for C 1 or C 2 functional interpolation have been proposed. In [14], [15] and [16] the graph of the interpolating function x → s(x) has been considered as the image of a particular parametric curve (x(t), y(t)) whose first component is required to be strictly increasing. The amplitudes of the curve tangent vectors at the extremes of each subinterval determined by the data act as tension parameters, providing an easy control of the shape of the curve and then of the graph of s. Thus, the tension parameters introduced in this way have a simple geometric interpretation. They have a completely local effect and can assume values in a continuous set. Moreover, the shape of the interpolant is easily controlled via its B´ezier control polygon. The above mentioned papers provide different results by considering the interpolating functions as parametric curves having components in particular spaces (piecewise cubics, piecewise quadratics) and studying the properties of the constructed interpolant in the various cases. Nevertheless, the schemes proposed in [14], [15] and [16] are local schemes so, in addition to the data positions, they require that first (and second) derivatives at the data sites, consistent with the behavior of the data, be prescribed to obtain a C 1 (C 2 ) shape preserving interpolant. In practical situations only the data positions are available; thus estimates of derivatives, which are accurate and consistent with the behavior of the data, are necessary (see for example [6], [11]). This is in general a difficult task and the estimated derivatives strongly influence the final quality of the interpolant. A classical alternative is to consider global schemes that do not require the additional knowledge of the derivatives. Global tension methods based on exponentials, rationals or variable degree polynomials ([5], [12], [21], and references

P. Lamberti, C. Manni / Shape-preserving functional interpolation

3

quoted therein) for constructing C 2 shape preserving interpolants are well known. In this paper we describe a global method for the construction of a C 2 shape preserving interpolating function s based on parametric cubic curves. As previously mentioned, in each subinterval determined by the data the graph of s can be interpreted as the image of a planar curve with cubic components whose shape can be immediately controlled by its B´ezier control polygon. The scheme requires the solution of a tridiagonal linear system, as for tension methods based on rationals, exponentials or variable degree polynomials. We show that it is always possible to select tension parameters so that s preserves positivity, monotonicity and/or convexity of the data. In this regard, we provide conditions which are sufficient for positivity and monotonicity and necessary and sufficient for convexity, in order to ensure that the constructed interpolant is shape preserving. A simple automatic algorithm, based on the mentioned conditions, is proposed for selecting values of the tension parameters that provide visually pleasing shape preserving interpolants. We also investigate the approximation order of the presented scheme. It turns out that the constructed interpolant is at least second order accurate if the data are taken from a C 2 function. Moreover, s is fourth order accurate for C 4 functions if the tension parameters are sufficiently close to the data spacing. Summarizing, the proposed scheme is a parametric tension method that allows an automatic choice of the tension parameters ensuring shape preservation and accuracy. It is worth mentioning that the proposed construction presents some analogies with the theory of ν–splines ([10], [20]) introduced in the area of G2 parametric cubics. However, we emphasize that in our method the parametric setting is a technical tool to introduce a space of interpolating functions whose shape can be easily controlled via tension parameters, having a nice geometric meaning, and the method is tailored for functional interpolation. Thus, it is not pertinent to compare the presented scheme with methods for interpolation of parametric curves, where no error estimates are provided since the latter are basically interested in controlling the shape of the curve. It turns out that the parametric setting is, at the same time, the main tool and main drawback of the proposed scheme. In fact, since s is obtained as a parametric curve with piecewise cubic components, its computational cost is comparable with that of classic functional schemes if we only require a plot or a tabulation of the interpolant, but a third order equation must be solved

4

P. Lamberti, C. Manni / Shape-preserving functional interpolation

to compute the value of s at a given fixed value x. However, due to the very simple structure of the x component, see Section 2, this can be efficiently done by classical methods. As an example, since the x component is a piecewise cubic increasing function, a few Newton–Raphson iterations suffice to compute the solution of the above mentioned equation to machine precision. The remainder of this paper is organized into 6 sections. In the next section we describe the construction of s assuming first derivatives at the data sites are given, while in Section 3 we show how to compute the required derivatives in order to ensure s is of class C 2 . Section 4 is devoted to the study of the asymptotic behavior of s as the tension parameters approach their limit configuration and in Section 5 we determine the approximation order of the interpolant as a function of the tension parameters. The shape preserving algorithm is presented in Section 6. We conclude in Section 7 with some numerical examples. Throughout the paper dashes denote derivatives with respect to the x variable, while k.k denotes the infinity norm of vectors, matrices or functions.

2.

Interpolation via parametric cubics

In this section we describe the use of parametric cubic curves for the construction of a function interpolating a set of data. Let the data (xi , fi , di ), i = 0, · · · , N,

(1)

be given, where fi and di are function values and derivatives. We put hi = xi+1 − xi , ∆fi = fi+1 − fi , i = 0, · · · , N − 1, and we assume hi > 0. (j) For i = 0, · · · , N − 1, let hi ≥ 0, j = 0, 1 be free parameters. We consider the parametric cubic curve (see also Figure 1) (

(0)

(1)

(0)

(0)

(0)

(1)

(1)

(1)

Xi (t; hi , hi ) = xi H0 (t) + xi+1 H1 (t) + hi H0 (t) + hi H1 (t), (0) (1) (0) (0) (0) (1) (1) (1) Yi (t; hi , hi ) = fi H0 (t) + fi+1 H1 (t) + hi di H0 (t) + hi di+1 H1 (t), (2) (j) where t ∈ [0, 1] and Hi (t), j = 0, 1 denote the classical Hermite polynomials, i.e. (0)

H0 (t) = 1 − t2 (3 − 2t),

(0)

H1 (t) = t2 (3 − 2t),

P. Lamberti, C. Manni / Shape-preserving functional interpolation (1)

5

(1)

H0 (t) = t − 2t2 + t3 ,

H1 (t) = −t2 (1 − t).

Hence (j) Hi (r)

(j)

dHi (t) = δir δ1j . dt

= δir δ0j ,

(3)

From (3) one can immediately verify that the curve (2) interpolates the points µ

xi+j fi+j

¶

, j = 0, 1

at the extremes of the interval [0, 1] and has tangent vectors µ

(j)

hi (j) hi di+j

¶

, j = 0, 1 (0)

(1)

at the same extremes. Thus, the parameters hi , hi determine the amplitude of the tangent vectors of the curve at the two end points of the interval. Moreover, it is not difficult to see that if (0)

(1)

0 < hi , hi

≤ hi ,

(4)

then d (0) (1) Xi (t; hi , hi ) > 0. dt

(5)

Thus, the first component in (2) implicitly defines a function t = t(x) so that dt = dx (0)

µ

dXi dt

¶−1

, xi ≤ x ≤ xi+1 .

(1)

Then Yi (t; hi , hi ) can be expressed as a function of x. Namely, defining (1)

(0)

s(x) := Yi (t(x); hi , hi ), xi ≤ x ≤ xi+1 , i = 0, · · · , N − 1,

(6)

from the interpolation properties of (2) we have s(xi ) = fi , (0)

s′ (xi +) = hi di

1 (0) hi

(1)

s′ (xi −) = hi−1 di

= di , (0)

(1)

1 (1) hi−1

= di .

Hence, for each pair of values hi , hi satisfying (4), expression (6) defines a function of class C 1 [x0 , xN ] that interpolates the data (1).

6

P. Lamberti, C. Manni / Shape-preserving functional interpolation (0)

(1)

As previously noted, the values hi , hi determine the magnitude of the tangent vectors to the curve (2) at the endpoints and they control the shape of (0) (0) s. To be more specific, since H0 (t) + H1 (t) = 1, we have (

(0)

Xi (t; 0, 0) = xi + hi H1 (t), (0) Yi (t; 0, 0) = fi + ∆fi H1 (t), (j)

i.e., the curve (2) reduces to the line through (xi , fi ), (xi+1 , fi+1 ) if hi (j) 0, 1. On the other hand, if hi = hi , j = 0, 1, we have

= 0, j =

Xi (t; hi , hi ) = xi + hi t and the function s defined in (6) is the classical Hermite cubic polynomial interpo(j) lating (xi+j , fi+j , di+j ), j = 0, 1. Thus, the parameters hi , j = 0, 1 act as tension (j) parameters stretching the curve from the classical Hermite cubic (hi = hi ) to (j) the line segment (hi = 0) (see Figures 1 and 2).5 Assume now that fi = f (xi ). The above described method allows us to locally construct a function s which interpolates f at the data sites xi . Obviously, the quality of the interpolant strongly depends on the values di . If di = f ′ (xi ), we are dealing with a classical Hermite interpolation problem and the method produces C 1 interpolating functions having interesting properties from both the graphical and analytical point of view (see [14], [15] ). However, often in practice only the data (xi , fi ), i = 0, · · · , N

(7)

are available and the required derivatives di must be computed from them. As a possible alternative to the arduous task of estimating derivatives, we can consider a global approach using parametric cubics to construct a C 2 function which interpolates the data (7). This will be done in the next section. In the following we assume that condition (4) holds.

3.

Constructing a C 2 interpolant

Let the data (7) be given. The goal of this section is to construct a function s defined by (2)–(6), which interpolates the data (7) and is of class C 2 . (j) For each choice of the tension parameters hi the previous construction ensures that interpolation conditions are fulfilled. Thus, our problem consists of finding (j) d0 , · · · , dN , depending on the parameters hi , j = 0, 1, i = 0, · · · , N − 1 such

P. Lamberti, C. Manni / Shape-preserving functional interpolation 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

7

−0.4 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−0.4 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

−0.4 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 1. Different behavior of the parametric cubic interpolant for different values of the (0) (1) (0) (1) (0) (1) tension parameters. Top to bottom hi = hi = 1; hi = 0.2, hi = 0.1; hi = hi = 0. Left (t, y(t)), right (x, s(x) = y(t(x))).

that s ∈ C 2 [x0 , xN ]. Denoting by dots the derivatives with respect to t, from (6) we have: y¨(t)x(t) ˙ −x ¨(t)y(t) ˙ d2 s(x) = . 2 3 dx (x(t)) ˙ Hence, imposing continuity of ulations we obtain

d2 s(x) at xi , i = 1, · · · , N − 1, after some manipdx2

ui di−1 + di + vi di+1 = zi

(8)

8

P. Lamberti, C. Manni / Shape-preserving functional interpolation

where (0)

ui =

(0)

hi−1 (hi )2 , σi (1)

(1)

h (hi−1 )2 , vi = i σi 3 3 (1) (0) zi = ∆fi (hi−1 )2 + ∆fi−1 (hi )2 , σi σi (1)

(0)

σi = ai (hi )2 + bi (hi−1 )2 , (0)

ai = 3hi−1 − hi−1 , (1)

bi = 3hi − hi . In order to uniquely determine the values di we need two additional equations that will be obtained by imposing boundary conditions. Classical boundary conditions are natural end conditions: (1)

d0 +

h0

d = (1) 1

3h0 − h0

3∆f0

(0)

, (1)

3h0 − h0

hN −1 3hN −1 −

dN −1 (0) hN −1

+ dN =

3∆fN −1 (0)

3hN −1 − hN −1 (9)

and end tangent conditions: d0 = f0′ ,

′ dN = fN ,

(10)

′ are given in input (see also Section 6). where f0′ , fN One can immediately verify that, for any choice of the previous mentioned boundary conditions, (8) provides a tridiagonal, diagonally dominant system

Ad = z.

(11)

Thus we can state the following (j)

Theorem 1. For any sequence hi , i = 0, · · · , N − 1, j = 0, 1 satisfying (4), there exists a unique s ∈ C 2 [x0 , xN ] defined via (2)–(6) which interpolates the data (7) and satisfies natural or end tangent boundary conditions. (j)

We notice that for hi = hi , j = 0, 1, system (11) reduces to the tridiagonal system for the computation of classical C 2 cubic splines.

P. Lamberti, C. Manni / Shape-preserving functional interpolation

4.

9

Asymptotic behavior

In this section we study the asymptotic behavior of the function s defined (j) by (2), (6) and (11) as the tension parameters hi approach zero. We begin with the following Lemma 2. The coefficient matrix A of system (11) is invertible, kA−1 k ≤ 2 and (j) kzk is bounded independently of hi , i = 0, · · · , N − 1, j = 0, 1. Moreover, the (j) values di , i = 0, · · · , N, obtained from (11) are bounded independently of hi . Proof. Since (4) holds, we have (0)

(0)

(1)

(0)

(1)

(1)

hi−1 (hi )2 + hi (hi−1 )2 hi−1 (hi )2 + hi (hi−1 )2 1 ≤ ≤ , σi σi 2 (1)

h0

3h0 −

(1) h0

1 ≤ , 2

(0)

hN −1 3hN −1 −

(0) hN −1

1 ≤ . 2

Then A is a strictly diagonally dominant matrix and it can be written as A = I+E (0) (1) where I is the identity matrix and E is such that kEk ≤ 12 , for all 0 < hi , hi ≤ hi . Thus kA−1 k ≤ 2. Finally, for i = 1, · · · , N − 1, |zi | =

¯ 3 ¯¯ (1) (0) ¯ ¯∆fi (hi−1 )2 + ∆fi−1 (hi )2 ¯ σi

½ ¾ ¾ ½ ´ 3 ³ |∆fi−1 | |∆fi | 3 |∆fi−1 | |∆fi | (1) 2 (0) 2 ≤ hi (hi−1 ) + hi−1 (hi ) max ≤ max , , σi hi−1 hi 2 hi−1 hi

and from (9) we have 3|∆f0 | 3h0 −

(1) h0

≤

3 |∆f0 | , 2 h0

3|∆fN −1 | 3hN −1 −

(0) hN −1

≤

3 |∆fN −1 | . 2 hN −1

So the claim follows from kdk ≤ kA−1 k kzk.

As mentioned in the Introduction, our main goal is to propose a method for the construction of a smooth interpolant whose shape can be easily controlled via tension parameters, in order to reproduce the shape of the data. In Section 2 it has been shown that the interpolant s approaches the piecewise linear function (j) interpolating the data (7) as the tension parameters hi tend to zero. However,

10

P. Lamberti, C. Manni / Shape-preserving functional interpolation

this is not sufficient to ensure a complete reproduction of the behavior of the data (see Section 6 for formal definitions). In order to achieve our goal, we also have to require that the derivatives di determined by (11) are consistent with the (j) shape of the data as the tension parameters hi tend to zero. This immediately follows from (9) and (10) for d0 , dN . Concerning the derivatives at the interior (0) (1) points we have the following theorem showing that, as hi−1 , hi tend to zero, di ∆fi−1 i tends to a convex combination of ∆f hi and hi−1 . Theorem 3. For i = 1, · · · , N − 1 (1)

(0)

hi (hi−1 )2

∆fi ∆fi−1 hi−1 (hi )2 + . (1) 2 (0) 2 h (0) 2 (1) 2 h (0) (1) i i−1 hi−1 , hi →0 hi (hi−1 ) + hi−1 (hi ) hi−1 (hi ) + hi (hi−1 ) (12) (1) (0) uniformly in hi−1 , hi . lim

di =

Proof. With some manipulation we obtain: (0)

0 ≤ ui ≤

hi−1 (0) → 0 as hi−1 → 0, 2hi−1 (1)

h 0 ≤ vi ≤ i → 0 2hi

(1)

as hi

→ 0,

(1)

(0)

hi (hi−1 )2

∆fi ∆fi−1 hi−1 (hi )2 lim zi = + (1) 2 (0) 2 h (0) 2 (1) 2 h (0) (1) i i−1 hi−1 , hi →0 hi (hi−1 ) + hi−1 (hi ) hi−1 (hi ) + hi (hi−1 ) (1)

(0)

uniformly in hi−1 , hi . Then, the result follows from (8) because dk , k = (j)

0, · · · , N, are bounded independently of hl , l = 0, · · · , N − 1, j = 0, 1. (0)

We note that if

hi

(1) hi−1

=

hi hi−1

the right part of (12) reduces to the following

classical formula for approximation of first derivative hi−1 ∆fi ∆fi−1 hi + . hi + hi−1 hi hi + hi−1 hi−1 5.

Approximation order

We are now interested in studying the approximation power of the constructed interpolanting curve. We recall (see Section 2) that the function s re(j) duces to a linear polynomial in the interval [xi , xi+1 ] if hi = 0, j = 0, 1. On

P. Lamberti, C. Manni / Shape-preserving functional interpolation

11

(j)

the other hand if hi = hi , j = 0, 1, s becomes an algebraic cubic interpolant in the same interval. Thus, we expect to obtain a result on the approximation order that reduces to O(h2 ) and O(h4 ) in the above cases for data belonging to a function f ∈ C 2 and C 4 , respectively. For sake of clarity it is convenient to denote the second component of the (0) (1) curve (2) by Yi (t; hi , hi , di , di+1 ) when the dependence on the parameters (j) hi , j = 0, 1 and the derivatives di , di+1 has to be emphasized. In addition, let us put fi′ = f ′ (xi ), i = 0, · · · , N and h=

max

0≤i≤N −1

hi .

In order to show the main result of this section, we need the following Lemma 4. Let di , i = 0, · · · , N be determined by (8) and (10). C 4 [x0 , xN ], then

If f ∈

| di − fi′ |≤ 2kzk, i = 0, · · · , N

(13)

where 1 ′′′ 7 ′′ (j) (j) kzk ≤ 2kf k max(hi − hi ) + kf k max hi (hi − hi ) + kf (4) kh3 , i,j i,j 2 24

(14)

i = 0, · · · , N − 1, j = 0, 1. Moreover, if f ∈ C 2 [x0 , xN ], then ′′

| di − fi′ |≤ 5kf kh, i = 0, · · · , N.

(15)

Proof. We want to estimate the quantities di − fi′ that satisfy the system (0)

(0)

(1)

(1)

h (hi−1 )2 hi−1 (hi )2 ′ ′ (di−1 − fi−1 ) + (di − fi′ ) + i (di+1 − fi+1 ) = zi, σi σi i = 1, · · · , N − 1 d0 − f0′ = 0,

′ dN − fN = 0,

where i 1h (1) (0) (0) (0) (1) (1) ′ ′ 3∆fi (hi−1 )2 + 3∆fi−1 (hi )2 − σi fi′ − hi−1 (hi )2 fi−1 − hi (hi−1 )2 fi+1 . zi = σi

12

P. Lamberti, C. Manni / Shape-preserving functional interpolation

′ , f′ If we assume f ∈ C 4 [x0 , xN ] and we expand ∆fi , ∆fi−1 , fi+1 i−1 at x = xi , then we obtain 1 ′′ (1) (1) (0) (0) zi = {fi [hi (hi − hi )(hi−1 )2 − hi−1 (hi−1 − hi−1 )(hi )2 σi 1 1 (1) 2 (0) (0) (1) + ((hi−1 ) − h2i−1 )(hi )2 + (h2i − (hi )2 )(hi−1 )2 ] 2 2 1 ′′′ (1) (1) (0) (0) + fi [h2i (hi − hi )(hi−1 )2 + h2i−1 (hi−1 − hi−1 )(hi )2 ] 2 1 (1) (1) + [(hi−1 )2 (3h4i f (4) (ξi ) − 4hi h3i f (4) (ηi )) 24 (0) (0) + (hi )2 (−3h4i−1 f (4) (ξ i ) + 4hi−1 h3i−1 f (4) (η i ))]}

where ξi , ηi ∈ [xi , xi+1 ], ξ i , η i ∈ [xi−1 , xi ]. Since (0)

(0)

(1)

(1)

(0)

(1)

σi ≥ max{(3hi−1 −hi−1 )(hi )2 , (3hi −hi )(hi−1 )2 } ≥ 2 max{hi−1 (hi )2 , hi (hi−1 )2 }, after some algebraic manipulations we obtain (14). Thus, the result (13) follows from Lemma 2. Now, repeating the same proof for f ∈ C 2 [x0 , xN ], we obtain 5 ′′ kzk ≤ kf kh 2 from which the result (15) follows immediately. (j)

We remark that | fk′ − dk |= O(h3 ), k = 0, · · · , N if hi = hi , j = 0, 1, i = 0, · · · , N − 1, as expected from the approximation properties of C 2 cubic interpolanting splines with end tangent conditions. Indeed it is now possible to prove the following Theorem 5. Let s be determined by (2), (6), (8) and (10). If f ∈ C 4 [x0 , xN ], then for any x ∈ [x0 , xN ] (j)

(j)

(j)

| f (x)−s(x) |≤ D0 max(hi −hi )+D1 max(hi −hi )h+D2 max(hi −hi )h2 +D3 h4 , i,j

i,j

i,j

(16) i = 0, · · · , N − 1, j = 0, 1, where 32 ′′ 8 ′′′ 4 D0 = kf ′ k, D1 = kf k, D2 = kf k, D3 = 3 27 27

µ

14 1 + kf (4) k. 384 81 ¶

Moreover, if f ∈ C 2 [x0 , xN ], we have | f (x) − s(x) |≤

379 ′′ 2 kf kh . 216

(17)

P. Lamberti, C. Manni / Shape-preserving functional interpolation

13

Proof. Let x ∈ [xi , xi+1 ] be given. Let us assume that f ∈ C 4 [x0 , xN ]. Then, (0) (1) from (5) there exist t, tˆ such that Xi (t; hi , hi ) = Xi (tˆ; hi , hi ) = x. We can split the error as follows (0)

(1)

| f (x) − s(x) | = | f (x) − Yi (t; hi , hi , di , di+1 ) | ′ ≤ | f (x) − Yi (tˆ; hi , hi , fi′ , fi+1 )| ′ ′ )| + | Yi (tˆ; hi , hi , fi′ , fi+1 ) − Yi (t; hi , hi , fi′ , fi+1 (0)

(1)

′ + | Yi (t; hi , hi , fi′ , fi+1 ) − Yi (t; hi , hi , di , di+1 ) |=: A + B + C. ′ ) is the Hermite cubic interpolating f As noticed in Section 2, Yi (t; hi , hi , fi′ , fi+1 at xi , xi+1 . Thus, it follows that

1 kf (4) kh4i . 384 Now, by applying Lagrange theorem, we have the following upper bound for B A≤

d ′ B = | (tˆ − t) Yi (t; hi , hi , fi′ , fi+1 )|t=τ | dt ¯ ¯ h ∆f d ¯ d (1) d (1) i¯¯ (0) i ′ H1 (τ ) + fi′ H0 (τ ) + fi+1 H1 (τ ) ¯ = ¯¯hi (tˆ − t) hi dt dt dt 7 ≤ kf ′ khi | tˆ − t | 2 where τ ∈ [0, 1]. By subtracting the two expressions for x − xi x − xi = Xi (tˆ; hi , hi ) − xi = hi tˆ and (0)

(1)

Xi (t; hi , hi ) − xi = (0)

(1)

(1)

(0)

(1)

(1)

(1)

hi H1 (t) + hi H0 (t) + hi H1 (t) + (hi − hi )H0 (t) + (hi − hi )H1 (t), we obtain (0) (1) (1) (1) hi (t − tˆ) = (hi − hi )H0 (t) + (hi − hi )H1 (t)

which gives B≤

14 ′ 28 (0) (1) (j) kf k[(hi − hi ) + (hi − hi )] ≤ kf ′ k max(hi − hi ). j 27 27

Finally from Lemma 4 we can write (0)

(1)

(1)

(1)

′ − hi di+1 )H1 (t) C = (hi fi′ − hi di )H0 (t) + (hi fi+1

14

P. Lamberti, C. Manni / Shape-preserving functional interpolation (0)

(0)

(1)

= [fi′ (hi − hi ) + hi (fi′ − di )]H0 (t) (1)

(1)

(1)

′ ′ +[fi+1 (hi − hi ) + hi (fi+1 − di+1 )]H1 (t) 8 (j) ≤ {kf ′ k max(hi − hi ) + 2kzkh}. j 27

Then the result (16) follows from (14). (0) (1) If f ∈ C 2 ([x0 , xN ]), there exist t, t˜ such that Xi (t; hi , hi ) = Xi (t˜; 0, 0) = (0) (0) x. In fact, since Xi (t; 0, 0) = xi + hi H1 (t) and H1 is invertible, we obtain immediately the invertability of Xi (t; 0, 0). Then we can split the error as follows (0)

(1)

| f (x) − s(x) | = | f (x) − Yi (t; hi , hi , di , di+1 ) | ≤ | f (x) − Yi (t˜; 0, 0, di , di+1 ) | (0) (1) + | Yi (t˜; 0, 0, di , di+1 ) − Yi (t; hi , hi , di , di+1 ) |=: E + F.

Since Yi (t; 0, 0, di , di+1 ) is the line through (xi , f (xi )), (xi+1 , f (xi+1 )), it follows that 1 ′′ E ≤ kf kh2i . 8 For F we note that from (0)

(0)

(1)

(1)

(1)

(0)

x = xi + hi H1 (t) + hi H0 (t) + hi H1 (t) = xi + hi H1 (t˜) we can write (0)

(0)

(0)

(1)

(1)

(1)

hi [H1 (t˜) − H1 (t)] = hi H0 (t) + hi H1 (t) and hence F ≤|

³ ∆f

i

hi

´

(0)

(1)

− fi′ hi H0 (t) + (0)

(1)

³ ∆f

i

hi

´

(1)

(1)

′ − fi+1 hi H1 (t) | (1)

(1)

′ − di+1 )hi H1 (t) | + | (fi′ − di )hi H0 (t) + (fi+1 4 ′′ ≤ (kf kh2i + 2 max | fk′ − dk | hi ). k 27

Thus, from Lemma 4 we obtain (17). We remark that, if f ∈ C 2 , the result of Theorem 5 still holds even if s is determined by natural end conditions.

P. Lamberti, C. Manni / Shape-preserving functional interpolation

6.

15

Preserving the shape of the data

In this section we show how the shape of s can be controlled according to the behavior of the data. For sake of clarity we state the following Definition 6. We say that the data (7) are (a) positive in [xi , xi+1 ], i = 0, · · · , N − 1, if fi , fi+1 > 0; (b) monotone increasing in [xi , xi+1 ], i = 1, · · · , N − 2, if fi−1 < fi < fi+1 < fi+2 ; (c) convex in [xi , xi+1 ], i = 1, · · · , N − 2, if ∆fi ∆fi+1 ∆fi−1 < < . hi−1 hi hi+1 As far as the end intervals are concerned, in order to define the “monotonicity” and/or the “convexity” of the data we have to take into account the considered boundary conditions. If natural end conditions (9) are used we state the following Definition 7. We say that the data (7) are (b1 ) monotone increasing in [x0 , x1 ] if f0 < f1 < f2 ; (c1 ) convex in [x0 , x1 ] if ∆f1 ∆f0 < . h0 h1 If tangent end conditions (10) are used the set of input data becomes ′ (xi , fi ), i = 0, · · · , N, f0′ , fN

and we state the following Definition 8. We say that the data (18) are (b1 ) monotone increasing in [x0 , x1 ] if f0 < f1 < f2 , and f0′ ≥ 0;

(18)

16

P. Lamberti, C. Manni / Shape-preserving functional interpolation

(c1 ) convex in [x0 , x1 ] if ∆f1 ∆f0 ∆f0 < and f0′ < . h0 h1 h0 Similar definitions can be stated in the interval [xN −1 , xN ] and for negative, decreasing and concave data. In each subinterval [xi , xi+1 ] the curve (2) can be rewritten in the B´ezier form (see for example [10]) µ

(0)

(1)

Xi (t; hi , hi ) (0) (1) Yi (t; hi , hi )

¶

=

à ! 3 X 3 l=0

l

(i)

Pl tl (1 − t)3−l

where (i) P0 (i) P1

=

(i) P0

=

à !

xi , fi

(i) P3

=

Ã

!

xi+1 , fi+1

à !

Ã

1 (1) 1 1 (0) 1 (i) (i) , P2 = P3 − hi + hi 3 3 di di+1

!

(19)

are the B´ezier control points and the polygonal line they form is the B´ezier control polygon of (2) (see Figure 2). It is well known that the shape of the curve (2) can be easily controlled via its B´ezier control points. To be more precise we have the following result that for sake of simplicity we state for positive, increasing and convex data only. (i)

(i)

Proposition 9. Let us denote by Pl,x , (Pl,y ), l = 0, · · · , 3, the x (y) components of the control points (19). For i = 0, · · · , N − 1, (a) if (i)

Pl,y > 0, l = 0, 1, 2, 3,

(20)

then s(x) > 0, x ∈ [xi , xi+1 ]; (b) if (i)

(i)

Pl+1,y − Pl,y ≥ 0, l = 0, 1, 2,

(21)

then s′ (x) ≥ 0, x ∈ [xi , xi+1 ]; (c) s′′ (x) ≥ 0, x ∈ [xi , xi+1 ] if and only if (i)

(i)

(i)

(i)

Pl+2,y − Pl+1,y Pl+2,x − Pl+1,x

≥

(i)

(i)

(i)

(i)

Pl+1,y − Pl,y

Pl+1,x − Pl,x

, l = 0, 1.

(22)

P. Lamberti, C. Manni / Shape-preserving functional interpolation 15

15

14

14

13

13

12

12

11

11

10

10

9

9

8

8

7

7

6

9

9.2

9.4

9.6

9.8

10

10.2

10.4

10.6

10.8

11

6

9

9.2

9.4

9.6

9.8

10

10.2

10.4

17

10.6

10.8

11

Figure 2. Curve (2) and its control polygon for different values of tension parameters.

Proof. The results concerning positivity and monotonicity immediately follow from the shape preserving properties of the B´ezier–Bernstein representation (see for example [8], [10]). As far as convexity is concerned we notice that if (22) holds then the B´ezier control polygon of the curve is convex and we can still refer to the shape preserving properties of the B´ezier–Bernstein representation. Vice versa, from the properties of the derivatives of B´ezier curves (see for example [10]) it is not difficult to verify that (22) are equivalent to s′′ (xi+l ) ≥ 0, l = 0, 1.

In addition, we remark that, due to the shape preserving properties of the B´ezier–Bernstein representation, (see [8]) the control polygon of the curve (2) exactly mimics the behavior of s. Thus, in each subinterval [xi , xi+1 ] s can have at most two local extrema and/or one inflection point as classical cubic polynomials. Let us consider the set of indices I = {0, 1, · · · , N − 1}. From the previous proposition and from the results of Section 4 we have the following Theorem 10. Let s be constructed via (2), (6), (8) with boundary conditions (9) or (10). Let R > 0 be given. For i = 0, · · · , N − 1, if the data are

18

P. Lamberti, C. Manni / Shape-preserving functional interpolation

(a) positive/negative in [xi , xi+1 ] and (0)

(1)

hi , hi

are small enough, then s is positive/negative in [xi , xi+1 ]; (b) monotone increasing/decreasing in [xi , xi+1 ] and (0)

(1)

hl , hk , l ∈ {i − 1, i} ∩ I, k ∈ {i, i + 1} ∩ I, are small enough, then s is monotone increasing/decreasing in [xi , xi+1 ]; (c) convex/concave in [xi , xi+1 ] and (0)

(1)

hk , hk , k ∈ {i − 1, i, i + 1} ∩ I, are small enough with (1)

hk

(0)

hk+1

(0)

< R,

hl

(1)

hl−1

< R, k ∈ {i − 1} ∩ I, l ∈ {i + 1} ∩ I,

(23)

then s is convex/concave in [xi , xi+1 ]. Proof. Let us assume, as an example, that natural end conditions are used. If the data are positive in [xi , xi+1 ], then fi > 0 and fi+1 > 0. Hence, from (0) (1) (19), (20) hold if hi , hi are small enough. Thus, from Proposition 9, s is positive in [xi , xi+1 ]. Let us assume now that the data are monotone increasing in [xi , xi+1 ]. From (0) (1) (9) and (12) if hl , hi , l ∈ {i − 1} ∩ I, are small enough, then di > 0. Similarly, (1) (0) if hi , hk , k ∈ {i + 1} ∩ I, are small enough, then di+1 > 0. Hence, from (19), (0) (1) (21) hold if hi , hi are small enough. Thus, from Proposition 9, s is monotone increasing in [xi , xi+1 ]. Finally, let us assume that the data are convex in [xi , xi+1 ]. From (9) and (0) (1) (0) (1) (12), if hl , hi , hi , hk , l ∈ {i − 1} ∩ I, k ∈ {i + 1} ∩ I, are small enough and (0) (1) i are small (23) hold, then di < ∆f hi < di+1 . So, from (19), (22) hold if hi , hi enough. Thus, from Proposition 9, s is convex in [xi , xi+1 ]. We can conclude in a similar way if end tangent conditions are used. The following algorithm performs the construction of the parametric cubic interpolant with shape preserving properties. It has been organized to sequentially treat positivity, monotonicity and convexity. However, as it can be seen from the proof of Theorem 10, the various constraints can be considered in any order.

P. Lamberti, C. Manni / Shape-preserving functional interpolation

19

Algorithm 11. 0.

let the data (7) or (18) be given (j)

1. 2. 3. 4. 4.1

put hi = hi , j = 0, 1, i = 0, · · · , N − 1, construct the coefficient matrix A and the known vector z of system (11) with boundary conditions (9) or (10) solve system (11) construct the B´ezier control points automatic shape control: for i = 0, · · · , N − 1, if the data are positive in [xi , xi+1 ] and conditions (20) are not satisfied: (1)

4.1.1 decrease hk , k ∈ {i − 1, i} ∩ I, (0)

4.2

4.1.2 decrease hl , l ∈ {i, i + 1} ∩ I, 4.1.3 go to step 1 if the data are monotone increasing in [xi , xi+1 ] and conditions (21) are not satisfied: (1)

4.2.1 decrease hk , k ∈ {i − 2, i − 1, i, i + 1, } ∩ I, (0)

4.3

4.2.2 decrease hl , l ∈ {i − 1, i, i + 1, i + 2} ∩ I, 4.2.3 go to step 1 if the data are convex in [xi , xi+1 ] and conditions (22) are not satisfied: (1)

4.3.1 decrease hk , k ∈ {i − 2, i − 1, i, i + 1, } ∩ I, (0)

4.3.2 decrease hl , l ∈ {i − 1, i, i + 1, i + 2} ∩ I, 4.3.3 go to step 1. (j)

The previous algorithm requires a reduction of the parameters hk in case some shape constraints are not satisfied. Of course, various strategies can be used to implement such a reduction. For sake of completeness, we detail in the following Algorithm 12 the strategy we have used if conditions (22) are not satisfied in a interval where the data are convex. Completely similar strategies have been used for positivity and monotonicity. We emphasize that in Algorithms 11 the term decrease means that a value not exceeding the previous one must be considered (see also Algorithm 12). We also point out that Algorithms 11 and 12 have been structured to ensure that the ratio between the tension parameters adjacent to each point xk (see (23)) remains bounded.

20

P. Lamberti, C. Manni / Shape-preserving functional interpolation

Algorithm 12. 0. 1. 1.1

let 0 < ε, ε1 < 1, ̺ > 0 be given let i ∈ {0, . . . , N − 1} be given if the data are convex in [xi , xi+1 ] and conditions (22) are not satisfied: (0)

1.1.1 hi 1.1.2 1.1.3

1.1.4

7.

(1) hi

(0)

(1)

= εhi , hl

(1)

(0) hl hi ), (0) (1) hk min (hk , ̺hi hi ),

= min (hl , ̺hi

(1) (0) = εhi , hk = i if di > ∆f hi : (0) (0) hl = ε1 hl , l ∈ {i − 1} ∩ I, (1) (1) (0) hk = min (hk , ̺hl hhkl ), k ∈ i if di+1 < ∆f hi : (1) (1) hk = ε1 hk , k ∈ {i + 1} ∩ I, (0) (0) (1) hl = min (hl , ̺hk hhkl ), k ∈

l ∈ {i − 1} ∩ I, k ∈ {i + 1} ∩ I,

{i − 2} ∩ I, l ∈ {i − 1} ∩ I,

{i + 1} ∩ I,

l ∈ {i + 2} ∩ I.

Numerical examples and conclusions

In this section we present the performances of the proposed method over a pair of classical examples. In both cases natural end conditions have been considered and the following values for the input parameters have been used ε = .9, ε1 = .99, ̺ = 10. The first set of data is the so called Akima data taken from [1] (see Table 1). Figures 3 and 4 show the behavior of the interpolant s (right column) and of its control polygon (left column) for different values of the tension parameters (j) (0) (1) hi . In Figure 3 we have hi = hi = hi , so we are dealing with the classical C 2 cubic spline interpolating the data. In order to remove extraneous oscillations from the graph of the interpolant, Algorithm 11 reduces the values of the tension parameters in some intervals (see Figure 4). In Figures 5 and 6 we depict the first and second derivative of the C 2 cubic spline and of s in order to show the shape preserving performances of s. Figure 7 shows a magnification of the second derivative of the above mentioned functions. On this regard, we notice that, according to Definitions 6, 7, 8, the data are neither monotone nor convex in the interval [0, 8]. Thus, only the positivity of the data can be preserved in this interval.

P. Lamberti, C. Manni / Shape-preserving functional interpolation 90

90

80

80

70

70

60

21

60

50 50 40 40 30 30 20 20 10 10 0 0

5

10

15

0

0

Figure 3. Akima data: 90

90

80

80

70

70

60

5

(0) hi hi

=

(1) hi hi

10

15

10

15

= 1.

60

50 50 40 40 30 30 20 20 10 10 0 0

5

Figure 4.

10

15

0

0

5

Akima data:

(0)

hi hi

(1)

= [1 1 1 1 1 .4135 .4783 1 .5314 1],

hi hi

= [1 1 1 1 1 .4305 .4413 1 .5314 1].

Table 1 Akima data. xi fi

0 10

2 10

3 10

5 10

6 10

8 10

9 10.5

11 15

12 50

14 60

15 85

Table 2 Radiochemical data. xi 7.99 8.09 8.19 8.7 9.2 10 12 15 20 fi 0 .276429e-6 .437498e-3 .169183 .469428 .943740 .998636 .999919 .999994

In the second example we consider the radiochemical data taken from [7] (see Table 2). Figure 8 shows the behavior of the interpolant s (right) and of the classical C 2 cubic spline interpolating the data (left).

22

P. Lamberti, C. Manni / Shape-preserving functional interpolation

45

45

40

40

35

35

30

30

25

25

20

20

15

15

10

10

5

5

0

0

−5

−5 0

5

10

15

0

5

10

15

Figure 5. Akima data. First derivative: C 2 cubic spline (left) and s (right). 100

100

50

50

0

0

−50

−50

−100

−100

0

5

10

15

0

5

10

15

Figure 6. Akima data. Second derivative: C 2 cubic spline (left) and s (right). 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

0

1

2

3

4

5

6

7

8

9

0

1

2

3

4

5

6

7

8

9

Figure 7. Akima data. Magnification of second derivative: C 2 cubic spline (left) and s (right).

P. Lamberti, C. Manni / Shape-preserving functional interpolation

23

1.2

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2

0.2

0

0

8

10

12

14

16

18

Figure 8. Radiochemical data. Left:

20

(0) hi hi

=

(0)

hi hi

= [.1351 .5905 1 1 .3487 .2288 .3487 .2824],

8

(1) hi hi (1) hi hi

10

12

14

16

18

20

= 1 Right: = [.1351 .5905 1 1 .3487 .2069 .3217 .2606].

We conclude this section with some computational remarks. As mentioned in the Introduction, the explicit analytic expression of s is not known. However, we remark that from (2) in each subinterval the graph of s lies over an algebraic curve of degree 3 and genus 0 (see for example [22]). The implicit equation of such a curve can be easily obtained by classical implicitization techniques based on resultants (see for example [23] and references quoted therein). The implicit equation can be in principle inverted in order to obtain the explicit analytic expression of s. Nevertheless, such a procedure is not useful from the computational point of view. Finally, we mention that the method we have proposed can be extended to produce C 2 shape preserving interpolants for bivariate tensor product data. This extension is currently under study. Acknowledgments The authors would like to thank the referees for their careful reading of the paper and their useful comments and remarks. References [1] H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, J. Assoc. Comput. Mech. 17 (1970) 589–602. [2] P. Costantini, On monotone and convex spline interpolation, Math. Comput. 46 (1986) 203–214.

24

P. Lamberti, C. Manni / Shape-preserving functional interpolation

[3] P. Costantini, An algorithm for computing shape–preserving interpolating splines of arbitrary degree, J. Comput. Appl. Math. 22 (1988) 89–136. [4] P. Costantini, Curve and surface construction using variable degree polynomial splines, Comp. Aided Geom. Design 17 (2000) 419–446. [5] R. Delbourgo and J. A.Gregory, C 2 rational quadratic spline interpolation to monotonic data, I.M.A. J. Numer. Anal. 3 (1983) 141–152. [6] R. L. Dougherty, A. Edelman and J. M.Hyman, Nonnegativity, monotonicity or convexitypreserving cubic and quintic Hermite interpolation, Math. Comp. 52 (1989) 471–494. [7] F. N. Fritsch and R. E. Carlson, Monotone piecewise cubic interpolation, S.I.A.M. J. Numer. Anal. 17 (1980) 238–246. [8] T. N. T. Goodman, Total positivity and the shape of curves, in: Total Positivity and its Applications, M. Gasca and C. A. Micchelli (eds.), Kluwer, Dordrecht, 1996, pp. 157–186. [9] C. Grandison, Behaviour of exponential splines as tensions increase without bound, J. Approx. Theory 89 (1997) 289–307. [10] J. Hoschek and D. Lasser, Fundamentals of Computer Aided Geometric Design, A. K. Peters, Ltd, Wellesley, Massachusetts, (1993). [11] H. T. Huynh, Accurate monotone cubic interpolation, S.I.A.M. J. Numer. Anal. 30 (1993) 57–101. [12] P. D. Kaklis and D. G. Pandelis, Convexity-preserving polynomial splines of non-uniform degree, I.M.A. J. Numer. Anal. 10 (1990) 223–234. [13] B. I. Kvasov, Algorithms for shape preserving local approximation with automatic selection of tension parameters, Comp. Aided Geom. Design 17 (2000) 17–37. [14] C. Manni, C 1 comonotone Hermite interpolation via parametric cubics, J. Comp. Appl. Math. 69 (1996) 143–157. [15] C. Manni, Parametric shape preserving Hermite interpolation by piecewise quadratics, in: Advanced Topics in Multivariate Approximation, F. Fontanella, K. Jetter and P. J. Laurent (eds.), World Scientific Publishing Co. Pte., Singapore, 1996, pp. 211–226. [16] C. Manni, On shape preserving C 2 Hermite interpolation, BIT 41 (2001) 127–148. [17] C. Manni and P. Sablonni`ere, Monotone interpolation of order 3 by C 2 cubic splines, I.M.A. J. Numer. Anal. 17 (1997) 305-320. [18] M. Marusic and M. Rogina, Sharp error bounds for interpolating splines in tension, J. Comp. Appl. Math. 61 (1995) 205-223. [19] B. Mulansky and J. W. Schmidt, Constructive methods in convex C 2 interpolation using quartic splines, Numer. Algorithms 12 (1996) 111–124. [20] G. M. Nielson, Some piecewise polynomial alternatives to spline under tension, in: Computer Aided Geometric Design, R. E. Barnill and R. F. Riesenfeld (eds.) Academic Press 1974, pp. 209–235. [21] D. G. Schweikert, An interpolation curve using a spline in tension, J. Math. Phys. 45 (1966) 312–317. [22] R.J. Walker, Algebraic Curves, Dover, New York (1962). [23] G. Wang and T. W. Sederberg, Verifying the implicitization formulae for degree n rational B´ezier curves, J. Comput. Math. 17 (1999) 33–40.