- Email: [email protected]
A bivariate rational cubic interpolating spline with biquadratic denominator
Applied Mathematics and Computation 264 (2015) 366–377
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A bivariate rational cubic interpolating spline with biquadratic denominator Youtian Tao a,b,c,∗, Dongyin Wang a a
Department of Mathematics, Chaohu College, Chaohu 238000, China School of Mathematics Sciences, USTC, Hefei 230026, China c Anhui Fuhuang Steel Structure, Chaohu 238076, China b
a r t i c l e
i n f o
Keywords: Bivariate rational interpolating spline Shape parameter Bounded property Error estimate Symmetry
a b s t r a c t A bivariate rational bicubic interpolating spline (BRIS) with biquadratic denominator and six shape parameters is constructed in a rectangle domain. The C1 continuous condition of BRIS discussed. BRIS is proved to be bounded and its error is estimated. In the case of the equally spaced knots, the matrix expression and symmetry of BRIS are presented. Some properties of the basis of BRIS are given. In the end, a numerical example is given to illustrat the effect of the shape parameters on the shape of BRIS surface. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The construction method of curves and surfaces and their mathematical description is a key issue in CAGD. There are many ways [5–7,9–11,19,22,24,25,27–29] to deal with this problem, such as the polynomial spline method, Bézier spline method and NURBS method. However, most of the polynomial spline methods are the interpolating cases, and their local shapes can not be modified for the interpolating surfaces while interpolating data is unchanged. NURBS and Bézier methods are the noninterpolating cases, that is to say, the constructed curve and surface do not fit with the given data and the given points are the control points. Therefore, when we construct the interpolating functions required for CAGD, we should consider the following cases: 1) The expressions of interpolating functions are simple and explicit. 2) The parameters of constructed curves and surfaces can be modified without changing the given data. Recently, there has been many works [3,4,8,12,18,20,21,26] about univariate rational spline interpolation with parameters. Some univariate rational interpolating splines are generalized to bivariate rational interpolating splines which expressions with parameters are simple and explicit. Some interesting results are presented. In [2,13–17,23,30], the authors constructed several bivariate interpolating splines over rectangular mesh, derived some properties such as the sufficient conditions of downconstrained and up-constrained for the shape control, the matrix expression, bounded property, stability, convexity control, the preserving positivity. In this paper, motivated by [2,15], we generate the bivariate rational interpolating spline with four shape parameters in [15] to the case with six shape parameters (BRIS). We present the definition of BRIS and discuss the effect of the parameter choice on the curve shape in Section 2. We give the sufficient condition that BRIS is C1 continuous in Section 3. We prove the bounded property and analyse the error of BRIS in Section 4. In the case of the equally spaced knots, we present the matrix expression of
∗
Corresponding author at: Chaohu College, Department of Mathematics, Jinzhai Rd., Chaohu 238000, China. Tel: +8615156520160. E-mail address: [email protected] (Y. Tao).
http://dx.doi.org/10.1016/j.amc.2015.04.100 0096-3003/© 2015 Elsevier Inc. All rights reserved.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
367
Table 1 Data points from [1]. i
1
2
3
4
5
6
7
8
9
xi fij
0 2
0.25 0.6
0.5 0.1
1 0.13
1.5 1
2 0.5
2.5 1.1
3 0.25
4 0.2
BRIS, prove that it is symmetry, and discuss some properties of the basis of BRIS in Section 5. We give a numerical example of BRIS to illustrat the effect of the shape parameters on the shape of BRIS surface in Section 6. 2. Rational interpolating spline Let D be the rectangular domain [a, b, c, d], {(xi , yi ), fij , i = 1, 2, . . . , n + 1; j = 1, 2, . . . , m + 1} be a given set of data points, where a = x1 < x2 < < xn+1 = b; c = y1 < y2 < < ym+1 = d. Set
fij = f (xi , yj ), hi = xi+1 − xi , lj = yj+1 − yj . Denote Dij = [xi , xi+1 yj , yj+1 ], i = 1, 2, . . . , n, j = 1, 2, . . . , m. For any point (x, y) Dij , let
θ = (x − xi )/ hi , η = (y − yj )/ lj . For each yj , j = 1, 2, . . . , m + 1, we construct the x-direction interpolation curve as follows
Pij∗ (x) =
p∗ij (x) q∗ij (x)
, i = 1, 2, . . . , n − 1,
(1)
where
p∗ij (x) =
αij∗ fij (1 − θ )3 + Vij∗ θ (1 − θ )2 + Wij∗ θ 2 (1 − θ ) + γij∗ fi+1,j θ 3
q∗ij (x) =
αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
and
βij∗ fij + αij∗ fi+1,j = (βij∗ + γij∗ )fi+1,j − γij∗ hi ∗i+1
Vij∗ = Wij∗
with αij∗ , βij∗ , γij∗ positive, and ∗i =
fi+1,j −fij . hi
We can prove that
Pij∗ (xi ) = fij , Pij∗ (xi+1 ) = fi+1,j , Pij∗ (xi ) = ∗ij , Pij∗ (xi+1 ) = ∗i+1,j . The interpolation (1) is called the rational cubic interpolating spline. It is clear that the interpolation is local in the interval [xi , xi+1 ] and depends on the data at three points {(xr , yj ), frj }, r = i, i + 1, i + 2 and the shape parameters αij∗ , βij∗ , γij∗ . It is interesting to note that the interpolation (1) becomes a standard cubic Hermite spline with the values of the shape parameters αij∗ = 1, βij∗ = 2, γij∗ = 1. we now illustrate the mathematical and graphical effects of the shape parameters αij∗ , βij∗ , γij∗ on the shape of a curve. The three free parameters can be exploited properly to modify the shape of curve according to the designer’s choice. We rewrite (1) as follows:
Pij∗ (x) =
αij∗ (1 − θ )2 (fij (1 − θ ) + fi+1,j θ ) + βij∗ θ (1 − θ )(fij (1 − θ ) + fi+1,j θ ) + γij∗ θ 2 (fi+1,j − hi (1 − θ )i+1,j ) αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(2)
By the computation for Eq. (2) we have the following formulas:
lim P ∗ (x) = fij (1 − θ ) + fi+1,j θ αij∗ →∞ ij
(3)
lim Pij∗ (x) = fij (1 − θ ) + fi+1,j θ
(4)
βij∗ →∞
lim P ∗ (x) = fi+1,j − (1 − θ )hi i+1,j γij∗ →∞ ij
(5)
We can see from Eqs. (3) and (4) that the increase in the shape parameter αij∗ or βij∗ reduces the rational interpolating spline (2) to the straight line fij (1 − θ ) + fi+1, j θ , and see from Eq. (5) that the increase in γij∗ reduces (2) to another straight line fi+1, j − (1 − θ )hi i+1, j . According to data points in Table 1 from [1], by choosing different shape parameters, we can observe the corresponding shape changes of the interpolating curves in Figs. 1–4. Each piecewise in Fig. 1 slopes heavily at the left-hand side. Each piecewise in Fig. 2 and 3 incline to a line segment. The curve in Fig. 4 is a cubic Hermite interpolating spline. We remark that in order to plot the 8th interpolating spline curve, we add the point {x10 , f10 } = {4.5, 0.2}.
368
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
2.0
1.5
1.0
0.5
0.0
0.5
0
1
2
3
4
Fig. 1. The rational spline interpolation Eq. (1) with α = 1, β = 1, γ = 100. ∗ ij
∗ ij
∗ ij
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 2. The rational spline interpolation Eq. (1) with α = 10, β = 1, γ = 1. ∗ ij
∗ ij
∗ ij
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 3. The rational spline interpolation Eq. (1) with αij∗ = 1, βij∗ = 10, γij∗ = 1.
We now use the x-direction interpolation Pij∗ (x), i = 1, 2, . . . , n − 1; j = 1, 2, . . . , m + 1 to construct the bivariate rational interpolating function spline (BRIS) in D. For each pair (i, j), i = 1, 2, . . . , n − 1; j = 1, 2, . . . , m − 1, we define a bivariate interpolating spline Pij (x, y) in Dij as follows
Pij (x, y) =
pij (x, y) , qij (x, y)
pij (x, y) =
∗ αij Pij∗ (x)(1 − η)3 + Vij η(1 − η)2 + Wij η2 (1 − η) + γij η3 Pi,j+1 (x)
qij (x, y) =
αij (1 − η)2 + βij η(1 − η) + γij η2
where
(6)
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
369
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 4. The rational spline interpolation Eq. (1) with αij∗ = 1, βij∗ = 2, γij∗ = 1.
and ∗ βij Pij∗ (x) + αij Pi,j+1 (x) ∗ Wij = (βij + γij )Pi,j+1 (x) − γij lj i,j+1 (x)
Vij =
∗ with α ij , β ij , γ ij positive, and ij (x) = (Pi,j+1 (x) − Pij∗ (x))/ lj . Pij (x, y) is called a bivariate bicubic rational interpolating spline with a bi-quadratic denominator (BRIS). It is clear that Pij (x, y) satisfies the interpolating conditions
Pij (xr , ys ) = f (xr , ys ), r = i, i + 1, s = j, j + 1. Additionally, one can easily prove that
∂ Pij (x, yj ) ∂ Pij (x, yj+1 ) = ij (x), = i,j+1 (x), ∂y ∂y
(7)
and the piecewise function (6) on Dij is unique for the nine given data points {(xr , ys ), frs }, r = i, i + 1, i + 2; s = j, j + 1, j + 2. We remark that if we fix βij = βij∗ = 2, Pij (x, y) becomes the case in [15]. If we furthermore choose αij∗ = 1, βij∗ = 2, γij∗ = 1, αij = 1, βij = 2, γij = 1, Pij (x, y) becomes a bicubic Hermite spline. 3. Smoothing conditions For each yj , j = 1, 2, . . . , m + 1, from the definition of Pij (x, y), we can see that it is C1 continuous along the x-direction. Similarly, we can see from (7) that Pij (x, y) have a continuous first-order partial derivative Pij (x, y)/y in [x1 , xn ; y1 , yn ]. So it is sufficient for Pij (x, y) C1 in [x1 , xn ; y1 , ym ] if Pij (xi +, y)/x = Pi−1, j (xi −, y)/x holds. We then have the following theorem. Theorem 1. The interpolating spline Pij (x, y) is C1 in the whole interpolating region D if the parameters α i,j = constant, β i,j = constant and γ i,j = constant for each j, j = 1, 2, . . . , m − 1 and all i = 1, 2, . . . , n − 1. Proof. Based on above discussion, we just only prove that
∂ Pij (xi +, y) ∂ Pi−1,j (xi −, y) = . ∂y ∂y We first compute that
(8)
dP ∗ (x) dP ∗ (x) ∂ Pij (x, y) dV dW 1 = . (1 − η)3 αij ij + η(1 − η)2 ij + η2 (1 − η) ij + γ 3 βij i,j+1 ∂x qij (y) dx dx dx dx
(9)
Since
Pir∗ (xi +) = ∗ir , r = j, j + 1, j + 2,
(10)
hence
Vij (xi +) = Wij (xi +) =
βij ∗ij + αij ∗i,j+1 (βij + γij )∗i,j+1 − γij
lj (∗ − ∗i,j+1 ) lj+1 i,j+2
Substituting Eqs. (10) and (11) into (9), we obtain that
(11)
370
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
∂ Pij (xi +, y) 1 (1 − η)3 αij ∗ij + η(1 − η)2 (βij ∗ij + αij ∗i,j+1 ) = ∂x αij (1 − η)2 + βij η(1 − η) + γij η2 +η2 (1 − η)
(βij + γij )∗i,j+1 − γij
lj lj+1
(∗i,j+2 − ∗i,j+1 ) + η3 γij ∗i,j+1
(12)
∗ Similarly, since Pi−1,r (xi −) = ∗ir , r = j, j + 1, j + 2, we can also compute that
∂ Pi−1,j (xi −, y) 1 = (1 − η)3 αi−1,j ∗ij ∂x αi−1,j (1 − η)2 + βi−1,j η(1 − η) + γi−1,j η2 +η(1 − η)2 (βi−1,j ∗i−1,j + αi−1,j ∗i,j+1 ) +η2 (1 − η)
(βi−1,j + γi−1,j )∗i,j+1 − γi−1,j
lj lj+1
(∗i,j+2 − ∗i,j+1 )
+η3 γi−1,j ∗i,j+1 .
(13)
By comparing (12) with (13), we can yield α i−1, j = α ij , β i−1, j = β ij and γ i−1,j = γ ij . Therefore Eq. (8) holds. 4. Bounded property and error analysis We rewrite Pij∗ (x) as follows
Pij∗ (x) = ω0 (θ , αij∗ , βij∗ , γij∗ )fij + ω1 (θ , αij∗ , βij∗ , γij∗ )fi+1,j + ω2 (θ , αij∗ , βij∗ , γij∗ )fi+2,j
(14)
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 , α (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(15)
αij∗ θ (1 − θ )2 + βij∗ + γij∗ 1 + hhi+1i θ 2 (1 − θ ) + γij∗ θ 3 ω1 (θ , αij∗ , βij∗ , γij∗ ) = , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(16)
where
ω0 (θ , αij∗ , βij∗ , γij∗ ) =
∗ ij
ω2 (θ , αij∗ , βij∗ , γij∗ ) = −
γij∗ hhi+1i θ 2 (1 − θ ) , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(17)
ωr (θ , αij∗ , βij∗ , γij∗ ), r = 0, 1, 2 are called basis functions of Pij∗ (x). It is easy to compute that 2
ωr (θ , αij∗ , βij∗ , γij∗ ) = 1.
(18)
r=0
For any positive αij∗ , βij∗ , γij∗ , it is obvious that
ω0 (θ , αij∗ , βij∗ , γij∗ ) > 0, ω1 (θ , αij∗ , βij∗ , γij∗ ) > 0, ω2 (θ , αij∗ , βij∗ , γij∗ ) < 0. Therefore, from (15), we have
ω0 (θ , αij∗ , βij∗ , γij∗ ) = ω0 (θ , αij∗ , βij∗ , γij∗ ) ≤
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 = 1 − θ ≤ 1. αij∗ (1 − θ )2 + βij∗ θ (1 − θ )
(19)
From (16), we have
ω1 (θ , αij∗ , βij∗ , γij∗ ) =
ω1 (θ , αij∗ , βij∗ , γij∗ )
γij∗ 1 + hhi+1i θ 2 (1 − θ ) =θ+ ∗ αij (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 hi ≤ θ + 1+ (1 − θ ) hi+1 hi ≤ 1+ hi+1
(20)
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
371
From (17), we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) =
γij∗ hhi+1i θ 2 (1 − θ ) hi h ≤ (1 − θ ) ≤ i . hi+1 αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 hi+1
(21)
Similarly, we rewrite Pij (x, y) as follows: ∗ ∗ Pij (x, y) = ω0 (η, αij , βij , γij )Pij∗ (x) + ω1 (η, αij , βij , γij )Pi,j+1 (x) + ω2 (η, αij , βij , γij )Pi,j+2 (x),
(22)
where
ω0 (η, αij , βij , γij ) =
αij (1 − η)3 + βij η(1 − η)2 , αij (1 − η)2 + βij η(1 − η) + γij η2
lj αij η(1 − η)2 + βij + γij∗ 1 + lj+1 η2 (1 − η) + γij η3 ω1 (η, αij , βij , γij ) = , αij (1 − η)2 + βij η(1 − η) + γij η2
ω2 (η, αij , βij , γij ) = −
lj γij lj+1 η2 (1 − η) . αij (1 − η)2 + βij η(1 − η) + γij η2
Similarly to Eqs. (19)–(21), it follows that
|ω0 (η, αij , βij , γij )| ≤ 1, |ω1 (η, αij , βij , γij )| ≤ 1 +
(23)
lj , lj+1
(24)
lj . lj+1
|ω2 (η, αij , βij , γij )| ≤
(25)
ωr (θ , αij∗ , βij∗ , γij )∗ ωs (η, αij , βij , γij ), r, s = 0, 1, 2 are called basis functions of BRIS Pij (x, y). It is easy to compute from (18) that 2 2
∗ ∗ ∗ ωs (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s ) = 1.
(26)
s=0 r=0
Denote M = max |fi+r,j+s |, then we present the bounded theorem of BRIS Pij (x, y). r,s=0,1,2
Theorem 2. For any positive shape parameters αij∗ , βij∗ , γij∗ , αij , βij and γ ij , Pij (x, y) are bounded, namely
|Pij (x, y)| ≤ 4M 1 +
hi hi+1
1+
lj
lj+1
.
(27)
Proof. From (14) and (22), we have
Pij (x, y) =
2 2
∗ ∗ ∗ ω0 (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s )fi+r,j+s .
s=0 r=0
From Eqs. (19)–(21) and (23)–(25), it follows that
|Pij (x, y)| ≤ M |ω0 (η, αij , βij , γij )|
2
ωr (θ , αij∗ , βij∗ , γij∗ ) r=0
+M |ω1 (η, αij , βij , γij )|
2
∗ ∗ ∗ , βi,j+1 , γi,j+1 ) ωr (θ , αi,j+1 r=0
+M |ω2 (η, αij , βij , γij )|
hi ≤ 4M 1 + hi+1
2
∗ ∗ ∗ , βi,j+2 , γi,j+2 ) ωr (θ , αi,j+2 r=0
lj 1+ . lj+1
Let (x, y) Dij and denote
∂ = max ∂ , ∂ = max ∂ , ∂ x (x,y)∈Dij ∂ x ∂ y (x,y)∈Dij ∂ y hˆ i = max(hi , hi+1 ), lˆj = max(lj , lj+1 ).
We then have the following error formula.
(28)
372
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
Theorem 3. Suppose f(x, y) has a continuous first-order partial derivative, then the following error formula of BRIS Pij (x, y) holds
|f (x, y) − Pij (x, y)| ≤ 8 1 +
hi hi+1
1+
lj lj+1
∂ + lˆj ∂ hˆ i ∂x ∂y
(29)
Proof. The Taylor formula of f(x, y) at point (xi+r , yj+s ), r, s = 0, 1, 2 is that
f (x, y) = f (xi+r , yj+s ) +
(x − xi+r )
∂ ∂ + (y − yj+s ) f (μr , νs ), (x, y) ∈ Dij . ∂x ∂y
where μr is between x and xi+r , ν s is between y and yj+s . It follows that
∂ ∂ + l , r, s = 0, 1 ∂x j ∂y ∂ ∂ |f (x, y) − f (xi+2 , yj+2 )| ≤ (hi + hi+1 ) + (lj + lj+1 ) ∂x ∂y
|f (x, y) − f (xi+r , yj+s )| ≤ hi
It implies that
∂ ∂ |f (x, y) − f (xi+r , yj+s )| ≤ 2 hˆ i + lˆj , r, s = 0, 1, 2. ∂x ∂y From (26) and (28), we have
f (x, y) − Pij (x, y) =
2 2
∗ ∗ ∗ ω0 (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s )(f (x, y) − fi+r,j+s )
s=0 r=0
then
2 2 ∂ ∂
∗ ∗ ∗ , βi,j+s , γi,j+s ) . |f (x, y) − Pij (x, y)| ≤ 2 hˆ i + lˆj |ω0 (η, αij , βij , γij )| ωr (θ , αi,j+s ∂x ∂ y s=0 r=0
From (27), (29) holds. 5. Matrix expression and some properties of basis In this section, we consider the equally spaced knots case, namely, hi = h, lj = l, i = 1, 2, . . . , n, j = 1, 2, . . . , m. Assume αij∗ , βij∗ and γij∗ , i = 1, 2, . . . , n − 1, j = 1, 2, . . . , m + 1 are constant respectively. 5.1. Matrix expression and the symmetry of BRIS In this case the basis of Pij∗ (x) satisfies that ∗ ∗ ∗ ∗ ∗ ∗ ωr (θ , αij∗ , βij∗ , γij∗ ) = ωr (θ , αi,j+1 , βi,j+1 , γi,j+1 ) = ωr (θ , αi,j+2 , βi,j+2 , γi,j+2 ), r = 0, 1, 2,
where
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 , α (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 α θ (1 − θ )2 + βij∗ + 2γij∗ θ 2 (1 − θ ) + γij∗ θ 3 ω1 (θ , αij∗ , βij∗ , γij∗ ) = , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 γij∗ θ 2 (1 − θ ) ω2 (θ , αij∗ , βij∗ , γij∗ ) = − ∗ . αij (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 ω0 (θ , αij∗ , βij∗ , γij∗ ) =
Hence we have that ∗ Pi,j+s (x) =
Denote
(30)
∗ ij ∗ ij
⎡
(31) (32)
fi,j+s
⎤
ω0 (θ , αij∗ , βij∗ , γij∗ ), ω1 (θ , αij∗ , βij∗ , γij∗ ), ω2 (θ , αij∗ , βij∗ , γij∗ ) ⎣ fi+1,j+s ⎦ , s = 0, 1, 2.
fi+2,j+s
= ω0 (θ , αij∗ , βij∗ , γij∗ ) ω1 (θ , αij∗ , βij∗ , γij∗ ) ω2 (θ , αij∗ , βij∗ , γij∗ ) , T Fs = fi,j+s , fi,j+s , fi,j+s , then ∗ Pi,j+s (x) = Fs .
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
373
Denote
F = F0 F1 F2 ,
∗ ∗ (x) Pi,j+2 (x) , A = Pij∗ (x) Pi,j+1
then A = F. Denote
H=
ω0 (η, αij , βij , γij ), ω1 (η, αij , βij , γij ), ω2 (η, αij , βij , γij ) ,
we can rewrite Pij (x, y) as following two forms
Pij (x, y) = HAT = HF T T
(33)
= AHT = FHT .
(34)
Eq. (33) is the matrix expression of Pij (x, y). It means from Eqs. (33) and (34) that the expression of BRIS Pij (x, y) are the same whatever direction Pij (x, y) begins with first. This property is called the symmetry of BRIS Pij (x, y). 5.2. Some properties of basis Except for properties Eqs. (18) and (26), the basis of BRIS have some other interesting properties. Denote
ar =
0
1
ωr (θ , αij∗ , βij∗ , γij∗ )dθ , r = 0, 1, 2,
which are called the integral weight coefficients of Pij∗ (x) in (1). From (18), it is clear that
2
r=0
ar = 1.
Property 1. If fij = fi+1, j = fi+2,j = 1, there holds that
xi+1 xi
Pij∗ (x)dx = hi
2
ar fi+r,j = hi .
r=0
From Eqs. (30) to (32), we obtain that
a0 − a2 =
0
a1 + 2a2 =
1
0
1
1 2
(1 − θ )dθ = , 1 2
θ dθ = .
From (32), we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) = −
γij∗ θ 2 (1 − θ ) > −(1 − θ ), α (1 − θ ) + βij∗ θ (1 − θ ) + γij∗ θ 2
−
∗ ij
1 0
2
1 6
θ (1 − θ )dθ = − .
It leads to the following property. Property 2. In any subinterval [xi , xi+1 ], no matter what positive number the shape parameters αij∗ , βij∗ and γij∗ are, the integral weight coefficients ar ’s satisfy
1 1 < a0 < , 3 2
1 5 < a1 < , 2 6
−
1 < a2 < 0. 6
(35)
If αij∗ , βij∗ > γij∗ > 0, we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) > −
γij∗ θ 2 (1 − θ ) θ 2 (1 − θ ) = − , ∗ ∗ γ (1 − θ )2 + γij θ (1 − θ ) + γij θ 2 1 − θ + θ2
−
∗ ij
1 0
1 θ 2 (1 − θ ) dθ = − 2 1 − θ + θ2
√ 3π . 9
Property 3. In any subinterval [xi , xi+1 ], If αij∗ , βij∗ > γij∗ > 0, the integral weight coefficients ar ’s satisfy
√ 1 3π < a0 < , 1− 9 2
√ 1 2 3π 1 < a1 < − + , 2 2 9
1 − 2
√ 3π < a2 < 0. 9
(36)
374
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377 Table 2 Data points. f(xi , yj )
x1 x2 x3 x4
=1 =2 =3 =4
y1 = 1
y2 = 2
y3 = 3
y4 = 4
4 2 3 1
2 3 2 4
3 2 4 2
1 4 2 3
Let D˜ = [0, 1; 0, 1]. Denote
ωrs (θ , αij∗ , βij∗ , γij∗ ; η, αij , βij , γij ) = ωr (θ , αij∗ , βij∗ , γij∗ )ωs (η, αij , βij , γij ),
ars =
D˜
ωrs (θ , αij∗ , βij∗ , γij∗ ; η, αij , βij , γij )dθ dη
ars are called the integral weight coefficients of BRIS Pij (x, y). From Eq. (26), one can yield that
2
r=0
2
s=0
ars = 1.
Property 4. If fi+r,j+s = 1, r, s = 0, 1, 2, then
D
Pij (x, y)dxdy = hi lj
2 2
ars fi+r,j+s = hi lj .
r=0 s=0
From Eq. (35) in Property 2, we obtain the following property. Property 5. In any subinterval Dij , no matter what positive number the parameters αij∗ , βij∗ , γij∗ , αij , βij , γij are, the integral weight coefficients ars ’s satisfy
1 1 1 5 1 1 < a00 < < a01 < − < a02 < 9 4 6 12 12 4 1 5 5 25 1 < a10 < < a11 < − < a12 < 0 6 12 4 36 36 5 1 1 1 < a20 < − < a21 < 0 0 < a22 < − 12 4 36 36 From Eq. (36), we have that Property 6. In any subinterval Dij , if αij∗ , βij∗ > γij∗ > 0, αij , βij > γij > 0, ars ’s satisfy
2 √ √ √ √ 1 1 1 3π 3π 3π 1 3π − < a01 < − + − < a02 < 0 < a00 < 9 4 2 18 4 9 4 18 2 √ √ √ √ 2 3π 2π 2 1 1 1 3π 3π 1 3π 1 − < a10 < − + < a11 < − + − < a12 < 0 − + 2 18 4 9 4 2 9 4 6 27 2 √ √ √ 2π 2 1 1 3π 3π 3π 1 − < a20 < 0 − + − < a21 < 0 0 < a12 < − 4 18 4 6 27 2 9
1−
6. Numerical examples Let D = [1, 4; 1, 4], and let the knots be equally spaced, the set of data points are given in Table 2, then BRIS Pij (x, y) can be constructed in [1, 3; 1, 3] for the given positive parameters αij∗ , βij∗ , γij∗ , αij , βij and γ ij . The corresponding graphes are shown in Fig. 5–8. In Fig. 5 each patch tilts heavily nearby two adjacent wedges. This phenomenon is similar to that in Fig. 1. In Fig. 6 and 7 the interpolating vertexes can be identified easily, and each patch is relatively plane because αij∗ βij∗ , γij∗ ; αij βij , γij in Fig. 6 and βij∗ αij∗ , γij∗ ; βij αij , γij in Fig. 7. The surface in Fig. 8 is a bicubic Hermite interpolating spline. In the proposed method we generate four shape control parameters in [15] to six parameters case. In [15] when γ ij α ij or γij∗ αij∗ , the interpolating surface tilts heavily. After adding two parameters βij∗ , βij and making βij∗ be close to or greater than γij∗ , β ij be close to or greater than γ ij , we can avoid the above inharmony of interpolating surface. A compared effect is shown in Fig. 9. In Fig. 9(a), according to [15], γij∗ = γij = 100 αij∗ = αij∗ = 1 while βij∗ = βij = 2, the corresponding surface tilts heavily. In Fig. 9(b), when we adjust the value of the parameters βij∗ and β ij to 200, the inharmonic phenomenon vanishes.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
6 3.0
z 4 2 2.5
0 2 1.0
2.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 5. BRIS with αij∗ = 1, βij∗ = 1, γij∗ = 10, αij = 1, βij = 1, γij = 100.
4 3.0
z
3 2.5 2 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 6. BRIS with αij∗ = 20, βij∗ = 1, γij∗ = 1, αij = 20, βij = 1, γij = 1.
4.0 3.5 z
3.0
3.0 2.5
2.5
2.0 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 7. BRIS with αij∗ = 1, βij∗ = 20, γij∗ = 1, αij = 1, βij = 20, γij = 1.
375
376
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
4 z
3.0 3 2.5
2 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 8. BRIS with αij∗ = 1, βij∗ = 2, γij∗ = 1, αij = 1, βij = 2, γij = 1.
4
6 3.0
z 4 2 2.5
0 2.0
1.0 1.5
y
z
3.0 3 2.5
2 2.0
1.0
y
1.5 1.5
2.0 x
1.5
2.0 x
2.5
2.5 3.0
1.0
∗ = (a) The case with αij = α∗ij = 1, βij = βij ∗ = 100 in [15] 2, γij = γij
3.0
1.0
∗ = (b) BRIS with αij = α∗ij = 1, βij = βij ∗ = 100 200, γij = γij
Fig. 9. Comparison the proposed six shape parameters case with the four parameters in [15] (a) The case with αij = αij∗ = 1, βij = βij∗ = 2, γij = γij∗ = 100 in [15] (b) BRIS with αij = αij∗ = 1, βij = βij∗ = 200, γij = γij∗ = 100.
Acknowledgments This work was supported by the NNSF of China (11472063) and the Provincial Natural Science Research Program of Higher Education Institutions of Anhui Province (KJ2013A194, KJ2013Z230). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
M. Abbas, A.A. Majid, J.M. Ali, Positivity-preserving C2 rational cubic spline interpolation, Sci. Asia 39 (2013) 208–213. M. Abbas, A.A. Majid, J.M. Ali, Positivity-preserving rational bi-cubic spline interpolation for 3D positive data, Appl. Math. Comput. 234 (2014) 460–476. F. Bao, Q. Sun, Q. Duan, Point control of the interpolating curve with a rational cubic spline, J. Vis. Commun. Image R. 20 (2009) 275–280. F. Bao, Q. Sun, J. Pan, Q. Duan, A blending interpolator with value control and minimal strain energy, Comput. Graph. 34 (2010) 119–124. P.E. Bézier, The Mathematical Basis of The UNISURF CAD System, Butterworth, London, 1986. C. de Boor, B-form basics, in: G. Farin (Ed.), Geometric Modeling, SIAM, Philadelphia, 1987, pp. 131–148. C.K. Chui, Multivariate Splines, SIAM, 1988. R. Delbourgo, Shape preserving interpolation to convex data by rational functions with quadratic numerator and linear denominator, IMA J. Numer. Anal. 9 (1989) 23–136. J. Deng, F. Chen, Y. Feng, Dimensions of spling spaces over t-meshes, J. Comput. Appl. Math. 194 (2006) 267–283. J. Deng, F. Chen, X. Li, Graph. Models 74 (2008) 76–86. P. Dierck, B. Tytgat, Generating the Bézier points of BETA-spline curve, Comput. Aided Geom. Des. 6 (1989) 279–291. Q. Duan, K. Djidjeli, W.G. Price, E.H. Twizell, The approximation properties of some rational cubic splines, Int. J. Comput. Math. 72 (1999) 155–166.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
377
Q. Duan, L. Wang, E.H. Twizell, A new bivariate rational interpolation based on function values, Inf. Sci. 166 (2004) 181–191. Q. Duan, L. Wang, E.H. Twizell, A new weighted rational cubic interpolation and its approximation, Appl. Math. Comput. 168 (2005) 990–1003. Q. Duan, H. Zhang, A. Liu, H. Li, A bivariate rational interpolation with a bi-quadratic denominator, J. Comput. Appl. Math. 195 (2006) 24–33. Q. Duan, Y. Zhang, E.H. Twizell, A bivariate rational interpolation and the properties, Appl. Math. Comput. 179 (2006) 190–199. Q. Duan, S. Li, F. Bao, E.H. Twizell, Hermite interpolation by piecewise rational surface, Appl. Math. Comput. 198 (2008) 59–72. Q. Duan, F. Bao, S. Du, E.H. Twizell, Local control of interpolating rational cubic spline curves, Comput. Aided Des. 41 (2009) 825–829. G. Farin, Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide, fifth ed., Morgan Kaufman, Menlo Park, 2002. J.A. Gregory, M. Sarfraz, P.K. Yuen, Interactive curve design using C2 rational splines, Comput. Graph. 18 (1994) 153–159. X. Han, Convexity preserving piecewise rational quartic interpolation, SIAM J. Numer. Anal. 46 (2008) 920–929. K. Konno, H. Chiyokura, An approach of designing and controlling free-form surfaces by using NURBS boundary Gregory patches, Comput. Aided Geom. Des. 13 (1996) 825–849. M.Z. Hussain, M. Sarfraz, Positivity-preserving interpolation of positive data by rational cubics, J. Comput. Appl. Math. 218 (2008) 446–458. R. Müler, Universal parametrization and interpolation on cubic surfaces, Comput. Aided Geom. Des. 19 (2002) 479–502. L. Piegl, On NURBS: a survey, IEEE Comput. Graph. Appl. 11 (1991) 55–71. M. Sarfraz, A C2 rational cubic spline which has linear denominator and shape control, Ann. Univ. Sci. Budapest 37 (1994) 53–62. T.W. Sederberg, J.M. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCS, ACM Trans. Graph. 22 (2003) 161–172. J. Tan, S. Tang, Composite schemes for multivariate blending rational interpolation, J. Comput. Appl. Math. 144 (2002) 263–275. R. Wang, Multivariate Spline Functions and Their Applications, Kluwer Academic Publishers, Beijing/New York/Dordrecht/Boston/London, 2001. Y. Zhang, Q. Duan, E.H. Twizell, Convexity control of a bivariate rational interpolating spline surfaces, Comput. Graph. 31 (2007) 679–687.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A bivariate rational cubic interpolating spline with biquadratic denominator Youtian Tao a,b,c,∗, Dongyin Wang a a
Department of Mathematics, Chaohu College, Chaohu 238000, China School of Mathematics Sciences, USTC, Hefei 230026, China c Anhui Fuhuang Steel Structure, Chaohu 238076, China b
a r t i c l e
i n f o
Keywords: Bivariate rational interpolating spline Shape parameter Bounded property Error estimate Symmetry
a b s t r a c t A bivariate rational bicubic interpolating spline (BRIS) with biquadratic denominator and six shape parameters is constructed in a rectangle domain. The C1 continuous condition of BRIS discussed. BRIS is proved to be bounded and its error is estimated. In the case of the equally spaced knots, the matrix expression and symmetry of BRIS are presented. Some properties of the basis of BRIS are given. In the end, a numerical example is given to illustrat the effect of the shape parameters on the shape of BRIS surface. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The construction method of curves and surfaces and their mathematical description is a key issue in CAGD. There are many ways [5–7,9–11,19,22,24,25,27–29] to deal with this problem, such as the polynomial spline method, Bézier spline method and NURBS method. However, most of the polynomial spline methods are the interpolating cases, and their local shapes can not be modified for the interpolating surfaces while interpolating data is unchanged. NURBS and Bézier methods are the noninterpolating cases, that is to say, the constructed curve and surface do not fit with the given data and the given points are the control points. Therefore, when we construct the interpolating functions required for CAGD, we should consider the following cases: 1) The expressions of interpolating functions are simple and explicit. 2) The parameters of constructed curves and surfaces can be modified without changing the given data. Recently, there has been many works [3,4,8,12,18,20,21,26] about univariate rational spline interpolation with parameters. Some univariate rational interpolating splines are generalized to bivariate rational interpolating splines which expressions with parameters are simple and explicit. Some interesting results are presented. In [2,13–17,23,30], the authors constructed several bivariate interpolating splines over rectangular mesh, derived some properties such as the sufficient conditions of downconstrained and up-constrained for the shape control, the matrix expression, bounded property, stability, convexity control, the preserving positivity. In this paper, motivated by [2,15], we generate the bivariate rational interpolating spline with four shape parameters in [15] to the case with six shape parameters (BRIS). We present the definition of BRIS and discuss the effect of the parameter choice on the curve shape in Section 2. We give the sufficient condition that BRIS is C1 continuous in Section 3. We prove the bounded property and analyse the error of BRIS in Section 4. In the case of the equally spaced knots, we present the matrix expression of
∗
Corresponding author at: Chaohu College, Department of Mathematics, Jinzhai Rd., Chaohu 238000, China. Tel: +8615156520160. E-mail address: [email protected] (Y. Tao).
http://dx.doi.org/10.1016/j.amc.2015.04.100 0096-3003/© 2015 Elsevier Inc. All rights reserved.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
367
Table 1 Data points from [1]. i
1
2
3
4
5
6
7
8
9
xi fij
0 2
0.25 0.6
0.5 0.1
1 0.13
1.5 1
2 0.5
2.5 1.1
3 0.25
4 0.2
BRIS, prove that it is symmetry, and discuss some properties of the basis of BRIS in Section 5. We give a numerical example of BRIS to illustrat the effect of the shape parameters on the shape of BRIS surface in Section 6. 2. Rational interpolating spline Let D be the rectangular domain [a, b, c, d], {(xi , yi ), fij , i = 1, 2, . . . , n + 1; j = 1, 2, . . . , m + 1} be a given set of data points, where a = x1 < x2 < < xn+1 = b; c = y1 < y2 < < ym+1 = d. Set
fij = f (xi , yj ), hi = xi+1 − xi , lj = yj+1 − yj . Denote Dij = [xi , xi+1 yj , yj+1 ], i = 1, 2, . . . , n, j = 1, 2, . . . , m. For any point (x, y) Dij , let
θ = (x − xi )/ hi , η = (y − yj )/ lj . For each yj , j = 1, 2, . . . , m + 1, we construct the x-direction interpolation curve as follows
Pij∗ (x) =
p∗ij (x) q∗ij (x)
, i = 1, 2, . . . , n − 1,
(1)
where
p∗ij (x) =
αij∗ fij (1 − θ )3 + Vij∗ θ (1 − θ )2 + Wij∗ θ 2 (1 − θ ) + γij∗ fi+1,j θ 3
q∗ij (x) =
αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
and
βij∗ fij + αij∗ fi+1,j = (βij∗ + γij∗ )fi+1,j − γij∗ hi ∗i+1
Vij∗ = Wij∗
with αij∗ , βij∗ , γij∗ positive, and ∗i =
fi+1,j −fij . hi
We can prove that
Pij∗ (xi ) = fij , Pij∗ (xi+1 ) = fi+1,j , Pij∗ (xi ) = ∗ij , Pij∗ (xi+1 ) = ∗i+1,j . The interpolation (1) is called the rational cubic interpolating spline. It is clear that the interpolation is local in the interval [xi , xi+1 ] and depends on the data at three points {(xr , yj ), frj }, r = i, i + 1, i + 2 and the shape parameters αij∗ , βij∗ , γij∗ . It is interesting to note that the interpolation (1) becomes a standard cubic Hermite spline with the values of the shape parameters αij∗ = 1, βij∗ = 2, γij∗ = 1. we now illustrate the mathematical and graphical effects of the shape parameters αij∗ , βij∗ , γij∗ on the shape of a curve. The three free parameters can be exploited properly to modify the shape of curve according to the designer’s choice. We rewrite (1) as follows:
Pij∗ (x) =
αij∗ (1 − θ )2 (fij (1 − θ ) + fi+1,j θ ) + βij∗ θ (1 − θ )(fij (1 − θ ) + fi+1,j θ ) + γij∗ θ 2 (fi+1,j − hi (1 − θ )i+1,j ) αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(2)
By the computation for Eq. (2) we have the following formulas:
lim P ∗ (x) = fij (1 − θ ) + fi+1,j θ αij∗ →∞ ij
(3)
lim Pij∗ (x) = fij (1 − θ ) + fi+1,j θ
(4)
βij∗ →∞
lim P ∗ (x) = fi+1,j − (1 − θ )hi i+1,j γij∗ →∞ ij
(5)
We can see from Eqs. (3) and (4) that the increase in the shape parameter αij∗ or βij∗ reduces the rational interpolating spline (2) to the straight line fij (1 − θ ) + fi+1, j θ , and see from Eq. (5) that the increase in γij∗ reduces (2) to another straight line fi+1, j − (1 − θ )hi i+1, j . According to data points in Table 1 from [1], by choosing different shape parameters, we can observe the corresponding shape changes of the interpolating curves in Figs. 1–4. Each piecewise in Fig. 1 slopes heavily at the left-hand side. Each piecewise in Fig. 2 and 3 incline to a line segment. The curve in Fig. 4 is a cubic Hermite interpolating spline. We remark that in order to plot the 8th interpolating spline curve, we add the point {x10 , f10 } = {4.5, 0.2}.
368
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
2.0
1.5
1.0
0.5
0.0
0.5
0
1
2
3
4
Fig. 1. The rational spline interpolation Eq. (1) with α = 1, β = 1, γ = 100. ∗ ij
∗ ij
∗ ij
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 2. The rational spline interpolation Eq. (1) with α = 10, β = 1, γ = 1. ∗ ij
∗ ij
∗ ij
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 3. The rational spline interpolation Eq. (1) with αij∗ = 1, βij∗ = 10, γij∗ = 1.
We now use the x-direction interpolation Pij∗ (x), i = 1, 2, . . . , n − 1; j = 1, 2, . . . , m + 1 to construct the bivariate rational interpolating function spline (BRIS) in D. For each pair (i, j), i = 1, 2, . . . , n − 1; j = 1, 2, . . . , m − 1, we define a bivariate interpolating spline Pij (x, y) in Dij as follows
Pij (x, y) =
pij (x, y) , qij (x, y)
pij (x, y) =
∗ αij Pij∗ (x)(1 − η)3 + Vij η(1 − η)2 + Wij η2 (1 − η) + γij η3 Pi,j+1 (x)
qij (x, y) =
αij (1 − η)2 + βij η(1 − η) + γij η2
where
(6)
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
369
2.0
1.5
1.0
0.5
0.0 0
1
2
3
4
Fig. 4. The rational spline interpolation Eq. (1) with αij∗ = 1, βij∗ = 2, γij∗ = 1.
and ∗ βij Pij∗ (x) + αij Pi,j+1 (x) ∗ Wij = (βij + γij )Pi,j+1 (x) − γij lj i,j+1 (x)
Vij =
∗ with α ij , β ij , γ ij positive, and ij (x) = (Pi,j+1 (x) − Pij∗ (x))/ lj . Pij (x, y) is called a bivariate bicubic rational interpolating spline with a bi-quadratic denominator (BRIS). It is clear that Pij (x, y) satisfies the interpolating conditions
Pij (xr , ys ) = f (xr , ys ), r = i, i + 1, s = j, j + 1. Additionally, one can easily prove that
∂ Pij (x, yj ) ∂ Pij (x, yj+1 ) = ij (x), = i,j+1 (x), ∂y ∂y
(7)
and the piecewise function (6) on Dij is unique for the nine given data points {(xr , ys ), frs }, r = i, i + 1, i + 2; s = j, j + 1, j + 2. We remark that if we fix βij = βij∗ = 2, Pij (x, y) becomes the case in [15]. If we furthermore choose αij∗ = 1, βij∗ = 2, γij∗ = 1, αij = 1, βij = 2, γij = 1, Pij (x, y) becomes a bicubic Hermite spline. 3. Smoothing conditions For each yj , j = 1, 2, . . . , m + 1, from the definition of Pij (x, y), we can see that it is C1 continuous along the x-direction. Similarly, we can see from (7) that Pij (x, y) have a continuous first-order partial derivative Pij (x, y)/y in [x1 , xn ; y1 , yn ]. So it is sufficient for Pij (x, y) C1 in [x1 , xn ; y1 , ym ] if Pij (xi +, y)/x = Pi−1, j (xi −, y)/x holds. We then have the following theorem. Theorem 1. The interpolating spline Pij (x, y) is C1 in the whole interpolating region D if the parameters α i,j = constant, β i,j = constant and γ i,j = constant for each j, j = 1, 2, . . . , m − 1 and all i = 1, 2, . . . , n − 1. Proof. Based on above discussion, we just only prove that
∂ Pij (xi +, y) ∂ Pi−1,j (xi −, y) = . ∂y ∂y We first compute that
(8)
dP ∗ (x) dP ∗ (x) ∂ Pij (x, y) dV dW 1 = . (1 − η)3 αij ij + η(1 − η)2 ij + η2 (1 − η) ij + γ 3 βij i,j+1 ∂x qij (y) dx dx dx dx
(9)
Since
Pir∗ (xi +) = ∗ir , r = j, j + 1, j + 2,
(10)
hence
Vij (xi +) = Wij (xi +) =
βij ∗ij + αij ∗i,j+1 (βij + γij )∗i,j+1 − γij
lj (∗ − ∗i,j+1 ) lj+1 i,j+2
Substituting Eqs. (10) and (11) into (9), we obtain that
(11)
370
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
∂ Pij (xi +, y) 1 (1 − η)3 αij ∗ij + η(1 − η)2 (βij ∗ij + αij ∗i,j+1 ) = ∂x αij (1 − η)2 + βij η(1 − η) + γij η2 +η2 (1 − η)
(βij + γij )∗i,j+1 − γij
lj lj+1
(∗i,j+2 − ∗i,j+1 ) + η3 γij ∗i,j+1
(12)
∗ Similarly, since Pi−1,r (xi −) = ∗ir , r = j, j + 1, j + 2, we can also compute that
∂ Pi−1,j (xi −, y) 1 = (1 − η)3 αi−1,j ∗ij ∂x αi−1,j (1 − η)2 + βi−1,j η(1 − η) + γi−1,j η2 +η(1 − η)2 (βi−1,j ∗i−1,j + αi−1,j ∗i,j+1 ) +η2 (1 − η)
(βi−1,j + γi−1,j )∗i,j+1 − γi−1,j
lj lj+1
(∗i,j+2 − ∗i,j+1 )
+η3 γi−1,j ∗i,j+1 .
(13)
By comparing (12) with (13), we can yield α i−1, j = α ij , β i−1, j = β ij and γ i−1,j = γ ij . Therefore Eq. (8) holds. 4. Bounded property and error analysis We rewrite Pij∗ (x) as follows
Pij∗ (x) = ω0 (θ , αij∗ , βij∗ , γij∗ )fij + ω1 (θ , αij∗ , βij∗ , γij∗ )fi+1,j + ω2 (θ , αij∗ , βij∗ , γij∗ )fi+2,j
(14)
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 , α (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(15)
αij∗ θ (1 − θ )2 + βij∗ + γij∗ 1 + hhi+1i θ 2 (1 − θ ) + γij∗ θ 3 ω1 (θ , αij∗ , βij∗ , γij∗ ) = , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(16)
where
ω0 (θ , αij∗ , βij∗ , γij∗ ) =
∗ ij
ω2 (θ , αij∗ , βij∗ , γij∗ ) = −
γij∗ hhi+1i θ 2 (1 − θ ) , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2
(17)
ωr (θ , αij∗ , βij∗ , γij∗ ), r = 0, 1, 2 are called basis functions of Pij∗ (x). It is easy to compute that 2
ωr (θ , αij∗ , βij∗ , γij∗ ) = 1.
(18)
r=0
For any positive αij∗ , βij∗ , γij∗ , it is obvious that
ω0 (θ , αij∗ , βij∗ , γij∗ ) > 0, ω1 (θ , αij∗ , βij∗ , γij∗ ) > 0, ω2 (θ , αij∗ , βij∗ , γij∗ ) < 0. Therefore, from (15), we have
ω0 (θ , αij∗ , βij∗ , γij∗ ) = ω0 (θ , αij∗ , βij∗ , γij∗ ) ≤
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 = 1 − θ ≤ 1. αij∗ (1 − θ )2 + βij∗ θ (1 − θ )
(19)
From (16), we have
ω1 (θ , αij∗ , βij∗ , γij∗ ) =
ω1 (θ , αij∗ , βij∗ , γij∗ )
γij∗ 1 + hhi+1i θ 2 (1 − θ ) =θ+ ∗ αij (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 hi ≤ θ + 1+ (1 − θ ) hi+1 hi ≤ 1+ hi+1
(20)
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
371
From (17), we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) =
γij∗ hhi+1i θ 2 (1 − θ ) hi h ≤ (1 − θ ) ≤ i . hi+1 αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 hi+1
(21)
Similarly, we rewrite Pij (x, y) as follows: ∗ ∗ Pij (x, y) = ω0 (η, αij , βij , γij )Pij∗ (x) + ω1 (η, αij , βij , γij )Pi,j+1 (x) + ω2 (η, αij , βij , γij )Pi,j+2 (x),
(22)
where
ω0 (η, αij , βij , γij ) =
αij (1 − η)3 + βij η(1 − η)2 , αij (1 − η)2 + βij η(1 − η) + γij η2
lj αij η(1 − η)2 + βij + γij∗ 1 + lj+1 η2 (1 − η) + γij η3 ω1 (η, αij , βij , γij ) = , αij (1 − η)2 + βij η(1 − η) + γij η2
ω2 (η, αij , βij , γij ) = −
lj γij lj+1 η2 (1 − η) . αij (1 − η)2 + βij η(1 − η) + γij η2
Similarly to Eqs. (19)–(21), it follows that
|ω0 (η, αij , βij , γij )| ≤ 1, |ω1 (η, αij , βij , γij )| ≤ 1 +
(23)
lj , lj+1
(24)
lj . lj+1
|ω2 (η, αij , βij , γij )| ≤
(25)
ωr (θ , αij∗ , βij∗ , γij )∗ ωs (η, αij , βij , γij ), r, s = 0, 1, 2 are called basis functions of BRIS Pij (x, y). It is easy to compute from (18) that 2 2
∗ ∗ ∗ ωs (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s ) = 1.
(26)
s=0 r=0
Denote M = max |fi+r,j+s |, then we present the bounded theorem of BRIS Pij (x, y). r,s=0,1,2
Theorem 2. For any positive shape parameters αij∗ , βij∗ , γij∗ , αij , βij and γ ij , Pij (x, y) are bounded, namely
|Pij (x, y)| ≤ 4M 1 +
hi hi+1
1+
lj
lj+1
.
(27)
Proof. From (14) and (22), we have
Pij (x, y) =
2 2
∗ ∗ ∗ ω0 (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s )fi+r,j+s .
s=0 r=0
From Eqs. (19)–(21) and (23)–(25), it follows that
|Pij (x, y)| ≤ M |ω0 (η, αij , βij , γij )|
2
ωr (θ , αij∗ , βij∗ , γij∗ ) r=0
+M |ω1 (η, αij , βij , γij )|
2
∗ ∗ ∗ , βi,j+1 , γi,j+1 ) ωr (θ , αi,j+1 r=0
+M |ω2 (η, αij , βij , γij )|
hi ≤ 4M 1 + hi+1
2
∗ ∗ ∗ , βi,j+2 , γi,j+2 ) ωr (θ , αi,j+2 r=0
lj 1+ . lj+1
Let (x, y) Dij and denote
∂ = max ∂ , ∂ = max ∂ , ∂ x (x,y)∈Dij ∂ x ∂ y (x,y)∈Dij ∂ y hˆ i = max(hi , hi+1 ), lˆj = max(lj , lj+1 ).
We then have the following error formula.
(28)
372
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
Theorem 3. Suppose f(x, y) has a continuous first-order partial derivative, then the following error formula of BRIS Pij (x, y) holds
|f (x, y) − Pij (x, y)| ≤ 8 1 +
hi hi+1
1+
lj lj+1
∂ + lˆj ∂ hˆ i ∂x ∂y
(29)
Proof. The Taylor formula of f(x, y) at point (xi+r , yj+s ), r, s = 0, 1, 2 is that
f (x, y) = f (xi+r , yj+s ) +
(x − xi+r )
∂ ∂ + (y − yj+s ) f (μr , νs ), (x, y) ∈ Dij . ∂x ∂y
where μr is between x and xi+r , ν s is between y and yj+s . It follows that
∂ ∂ + l , r, s = 0, 1 ∂x j ∂y ∂ ∂ |f (x, y) − f (xi+2 , yj+2 )| ≤ (hi + hi+1 ) + (lj + lj+1 ) ∂x ∂y
|f (x, y) − f (xi+r , yj+s )| ≤ hi
It implies that
∂ ∂ |f (x, y) − f (xi+r , yj+s )| ≤ 2 hˆ i + lˆj , r, s = 0, 1, 2. ∂x ∂y From (26) and (28), we have
f (x, y) − Pij (x, y) =
2 2
∗ ∗ ∗ ω0 (η, αij , βij , γij )ωr (θ , αi,j+s , βi,j+s , γi,j+s )(f (x, y) − fi+r,j+s )
s=0 r=0
then
2 2 ∂ ∂
∗ ∗ ∗ , βi,j+s , γi,j+s ) . |f (x, y) − Pij (x, y)| ≤ 2 hˆ i + lˆj |ω0 (η, αij , βij , γij )| ωr (θ , αi,j+s ∂x ∂ y s=0 r=0
From (27), (29) holds. 5. Matrix expression and some properties of basis In this section, we consider the equally spaced knots case, namely, hi = h, lj = l, i = 1, 2, . . . , n, j = 1, 2, . . . , m. Assume αij∗ , βij∗ and γij∗ , i = 1, 2, . . . , n − 1, j = 1, 2, . . . , m + 1 are constant respectively. 5.1. Matrix expression and the symmetry of BRIS In this case the basis of Pij∗ (x) satisfies that ∗ ∗ ∗ ∗ ∗ ∗ ωr (θ , αij∗ , βij∗ , γij∗ ) = ωr (θ , αi,j+1 , βi,j+1 , γi,j+1 ) = ωr (θ , αi,j+2 , βi,j+2 , γi,j+2 ), r = 0, 1, 2,
where
αij∗ (1 − θ )3 + βij∗ θ (1 − θ )2 , α (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 α θ (1 − θ )2 + βij∗ + 2γij∗ θ 2 (1 − θ ) + γij∗ θ 3 ω1 (θ , αij∗ , βij∗ , γij∗ ) = , αij∗ (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 γij∗ θ 2 (1 − θ ) ω2 (θ , αij∗ , βij∗ , γij∗ ) = − ∗ . αij (1 − θ )2 + βij∗ θ (1 − θ ) + γij∗ θ 2 ω0 (θ , αij∗ , βij∗ , γij∗ ) =
Hence we have that ∗ Pi,j+s (x) =
Denote
(30)
∗ ij ∗ ij
⎡
(31) (32)
fi,j+s
⎤
ω0 (θ , αij∗ , βij∗ , γij∗ ), ω1 (θ , αij∗ , βij∗ , γij∗ ), ω2 (θ , αij∗ , βij∗ , γij∗ ) ⎣ fi+1,j+s ⎦ , s = 0, 1, 2.
fi+2,j+s
= ω0 (θ , αij∗ , βij∗ , γij∗ ) ω1 (θ , αij∗ , βij∗ , γij∗ ) ω2 (θ , αij∗ , βij∗ , γij∗ ) , T Fs = fi,j+s , fi,j+s , fi,j+s , then ∗ Pi,j+s (x) = Fs .
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
373
Denote
F = F0 F1 F2 ,
∗ ∗ (x) Pi,j+2 (x) , A = Pij∗ (x) Pi,j+1
then A = F. Denote
H=
ω0 (η, αij , βij , γij ), ω1 (η, αij , βij , γij ), ω2 (η, αij , βij , γij ) ,
we can rewrite Pij (x, y) as following two forms
Pij (x, y) = HAT = HF T T
(33)
= AHT = FHT .
(34)
Eq. (33) is the matrix expression of Pij (x, y). It means from Eqs. (33) and (34) that the expression of BRIS Pij (x, y) are the same whatever direction Pij (x, y) begins with first. This property is called the symmetry of BRIS Pij (x, y). 5.2. Some properties of basis Except for properties Eqs. (18) and (26), the basis of BRIS have some other interesting properties. Denote
ar =
0
1
ωr (θ , αij∗ , βij∗ , γij∗ )dθ , r = 0, 1, 2,
which are called the integral weight coefficients of Pij∗ (x) in (1). From (18), it is clear that
2
r=0
ar = 1.
Property 1. If fij = fi+1, j = fi+2,j = 1, there holds that
xi+1 xi
Pij∗ (x)dx = hi
2
ar fi+r,j = hi .
r=0
From Eqs. (30) to (32), we obtain that
a0 − a2 =
0
a1 + 2a2 =
1
0
1
1 2
(1 − θ )dθ = , 1 2
θ dθ = .
From (32), we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) = −
γij∗ θ 2 (1 − θ ) > −(1 − θ ), α (1 − θ ) + βij∗ θ (1 − θ ) + γij∗ θ 2
−
∗ ij
1 0
2
1 6
θ (1 − θ )dθ = − .
It leads to the following property. Property 2. In any subinterval [xi , xi+1 ], no matter what positive number the shape parameters αij∗ , βij∗ and γij∗ are, the integral weight coefficients ar ’s satisfy
1 1 < a0 < , 3 2
1 5 < a1 < , 2 6
−
1 < a2 < 0. 6
(35)
If αij∗ , βij∗ > γij∗ > 0, we have
ω2 (θ , αij∗ , βij∗ , γij∗ ) > −
γij∗ θ 2 (1 − θ ) θ 2 (1 − θ ) = − , ∗ ∗ γ (1 − θ )2 + γij θ (1 − θ ) + γij θ 2 1 − θ + θ2
−
∗ ij
1 0
1 θ 2 (1 − θ ) dθ = − 2 1 − θ + θ2
√ 3π . 9
Property 3. In any subinterval [xi , xi+1 ], If αij∗ , βij∗ > γij∗ > 0, the integral weight coefficients ar ’s satisfy
√ 1 3π < a0 < , 1− 9 2
√ 1 2 3π 1 < a1 < − + , 2 2 9
1 − 2
√ 3π < a2 < 0. 9
(36)
374
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377 Table 2 Data points. f(xi , yj )
x1 x2 x3 x4
=1 =2 =3 =4
y1 = 1
y2 = 2
y3 = 3
y4 = 4
4 2 3 1
2 3 2 4
3 2 4 2
1 4 2 3
Let D˜ = [0, 1; 0, 1]. Denote
ωrs (θ , αij∗ , βij∗ , γij∗ ; η, αij , βij , γij ) = ωr (θ , αij∗ , βij∗ , γij∗ )ωs (η, αij , βij , γij ),
ars =
D˜
ωrs (θ , αij∗ , βij∗ , γij∗ ; η, αij , βij , γij )dθ dη
ars are called the integral weight coefficients of BRIS Pij (x, y). From Eq. (26), one can yield that
2
r=0
2
s=0
ars = 1.
Property 4. If fi+r,j+s = 1, r, s = 0, 1, 2, then
D
Pij (x, y)dxdy = hi lj
2 2
ars fi+r,j+s = hi lj .
r=0 s=0
From Eq. (35) in Property 2, we obtain the following property. Property 5. In any subinterval Dij , no matter what positive number the parameters αij∗ , βij∗ , γij∗ , αij , βij , γij are, the integral weight coefficients ars ’s satisfy
1 1 1 5 1 1 < a00 < < a01 < − < a02 < 9 4 6 12 12 4 1 5 5 25 1 < a10 < < a11 < − < a12 < 0 6 12 4 36 36 5 1 1 1 < a20 < − < a21 < 0 0 < a22 < − 12 4 36 36 From Eq. (36), we have that Property 6. In any subinterval Dij , if αij∗ , βij∗ > γij∗ > 0, αij , βij > γij > 0, ars ’s satisfy
2 √ √ √ √ 1 1 1 3π 3π 3π 1 3π − < a01 < − + − < a02 < 0 < a00 < 9 4 2 18 4 9 4 18 2 √ √ √ √ 2 3π 2π 2 1 1 1 3π 3π 1 3π 1 − < a10 < − + < a11 < − + − < a12 < 0 − + 2 18 4 9 4 2 9 4 6 27 2 √ √ √ 2π 2 1 1 3π 3π 3π 1 − < a20 < 0 − + − < a21 < 0 0 < a12 < − 4 18 4 6 27 2 9
1−
6. Numerical examples Let D = [1, 4; 1, 4], and let the knots be equally spaced, the set of data points are given in Table 2, then BRIS Pij (x, y) can be constructed in [1, 3; 1, 3] for the given positive parameters αij∗ , βij∗ , γij∗ , αij , βij and γ ij . The corresponding graphes are shown in Fig. 5–8. In Fig. 5 each patch tilts heavily nearby two adjacent wedges. This phenomenon is similar to that in Fig. 1. In Fig. 6 and 7 the interpolating vertexes can be identified easily, and each patch is relatively plane because αij∗ βij∗ , γij∗ ; αij βij , γij in Fig. 6 and βij∗ αij∗ , γij∗ ; βij αij , γij in Fig. 7. The surface in Fig. 8 is a bicubic Hermite interpolating spline. In the proposed method we generate four shape control parameters in [15] to six parameters case. In [15] when γ ij α ij or γij∗ αij∗ , the interpolating surface tilts heavily. After adding two parameters βij∗ , βij and making βij∗ be close to or greater than γij∗ , β ij be close to or greater than γ ij , we can avoid the above inharmony of interpolating surface. A compared effect is shown in Fig. 9. In Fig. 9(a), according to [15], γij∗ = γij = 100 αij∗ = αij∗ = 1 while βij∗ = βij = 2, the corresponding surface tilts heavily. In Fig. 9(b), when we adjust the value of the parameters βij∗ and β ij to 200, the inharmonic phenomenon vanishes.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
6 3.0
z 4 2 2.5
0 2 1.0
2.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 5. BRIS with αij∗ = 1, βij∗ = 1, γij∗ = 10, αij = 1, βij = 1, γij = 100.
4 3.0
z
3 2.5 2 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 6. BRIS with αij∗ = 20, βij∗ = 1, γij∗ = 1, αij = 20, βij = 1, γij = 1.
4.0 3.5 z
3.0
3.0 2.5
2.5
2.0 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 7. BRIS with αij∗ = 1, βij∗ = 20, γij∗ = 1, αij = 1, βij = 20, γij = 1.
375
376
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377
4 z
3.0 3 2.5
2 2.0
1.0
y
1.5 1.5
2.0 x 2.5 3.0
1.0
Fig. 8. BRIS with αij∗ = 1, βij∗ = 2, γij∗ = 1, αij = 1, βij = 2, γij = 1.
4
6 3.0
z 4 2 2.5
0 2.0
1.0 1.5
y
z
3.0 3 2.5
2 2.0
1.0
y
1.5 1.5
2.0 x
1.5
2.0 x
2.5
2.5 3.0
1.0
∗ = (a) The case with αij = α∗ij = 1, βij = βij ∗ = 100 in [15] 2, γij = γij
3.0
1.0
∗ = (b) BRIS with αij = α∗ij = 1, βij = βij ∗ = 100 200, γij = γij
Fig. 9. Comparison the proposed six shape parameters case with the four parameters in [15] (a) The case with αij = αij∗ = 1, βij = βij∗ = 2, γij = γij∗ = 100 in [15] (b) BRIS with αij = αij∗ = 1, βij = βij∗ = 200, γij = γij∗ = 100.
Acknowledgments This work was supported by the NNSF of China (11472063) and the Provincial Natural Science Research Program of Higher Education Institutions of Anhui Province (KJ2013A194, KJ2013Z230). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
M. Abbas, A.A. Majid, J.M. Ali, Positivity-preserving C2 rational cubic spline interpolation, Sci. Asia 39 (2013) 208–213. M. Abbas, A.A. Majid, J.M. Ali, Positivity-preserving rational bi-cubic spline interpolation for 3D positive data, Appl. Math. Comput. 234 (2014) 460–476. F. Bao, Q. Sun, Q. Duan, Point control of the interpolating curve with a rational cubic spline, J. Vis. Commun. Image R. 20 (2009) 275–280. F. Bao, Q. Sun, J. Pan, Q. Duan, A blending interpolator with value control and minimal strain energy, Comput. Graph. 34 (2010) 119–124. P.E. Bézier, The Mathematical Basis of The UNISURF CAD System, Butterworth, London, 1986. C. de Boor, B-form basics, in: G. Farin (Ed.), Geometric Modeling, SIAM, Philadelphia, 1987, pp. 131–148. C.K. Chui, Multivariate Splines, SIAM, 1988. R. Delbourgo, Shape preserving interpolation to convex data by rational functions with quadratic numerator and linear denominator, IMA J. Numer. Anal. 9 (1989) 23–136. J. Deng, F. Chen, Y. Feng, Dimensions of spling spaces over t-meshes, J. Comput. Appl. Math. 194 (2006) 267–283. J. Deng, F. Chen, X. Li, Graph. Models 74 (2008) 76–86. P. Dierck, B. Tytgat, Generating the Bézier points of BETA-spline curve, Comput. Aided Geom. Des. 6 (1989) 279–291. Q. Duan, K. Djidjeli, W.G. Price, E.H. Twizell, The approximation properties of some rational cubic splines, Int. J. Comput. Math. 72 (1999) 155–166.
Y. Tao, D. Wang / Applied Mathematics and Computation 264 (2015) 366–377 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
377
Q. Duan, L. Wang, E.H. Twizell, A new bivariate rational interpolation based on function values, Inf. Sci. 166 (2004) 181–191. Q. Duan, L. Wang, E.H. Twizell, A new weighted rational cubic interpolation and its approximation, Appl. Math. Comput. 168 (2005) 990–1003. Q. Duan, H. Zhang, A. Liu, H. Li, A bivariate rational interpolation with a bi-quadratic denominator, J. Comput. Appl. Math. 195 (2006) 24–33. Q. Duan, Y. Zhang, E.H. Twizell, A bivariate rational interpolation and the properties, Appl. Math. Comput. 179 (2006) 190–199. Q. Duan, S. Li, F. Bao, E.H. Twizell, Hermite interpolation by piecewise rational surface, Appl. Math. Comput. 198 (2008) 59–72. Q. Duan, F. Bao, S. Du, E.H. Twizell, Local control of interpolating rational cubic spline curves, Comput. Aided Des. 41 (2009) 825–829. G. Farin, Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide, fifth ed., Morgan Kaufman, Menlo Park, 2002. J.A. Gregory, M. Sarfraz, P.K. Yuen, Interactive curve design using C2 rational splines, Comput. Graph. 18 (1994) 153–159. X. Han, Convexity preserving piecewise rational quartic interpolation, SIAM J. Numer. Anal. 46 (2008) 920–929. K. Konno, H. Chiyokura, An approach of designing and controlling free-form surfaces by using NURBS boundary Gregory patches, Comput. Aided Geom. Des. 13 (1996) 825–849. M.Z. Hussain, M. Sarfraz, Positivity-preserving interpolation of positive data by rational cubics, J. Comput. Appl. Math. 218 (2008) 446–458. R. Müler, Universal parametrization and interpolation on cubic surfaces, Comput. Aided Geom. Des. 19 (2002) 479–502. L. Piegl, On NURBS: a survey, IEEE Comput. Graph. Appl. 11 (1991) 55–71. M. Sarfraz, A C2 rational cubic spline which has linear denominator and shape control, Ann. Univ. Sci. Budapest 37 (1994) 53–62. T.W. Sederberg, J.M. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCS, ACM Trans. Graph. 22 (2003) 161–172. J. Tan, S. Tang, Composite schemes for multivariate blending rational interpolation, J. Comput. Appl. Math. 144 (2002) 263–275. R. Wang, Multivariate Spline Functions and Their Applications, Kluwer Academic Publishers, Beijing/New York/Dordrecht/Boston/London, 2001. Y. Zhang, Q. Duan, E.H. Twizell, Convexity control of a bivariate rational interpolating spline surfaces, Comput. Graph. 31 (2007) 679–687.