A weighted bivariate blending rational interpolation based on function values

A weighted bivariate blending rational interpolation based on function values

Applied Mathematics and Computation 217 (2011) 4644–4653 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homep...

633KB Sizes 5 Downloads 143 Views

Applied Mathematics and Computation 217 (2011) 4644–4653

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A weighted bivariate blending rational interpolation based on function values Wei-Xian Huang, Guo-Jin Wang ⇑ Department of Mathematics, Zhejiang University, Hangzhou 310027, China State Key Laboratory of CAD and CG, Zhejiang University, Hangzhou 310027, China

a r t i c l e

i n f o

Keywords: Computer aided geometric design Rational spline interpolation C1-continuous Bivariate interpolation Cubic spline

a b s t r a c t This paper presents a new weighted bivariate blending rational spline interpolation based on function values. This spline interpolation has the following advantages: firstly, it can modify the shape of the interpolating surface by changing the parameters under the condition that the values of the interpolating nodes are fixed; secondly, the interpolating function is C1-continuous for any positive parameters; thirdly, the interpolating function has a simple and explicit mathematical representation; fourthly, the interpolating function only depends on the values of the function being interpolated, so the computation is simple. In addition, this paper discusses some properties of the interpolating function, such as the bases of the interpolating function, the matrix representation, the bounded property, the error between the interpolating function and the function being interpolated. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction Spline interpolation is a useful tool for designing curves and surfaces in computer aided geometric design, hence in the past decade, there has been considerable interest in spline interpolation. For instance, literature [1] gives a method that interpolating the spatial points by a uniform cubic B-spline, and literatures [2,3] realize the smooth interpolation for surfaces. But this kind of spline interpolation has one disadvantage that given the interpolating data, the interpolating function is unique. In order to overcome this disadvantage, in recent years the rational spline, especially the rational cubic spline with linear and quadratic denominators, and its application to shape control, have received attention in literatures [4–10]. This type of interpolation has the advantage that it can modify the shape of curves and surfaces by changing the parameters for the given interpolating data. Literature [4] presents a bivariate rational interpolation based on function values with linear denominators. Literature [5] presents a univariate rational cubic interpolation based on function values with linear denominators. Literature [6] gives a bivariate rational interpolation based on function values and derivatives with linear denominators. Literature [7] discusses the bounded property and the shape control of a bivariate rational interpolating surface. Literature [8] combines the methods in [4,6], and gets a new spline interpolation with more parameters, so it is more effective to modify the shape of curves and surfaces. Also literatures [9,10] discusses the rational interpolation with quadratic denominators and properties of it. The bivariate rational interpolation in literatures [4,6–8], in order to maintain smoothness of the interpolating surface, the parameters of y-direction must be limited, hence it is inconvenient to modify the shape of surfaces. For this reason, this paper presents a new weighted bivariate rational bicubic spline interpolation based on function values, which is C1-continuous for any positive parameters. Besides, this kind of interpolation only depends on the values of the function to be ⇑ Corresponding author at: Institute of Computer Images and Graphics, Zhejiang University, Hangzhou 310027, China. E-mail address: [email protected] (G.-J. Wang). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.11.016

4645

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

interpolated, without partial derivatives, which is convenient in the engineering. In addition, the paper discusses some properties of this type of interpolation, such as the bases, the matrix representation, the bounded property, the error, and gives some examples. The rest of this paper is organized as follows. In the next section we present a new weighted bivariate rational interpolation based on function values. In Section 3, we present the bases of the interpolation and the matrix representation of the interpolating function. In Section 4, we discuss the bounded property of this type of interpolation. In Section 5, we present two methods for the error estimation of the interpolation. Finally, in Section 6, we present some examples of this kind of interpolation and conclude in Section 7. 2. Interpolation Let X : ½a; b; c; d be the plane region, and fðxi ; yi ; fi;j Þ; i ¼ 1; 2; . . . ; n; n þ 1; j ¼ 1; 2; . . . ; m; m þ 1g be a given set of data points, where a ¼ x1 < x2 <    < xn < xnþ1 ¼ b; c ¼ y1 < y2 <    < ym < ymþ1 ¼ d are the knot spacings. Let yy i hi ¼ xiþ1  xi ; lj ¼ yjþ1  yj , and for any point ðx; yÞ 2 ½xi ; xiþ1 ; yj ; yjþ1  in the xy-plane, let h ¼ xx ; g ¼ lj j . First, for each hi y ¼ yj ; j ¼ 1; 2; . . . ; m þ 1, we construct the x-direct interpolating curve [4]:

pi;j ðxÞ ¼ p ðxÞj½xi ;xiþ1  ¼

ð1  hÞ3 ai;j fi;j þ hð1  hÞ2 V i;j þ h2 ð1  hÞW i;j þ h3 fiþ1;j ; ð1  hÞai;j þ h

i ¼ 1; 2; . . . ; n  1;

where V i;j ¼ ðai;j þ 1Þfi;j þ ai;j fiþ1;j ; W i;j ¼ ðai;j þ 2Þfiþ1;j  hi Diþ1;j , with ai;j > 0 and Di;j ¼ rational cubic interpolation based on function values which satisfies:

pi;j ðxi Þ ¼ fi;j ;

pi;j ðxiþ1 Þ ¼ fiþ1;j ;

@pi;j ðxi Þ @x

¼ Di;j ;

@pi;j ðxiþ1 Þ @x

fiþ1;j fi;j . hi

ð1Þ

This interpolation is called the

¼ Diþ1;j :

Clearly, the interpolating function pi;j ðxÞ on ½xi ; xiþ1  is unique for the given data ðxr ; f ðxr ; yj ÞÞ; r ¼ i; i þ 1; i þ 2 and the parameter ai;j . Based on the x-direct interpolation pi;j ðxÞ, we now construct the bivariate C1-continuous rational interpolation as follows:

p1i;j ðx; yÞ ¼ p1 ðx; yÞj½xi ;xiþ1 ;yj ;yjþ1  ¼ ð1  gÞ3 pi;j ðxÞ þ gð1  gÞ2 V i;j þ g2 ð1  gÞW i;j þ g3 pi;jþ1 ðxÞ; i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1;

ð2Þ

where V i;j ¼ 2pi;j ðxÞ þ pi;jþ1 ðxÞ; W i;j ¼ 3pi;jþ1 ðxÞ  lj Di;jþ1 ðxÞ, with Di;j ðxÞ ¼ ðpi;jþ1 ðxÞ  pi;j ðxÞÞ=lj . p1i;j ðx; yÞ is called the bivariate rational bicubic interpolation based on function values which satisfies:

p1i;j ðxr ; ys Þ ¼ f ðxr ; ys Þ ¼ fr;s ;

@p1i;j ðx; ys Þ @y

¼ Di;s ðxÞ;

r ¼ i; i þ 1;

s ¼ j; j þ 1:

Obviously, the interpolating function p1i;j ðx; yÞ on ½xi ; xiþ1 ; yj ; yjþ1  is unique ðxr ; ys ; f ðxr ; ys ÞÞ; r ¼ i; i þ 1; i þ 2; s ¼ j; j þ 1; j þ 2 and the parameters ai;s ; s ¼ j; j þ 1; j þ 2.

for

the

given

data

Lemma 1. For any positive parameter ai;j ; ði ¼ 1; 2 . . . ; n  1; j ¼ 1; 2; . . . ; m  1Þ, the interpolating function p1 ðx; yÞ on the region ½x1 ; xn ; y1 ; yn  is C1-continuous. Proof. By Theorem 1 in literature [4], obviously, Lemma 1 holds. h In the above we firstly construct x-direct interpolation. Similarly, for each x ¼ xi ; i ¼ 1; 2; . . . ; n þ 1, we can construct ydirect interpolation as follows:

qi;j ðyÞ ¼ q ðyÞj½yj ;yjþ1  ¼

ð1  gÞ3 bi;j fi;j þ gð1  gÞ2 T i;j þ g2 ð1  gÞZ i;j þ g3 fi;jþ1 ; ð1  gÞbi;j þ g

j ¼ 1; 2; . . . ; m  1:

ð3Þ

where T i;j ¼ ðbi;j þ 1Þfi;j þ bi;j fi;jþ1 ; Z i;j ¼ ðbi;j þ 2Þfi;jþ1  lj di;jþ1 , with di;j ¼ ðfi;jþ1  fi;j Þ=lj , and bi;j > 0. Then we construct the bivariate rational interpolation as follows:

p2i;j ðx; yÞ ¼ p2 ðx; yÞj½xi ;xiþ1 ;yj ;yjþ1  ¼ ð1  hÞ3 qi;j ðyÞ þ hð1  hÞ2 T i;j þ h2 ð1  hÞZ i;j þ h3 qiþ1;j ðyÞ; i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1;

ð4Þ

where T i;j ¼ 2qi;j ðyÞ þ qiþ1;j ðyÞ; Z i;j ¼ 3qiþ1;j ðyÞ  hidiþ1;j ðyÞ, with di;j ðyÞ ¼ ðqiþ1;j ðyÞ  qi;j ðyÞÞ=hi . p2i;j ðx; yÞ is also called the bivariate rational bicubic interpolation based on function values which satisfies:

p2i;j ðxr ; ys Þ ¼ f ðxr ; ys Þ ¼ fr;s ;

@p2i;j ðxr ; yÞ @x

¼ dr;j ðyÞ;

r ¼ i; i þ 1;

s ¼ j; j þ 1:

4646

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

Clearly, the interpolating function p2i;j ðx; yÞ on ½xi ; xiþ1 ; yj ; yjþ1  is unique ys ÞÞ; r ¼ i; i þ 1; i þ 2; s ¼ j; j þ 1; j þ 2 and the parameters br;j ; r ¼ i; i þ 1; i þ 2. For k 2 ½0; 1, let

pi;j ðx; yÞ ¼ pðx; yÞj½xi ;xiþ1 ;yj ;yjþ1  ¼ kp1i;j ðx; yÞ þ ð1  kÞp2i;j ðx; yÞ;

i ¼ 1; 2; . . . ; n  1;

for

the

given

j ¼ 1; 2; . . . ; m  1:

data

ðxr ; ys ; f ðxr ;

ð5Þ

Then p(x, y) is called the weighted bivariate rational bicubic interpolation which satisfies:

@pi;j ðxr ; ys Þ ¼ Dr;s ¼ di;j ðys Þ; @x

pi;j ðxr ; ys Þ ¼ fr;s ;

@pi;j ðxr ; ys Þ ¼ dr;s ¼ Di;j ðxr Þ; @y

ð6Þ

for r ¼ i; i þ 1; s ¼ j; j þ 1. Hence by Lemma 1 we can obtain the following theorem: Theorem 1. For given ai;j > 0; bi;j > 0; i ¼ 1; 2; . . . ; n; n  1; j ¼ 1; 2; . . . ; m; m  1, and ðxi ; yi ; fi;j Þ; i ¼ 1; 2; . . . ; n; n þ 1; j ¼ 1; 2; . . . ; m; m þ 1, and k 2 ½0; 1, the interpolating function p(x, y) which is defined by Eq. (5) and satisfies Eq. (6) is C1-continuous on the region ½x1 ; xn ; y1 ; yn . 3. The bases and the matrix representation Let 1;1

1;1

bi;j ¼ bi;j ðai;j ; hÞ ¼

1;2

ð1  hÞ2 ðai;j þ hÞ ; ð1  hÞai;j þ h

1;2

biþ1;j ¼ biþ1;j ðai;j ; hÞ ¼

1;3

1;3

biþ2;j ¼ biþ2;j ðai;j ; hÞ ¼

hð1  hÞai;j þ 2h2  h3 þ h2 ð1  hÞ hhiþ1i ð1  hÞai;j þ h h2 ð1  hÞ hhiþ1i ð1  hÞai;j þ h

;

;

2 1;1 c1;1 i;j ¼ ci;j ðgÞ ¼ ð1  gÞ ð1 þ gÞ;

1;2 2 2 c1;2 i;jþ1 ¼ c i;jþ1 ðgÞ ¼ gð1 þ g  g Þ þ g ð1  gÞ

1;3 2 c1;3 i;jþ2 ¼ c i;jþ2 ðgÞ ¼ g ð1  gÞ

2;1

2;1

bi;j ¼ bi;j ðbi;j ; gÞ ¼

2;2

2;2

2;3

2;3

lj ; ljþ1

lj ; ljþ1

ð1  gÞ2 ðbi;j þ gÞ ; ð1  gÞbi;j þ g

bi;jþ1 ¼ bi;jþ1 ðbi;j ; gÞ ¼

bi;jþ2 ¼ bi;jþ2 ðbi;j ; gÞ ¼

lj gð1  gÞbi;j þ 2g2  g3 þ g2 ð1  gÞ ljþ1 ; ð1  gÞbi;j þ g

g2 ð1  gÞ l

lj jþ1

ð1  gÞbi;j þ g

;

2 2;1 c2;1 i;j ¼ ci;j ðhÞ ¼ ð1  hÞ ð1 þ hÞ;

2;2 2 2 c2;2 iþ1;j ¼ c iþ1;j ðhÞ ¼ hð1 þ h  h Þ þ h ð1  hÞ

2;3 2 c2;3 iþ2;j ¼ c iþ2;j ðhÞ ¼ h ð1  hÞ

hi ; hiþ1

hi ; hiþ1

then the interpolating function pi;j ðx; yÞ defined by Eq. (5) can be written as follows:

pi;j ðx; yÞ ¼

jþ2 iþ2 X X r¼i

s¼j

ar;s fr;s ;

i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1;

ð7Þ

4647

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

where 1;1

2;1

2;1 ai;j ¼ bi;j ðai;j ; hÞc1;1 i;j ðgÞk þ bi;j ðbi;j ; gÞc i;j ðhÞð1  kÞ; 1;2

2;1

1;3

2;1

2;2 aiþ1;j ¼ biþ1;j ðai;j ; hÞc1;1 i;j ðgÞk þ biþ1;j ðbiþ1;j ; gÞc iþ1;j ðhÞð1  kÞ; 2;3 aiþ2;j ¼ biþ2;j ðai;j ; hÞc1;1 i;j ðgÞk þ biþ2;j ðbiþ2;j ; gÞc iþ2;j ðhÞð1  kÞ; 1;1

2;2

1;2 2;1 ðgÞk þ bi;jþ1 ðbi;j ; gÞci;j ðhÞð1  kÞ; ai;jþ1 ¼ bi;jþ1 ðai;jþ1 ; hÞci;jþ1 1;2

2;2

1;3

2;2

2;2 aiþ1;jþ1 ¼ biþ1;jþ1 ðai;jþ1 ; hÞc1;2 i;jþ1 ðgÞk þ biþ1;jþ1 ðbiþ1;j ; gÞc iþ1;j ðhÞð1  kÞ; 2;3 aiþ2;jþ1 ¼ biþ2;jþ1 ðai;jþ1 ; hÞc1;2 i;jþ1 ðgÞk þ biþ2;jþ1 ðbiþ2;j ; gÞc iþ2;j ðhÞð1  kÞ; 1;1

2;3

1;3 2;1 ðgÞk þ bi;jþ2 ðbi;j ; gÞci;j ðhÞð1  kÞ; ai;jþ2 ¼ bi;jþ2 ðai;jþ2 ; hÞci;jþ2 1;2

2;3

1;3

2;3

2;2 aiþ1;jþ2 ¼ biþ1;jþ2 ðai;jþ2 ; hÞc1;3 i;jþ2 ðgÞk þ biþ1;jþ2 ðbiþ1;j ; gÞc iþ1;j ðhÞð1  kÞ; 2;3 aiþ2;jþ2 ¼ biþ2;jþ2 ðai;jþ2 ; hÞc1;3 i;jþ2 ðgÞk þ biþ2;jþ2 ðbiþ2;j ; gÞc iþ2;j ðhÞð1  kÞ:

far;s ; r ¼ i; i þ 1; i þ 2; s ¼ j; j þ 1; j þ 2g is called the bases of the interpolating function which satisfies: jþ2 iþ2 X X r¼i

ar;s ¼ 1:

ð8Þ

s¼j

Here suppose that ai;j ¼ ai;jþ1 ¼ ai;jþ2 ; bi;j ¼ biþ1;j ¼ biþ2;j , then it is easy to obtain the matrix representation for the interpolating function. Let

0

fi;j B F i;j ¼ @ fiþ1;j fiþ2;j B1i;j ¼ B2i;j ¼

 

fi;jþ1 fiþ1;jþ1 fiþ2;jþ1

1;1

biþ1;j

2;1

bi;jþ1

bi;j bi;j

fi;jþ2

1;2

2;2

1

C fiþ1;jþ2 A; fiþ2;jþ2

 1;3 biþ2;j ; 2;3

bi;jþ2

T

C 1i;j ¼ ;



C 2i;j ¼

c1;1 i;j



c1;2 i;jþ1

2;1 ci;j

c2;2 iþ1;j

c1;3 i;jþ2

T

;

 2;3 ciþ2;j ;

then the interpolation pi;j ðx; yÞ defined by Eq. (5) can be represented by matrices as follows:

pi;j ðx; yÞ ¼ kB1i;j F i;j C 1i;j þ ð1  kÞC 2i;j F i;j B2i;j ;

i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1:

ð9Þ

4. The bounded property of the interpolation r¼iþ2;s¼jþ2 Let M ¼ maxr¼i;s¼j jfr;s j, we can obtain the following theorem:

Theorem 2. For any positive parameters ai;s ; br;j ; s ¼ j; j þ 1; j þ 2; r ¼ i; i þ 1; i þ 2, the upper boundary of the interpolating function defined by Eq. (5) on ½xi ; xiþ1 ; yj ; yjþ1  is determined by the following inequality:

 2  2  2  2 !1=2 8lj hi 8hi lj jpi;j ðx; yÞj 6 M 1þ 1þ þ 1þ 1þ ; 27ljþ1 2hiþ1 27hiþ1 2ljþ1 ¼ 1; 2; . . . ; m  1:

i ¼ 1; 2; . . . ; n  1;

j ð10Þ

Especially when hi ¼ hiþ1 ; lj ¼ ljþ1 , Then

pffiffiffi 35 2M : jpi;j ðx; yÞj 6 18

ð11Þ

4648

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

Proof. By the bases representation (Eq. (7)) of the interpolation pi;j ðx; yÞ and k 2 ½0; 1; h 2 ½0; 1; g 2 ½0; 1, it follows that:

  jþ2 jþ2 jþ2  X X iþ2 X iþ2 X iþ2 X X   jpi;j ðx; yÞj ¼  ar;s fr;s  6 jar;s fr;s j 6 M jar;s j;   r¼i s¼j r¼i s¼j r¼i s¼j 1;riþ1

By the definitions of c1;sjþ1 and br;s i;s

; r ¼ i; i þ 1; i þ 2; s ¼ j; j þ 1; j þ 2, it is easy to obtain that:

jþ2   X lj 8lj  1;sjþ1  61þ ; ci;s  ¼ 1 þ 2g2 ð1  gÞ l 27l jþ1 jþ1 s¼j

 iþ2  X  1;riþ1  br;s ¼1þ r¼i

2h2 ð1  hÞhi hi 61þ : ½ð1  hÞai;s þ hhiþ1 2hiþ1

Similarly, iþ2   X hi 8hi  2;riþ1  61þ ; cr;j  ¼ 1 þ 2h2 ð1  hÞ h 27h iþ1 iþ1 r¼i

jþ2   X  2;sjþ1  br;s ¼1þ s¼j

2g2 ð1  gÞlj lj 61þ ; ½ð1  gÞbr;j þ gljþ1 2ljþ1

for r ¼ i; i þ 1; i þ 2; s ¼ j; j þ 1; j þ 2. Then jþ2 iþ2 X X r¼i

s¼j

      8lj hi 8hi lj 1þ þ ð1  kÞ 1 þ 1þ jar;s j 6 k 1 þ 27ljþ1 2hiþ1 27hiþ1 2ljþ1 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi      2  2ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 8lj hi 8hi lj 2 2 1þ þ 1þ 1þ 1þ 6 k þ ð1  kÞ 27ljþ1 2hiþ1 27hiþ1 2ljþ1 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  2  2  2ffi 8lj hi 8hi lj : 1þ þ 1þ 1þ 6 1þ 27ljþ1 2hiþ1 27hiþ1 2ljþ1

Hence Eq. (10) holds, and Eq. (11) can be derived directly from Eq. (10).

h

5. The error of the interpolation Here we give two methods for the error estimation of the interpolation. The first is based on the Taylor series, and the second is based on the Peano–Kernel. 5.1. The error estimation based on Taylor series Suppose ðx; yÞ 2 ½xi ; xiþ1 ; yj ; yjþ1 , and f ðx; yÞ 2 C 1 is the function to be interpolated, by Taylor series we have:

  @ @ f ðn1 ; n2 Þ; f ðx; yÞ ¼ f ðxr ; ys Þ þ ðx  xr Þ þ ðy  ys Þ @x @y

r ¼ i; i þ 1; i þ 2;

s ¼ j; j þ 1; j þ 2;

where n1 2 ½xi ; xiþ2 ; n2 2 ½yj ; yjþ2 . Let

    @f      ¼ max @f ðx; yÞ; @x x2½xi ;xiþ2   @x  Then

     @f      ¼ max @f ðx; yÞ: @y y2½yj ;yjþ2   @y 

    @f   @f  jf ðx; yÞ  f ðxr ; ys Þj 6 ðhi þ hiþ1 Þ  þ ðlj þ ljþ1 Þ ; @x @y

r ¼ i; i þ 1; i þ 2;

s ¼ j; j þ 1; j þ 2:

Suppose pi;j ðx; yÞ is defined by Eq. (5), and it can be rewritten by Eq. (7). By Eq. (8) it follows that:

f ðx; yÞ  pi;j ðx; yÞ ¼

jþ2 iþ2 X X r¼i

ar;s ðf ðx; yÞ  fr;s Þ;

s¼j

hence

    iþ2 jþ2  @f   @f  X X jf ðx; yÞ  pi;j ðx; yÞj 6 ðhi þ hiþ1 Þ  þ ðlj þ ljþ1 Þ  jar;s j: @x @y r¼i s¼j According to Theorem 2 and the proof of it, the following theorem can be obtained:

ð12Þ

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

4649

Theorem 3. Suppose pi;j ðx; yÞ is defined by Eq. (5) , and f ðx; yÞ 2 C 1 is the function to be interpolated, for any positive parameters ai;s ; br;j ; s ¼ j; j þ 1; j þ 2; r ¼ i; i þ 1; i þ 2, the error of the interpolation on ½xi ; xiþ1 ; yj ; yjþ1  is determined by the following inequality:

     @f   @f  jf ðx; yÞ  pi;j ðx; yÞj 6 ðhi þ hiþ1 Þ  þ ðlj þ ljþ1 Þ  @x @y " 2  2  2  2 #1=2 8lj hi 8hi lj  1þ 1þ þ 1þ 1þ ; 27ljþ1 2hiþ1 27hiþ1 2ljþ1 i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1:

ð13Þ

Especially when hi ¼ hiþ1 ; lj ¼ ljþ1 ,

jf ðx; yÞ  pi;j ðx; yÞj 6

pffiffiffi      @f   @f  35 2 hi   þ lj   : @x @y 9

ð14Þ

5.2. The error estimation based on the Peano–Kernel Suppose ðx; yÞ 2 ½xi ; xiþ1 ; yj ; yjþ1 , and f ðx; yÞ 2 C 2 is the function to be interpolated, pi;j ðx; yÞ is defined by Eq. (5). Here let  

 

@p

@f ðx;yÞ

¼ maxx2½x ;x  @pi;j ðx;yÞ;

@f ¼ maxx2½x ;x   . @x

i

iþ1

@x

@x

i

@x

iþ1

According to Theorem 4.1.1 in literature [11], we have

jf ðx; yÞ  pi;j ðx; yÞj 6 jf ðx; yÞ  f ðxi ; yÞj þ jf ðxi ; yÞ  pi;j ðxi ; yÞj þ jpi;j ðxi ; yÞ  pi;j ðx; yÞj          Z yjþ2 2  @f  @p @f  @p @ f ðxi ; sÞ   R ½ðy  s Þ d s 6 hi   þ   þ jf ðxi ; yÞ  pi;j ðxi ; yÞj ¼ hi   þ   þ   y þ 2   yj @x @x @x @x @y  Z      2   y jþ2 @f  @p @ f ðxi ; yÞ 6 hi   þ   þ  jRy ½ðy  sÞþ jds;   @y2  yj @x @x

2

2 



  where Ry ½ðy  sÞþ  ¼ kR1y ½ðy  sÞþ  þ ð1  kÞR2y ½ðy  sÞþ ; @ f@yðx2i ;yÞ ¼ maxy2½yj ;yjþ2  @ f@yðx2i ;yÞ, with

8   8 1 > > ðy  sÞ  c1;2 ðyjþ1  sÞ þ c1;3 ðyjþ2  sÞ ; yj < s < y; > i;jþ1 i;jþ2 > > < r 1 ðsÞ; yj < s < y; <   1 1;2 1;3 Ry ½ðy  sÞþ  ¼  c ðyjþ1  sÞ þ c ðyjþ2  sÞ ; ¼ r 12 ðsÞ; y < s < yjþ1 ; y < s < yjþ1 ; i;jþ1 i;jþ2 > > : 1 > > > r 3 ðsÞ; yjþ1 < s < yjþ2 ; : c1;3 ðy  sÞ; y < s < y ; jþ1 jþ2 i;jþ2 jþ2 8   2;2 2;3 8 2 > > > ðy  sÞ  bi;jþ1 ðyjþ1  sÞ þ bi;jþ2 ðyjþ2  sÞ ; yj < s < y; > > <  < r 1 ðsÞ; yj < s < y;  2;3 r 2 ðsÞ; y < s < yjþ1 ; R2y ½ðy  sÞþ  ¼  b2;2 ¼ y < s < yjþ1 ; i;jþ1 ðyjþ1  sÞ þ bi;jþ2 ðyjþ2  sÞ ; > > 2 > : > > r 23 ðsÞ; yjþ1 < s < yjþ2 : : b2;3 ðy  sÞ; yjþ1 < s < yjþ2 ; i;jþ2 jþ2 2g½ð1gÞbi;j þg 2 i;j þ2gg

Let wðg; bi;j Þ ¼ ð1gÞb

Z

y

yi

þ

gðg2 2g1Þ½ð1gÞbi;j þ2gg2  þ 4g. ð1gÞbi;j þg

By simple integral calculation, we can obtain that

h

i h i  1  o kr ðsÞ þ ð1  kÞr 2 ðsÞds ¼ g þ kc1;2 þ ð1  kÞb2;2 ð2  gÞ þ kc1;3 þ ð1  kÞb2;3 ð2  gÞ þ 2ljþ1 =lj gl2 =2 i;jþ1 i;jþ2 j 1 1 i;jþ1 i;jþ2 n

2

¼ B1 ðg; k; bi;j ; ljþ1 =lj Þlj : ð15Þ Z y

yjþ1

 1  kr ðsÞ þ ð1  kÞr 2 ðsÞds 6 2 2

Z y

yjþ1

jkr 12 ðsÞjds þ

Z

yjþ1 y

jð1  kÞr22 ðsÞjds

 2 gð1  6g2 þ 6g3 þ 2g4  4g5 þ g6 Þ 1;3 2ð1  2g þ g3 Þljþ1 lj ¼k þ c þ ð1 i;jþ2 ð1 þ g  g2 Þ ð1 þ g  g2 Þlj 2 " # 2 2 2 2;3 2ð1  gÞ bi;j  2g ð1  gÞ ljþ1 lj  kÞ wðg; bi;j Þ þ bi;jþ2 lj 2 ð1  gÞbi;j þ 2g  g2 2

¼ B2 ðg; k; bi;j ; ljþ1 =lj Þlj ;

ð16Þ

4650

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

Z

yjþ2

yjþ1

  2;3 2 2 jkr 13 ðsÞ þ ð1  kÞr 23 ðsÞjds ¼ kc1;3 i;jþ2  ð1  kÞbi;jþ2 ljþ1 =2 ¼ B3 ðg; k; bi;j ; ljþ1 =lj Þlj ;

ð17Þ

Then the following theorem can be derived: Theorem 4. Suppose pi;j ðx; yÞ is defined by Eq. (5) , and f ðx; yÞ 2 C 2 is the function to be interpolated, for any positive parameters ai;s ; br;j ; s ¼ j; j þ 1; j þ 2; r ¼ i; i þ 1; i þ 2, the error of the interpolation on the subregion ½xi ; xiþ1 ; yj ; yjþ1  is determined by the following inequality:

      2  @f  @p @ f ðxi ; yÞ x 2 jf ðx; yÞ  pi;j ðx; yÞj 6 hi   þ   þ  Bi;j lj ; 2  @y  @x @x

i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1;

ð18Þ

where Bxi;j ¼ max0
2

2 

@ f ðx ;yÞ

@ f ðx ;yÞ Similarly, let @yiþ1

¼ maxy2½yj ;yjþ2   @yiþ1 , we have 2 2 Theorem 5. Suppose pi;j ðx; yÞ is defined by Eq. (5) , and f ðx; yÞ 2 C 2 is the function to be interpolated, for any positive parameters ai;s ; br;j ; s ¼ j; j þ 1; j þ 2; r ¼ i; i þ 1; i þ 2, the error of the interpolation on the subregion ½xi ; xiþ1 ; yj ; yjþ1  is determined by the following inequality:

      2  @f  @p @ f ðxiþ1 ; yÞ x 2 jf ðx; yÞ  pi;j ðx; yÞj 6 hi   þ   þ  Biþ1;j lj ; 2   @x @x @y

i ¼ 1; 2; . . . ; n  1;

j ¼ 1; 2; . . . ; m  1;

ð19Þ

where Bxiþ1;j ¼ max0


 

 

2

2 



@p ðx;yÞ @f

@f ðx;yÞ

@ f ðx;y Þ

@ f ðx;y Þ ¼ maxy2½yj ;yjþ1   i;j@y ; @y Let @p

¼ maxy2½yj ;yjþ1   @y  and @x2 s ¼ maxx2½xi ;xiþ2   @x2 s ; s ¼ j; j þ 1. Then by the symme@y

try, we can obtain the following theorem: Theorem 6. Suppose pi;j ðx; yÞ is defined by Eq. (5) , and f ðx; yÞ 2 C 2 is the function to be interpolated, for any positive parameters

ai;s ; br;j ; s ¼ j; j þ 1; j þ 2; r ¼ i; i þ 1; i þ 2, the error of the interpolation on the subregion ½xi ; xiþ1 ; yj ; yjþ1  is determined by the following inequality:

      2   @f  @p @ f ðx; ys Þ y 2         jf ðx; yÞ  pi;j ðx; yÞj 6 lj   þ   þ  B h ;  @x2  i;s i @y @y

s ¼ j; j þ 1;

i ¼ 1; 2; . . . ; n  1;

¼ 1; 2; . . . ; m  1;

j ð20Þ

Byi;s

where ¼ max0
B3 ðg; k; bi;j ; 1Þ 6 0:14124;

ð21Þ

B2 ðg; k; bi;j ; 1Þ 6 0:27386;

ð22Þ

B1 ðg; k; bi;j ; 1Þ 6 0:33799;

ð23Þ

Proof

" B3 ðg; k; bi;j ; 1Þ ¼ 1=2 kðg2  g3 Þ þ ð1  kÞ

#



g2  g3 g2  g3 6 1=2 kðg2  g3 Þ þ ð1  kÞ ð1  gÞbi;j þ g g



pffiffiffiffiffiffiffiffiffiffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1=2ðg  g2 Þðgk þ ð1  kÞÞ 6 1=2ðg  g2 Þ 1 þ g2 k2 þ ð1  kÞ2 6 1=2ðg  g2 Þ 1 þ g2 6 0:14124; wðg; bi;j Þ 6

4g  7g2 þ 4g3  g4 ; 2g

therefore it is easy to know that:

4651

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653



gð1  6g2 þ 6g3 þ 2g4  4g5 þ g6 Þ 2ð1  2g þ g3 Þðg2 þ g3 Þ þ 1 þ g  g2 1 þ g  g2 ! 2 2 3 4 2 4g  7g þ 4g  g 2g ð1  gÞ þ 1=2ð1  kÞ þ 2g 2g  3 5 6 1  2g þ 2g  2g þ g 4  5g þ g3 ¼ 1=2 kg þ ð1  kÞ g 1 þ g  g2 2g

B2 ðg; k; bi;j ; 1Þ 6 1=2k



sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  2g þ 2g3  2g5 þ g6 Þ2 ð4  5g þ g3 Þ2 2 2 6 1=2 k þ ð1  kÞ g þ ð1 þ g  g2 Þ2 ð2  gÞ2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  2g þ 2g3  2g5 þ g6 Þ2 ð4  5g þ g3 Þ2 6 0:27386; þ 6 ð1=2Þg ð1 þ g  g2 Þ2 ð2  gÞ2

B1 ðg; k; bi;j ; 1Þ 6 1=2½g2 þ kgððg þ 2g2  2g3 Þð2  gÞ þ ðg2 þ g3 Þð4  gÞÞ þ ð1  kÞgð3g  2g2 Þð2  gÞ ¼ 1=2½g2 þ kgð2g  g2  g3 þ g4 Þ þ ð1  kÞgð3g  2g2 Þð2  gÞ  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 1=2 g2 þ k2 þ ð1  kÞ2 g2 ð2  g  g2 þ g3 Þ2 þ ð2g2  7g þ 6Þ2  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  6 1=2 g2 þ g2 ð2  g  g2 þ g3 Þ2 þ ð2g2  7g þ 6Þ2 6 0:33799:

Remark 1. All the maximum values are calculated in Maple 12. Remark 2. According to the result of Theorem 7, the error estimation of the interpolation defined by Eq. (18) is stable. Similarly, the error estimations defined by Eqs. (19) and (20) are also stable.

6. Numerical examples qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Suppose the interpolated function is f ðx; yÞ ¼ 1  ð1  xÞ2  ð1  yÞ2 ; ðx; yÞ 2 ½0:5; 1:5; 0:5; 1:5, let hi ¼ lj ¼ 1=5; i ¼ 1; 2; 3; 4; 5; j ¼ 1; 2; 3; 4; 5, then xi ¼ 0:5 þ ði  1Þ=5; yj ¼ 0:5 þ ðj  1Þ=5; i ¼ 1; 2; 3; 4; 5; 6; j ¼ 1; 2; 3; 4; 5; 6. Also let ai;j ¼ 0:5þ ði  1Þ=2 þ ðj  1Þ=2; bi;j ¼ 0:5 þ ði  1Þ=2 þ ðj  1Þ=2; k ¼ 0:5. Then the surface f(x, y) is shown by the left picture of Fig. 1, the surface of interpolating function is shown by the central picture of Fig. 1, and the surface of the error is shown by the right picture of Fig. 1. Changing the parameters, let ai;j ¼ 1555 þ 2i þ 2j; bi;j ¼ 3555 þ 3i þ 3j, then the surface of the interpolation, written as p(x, y), and the surface of the error are shown by Fig. 2. From Fig. 1, it can be seen that the error of the interpolation is smaller than ±8  103. This means that under the condition of C1-continuous, a good approximation can be achieved. From Fig. 2, we can know that the error of the interpolation is smaller than ±15  103. This means that with large parameters, the error changes little, hence the interpolation is stable with parameters, therefore it can be used to modify the shape of the interpolating surface appropriately. In the following, we present a simple example to demonstrate how to construct the interpolating surface by solving the parameters ai;j ; bi;j ; k. For simplicity, we assume that ai;j ¼ bi;j  a; k ¼ 0:5. Suppose the interpolated function is f ðx; yÞ ¼

−3

x 10 0.95

0.95

2

0.9

0.9

0

0.85

0.85

−2

0.8

0.8

0.75

0.75

0.7 1.5 1

y

0.5

0.4

0.6

0.8

x

1

1.2

1.4

z

4

z

1

z

1

−4 −6 −8 1.5

0.7 1.5 1

y

0.5

0.4

0.6

0.8

x

1

1.2

1.4 1

y

0.5

0.4

0.6

0.8

1

1.2

x

Fig. 1. The surfaces when ai;j ¼ bi;j ¼ 0:5 þ ði  1Þ=2 þ ðj  1Þ=2: (left) function’s surface, (center) interpolating surface, (right) error’s surface.

1.4

4652

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

−3

x 10 5

1 0.95

0

z

z

0.9 0.85 0.8

−5 −10

0.75 0.7 1.5

y

1 0.5

0.8

0.6

0.4

1

1.2

1.4

−15 1.5

y

x

1 0.5

0.4

0.6

0.8

1

1.2

1.4

x

Fig. 2. The surfaces when ai;j ¼ 1555 þ 2i þ 2j; bi;j ¼ 3555 þ 3i þ 3j: (left) interpolating surface, (right) error’s surface.

0.08

z

0.06 0.04 0.02 0 0.7

0.65

0.6

x

0.55

0.5 0.7

0.65

0.6

0.55

0.5

y

Fig. 3. Central point interpolating surface when s = 0.0245.

ð1  xÞ2 ð1  yÞ2 ; ðx; yÞ 2 ½0:5; 0:9; 0:5; 0:9, and let hi ¼ lj ¼ 1=5; i ¼ 1; 2; j ¼ 1; 2, then xi ¼ 0:5 þ ði  1Þ=5; yj ¼ 0:5þ ðj  1Þ=5; i ¼ 1; 2; 3; j ¼ 1; 2; 3. Suppose we need to determine the parameter a such that the interpolating surface interpolates the central point (0.6, 0.6, s), s is some real number. The matrices in Eq. (9) at the domain point (0.6, 0.6) can be expressed as follows:

0

F 1;1

1 f ð0:5; 0:5Þ f ð0:5; 0:7Þ f ð0:5; 0:9Þ B C ¼ @ f ð0:7; 0:5Þ f ð0:7; 0:7Þ f ð0:7; 0:9Þ A; f ð0:9; 0:5Þ f ð0:9; 0:7Þ f ð0:9; 0:9Þ

B11;1 ¼ ðB21;1 ÞT ¼ ðC 11;1 ÞT ¼ C 21;1 ¼

0:5aþ0:25 aþ1

3 8

3 4

1 8

0:5aþ1 aþ1



0:25 aþ1



;

:

Hence, we can obtain the following equation:

p1;1 ð0:6; 0:6Þ ¼ B11;1 F 1;1 C 11;1 þ C 21;1 F 1;1 B21;1 ¼ 2s: If we take s = f(0.6, 0.6), then the parameter can be solved as a = 1; If we take s = 0.0245, then the parameter can be solved as a = 0.1 (see Fig. 3). 7. Conclusion This paper presents a new weighted bivariate rational bicubic interpolation based on function values which is defined by Theorem 1, and discusses some properties of it under some conditions, see Theorems 2–7. In the numerical examples, we can see that the error of the interpolation is much small and the interpolation is stable with parameters. Besides, we can modify the shape of the interpolating surface appropriately, hence it is useful in computer aided geometric design. The other properties of the interpolating function (such as the relationship between the shape and parameters) will be presented in the future research. Moreover, the upper boundary in Theorem 7 would not be the best, so it will also be presented in the future research. Acknowledgements This work was supported by the National Nature Science Foundations of China under Grant Nos. 60873111 and 60933007. We thank to all referees for their comments which helped us to improve the paper.

W.-X. Huang, G.-J. Wang / Applied Mathematics and Computation 217 (2011) 4644–4653

4653

References [1] B. Barsky, D.P. Greenberg, Determining a set of B-spline control vertices to generate an interpolating surface, Comput. Graph. Image Proc. 14 (1980) 203–226. [2] B.K. Choi, H.Y. Shin, W.S. Yoo, Visually smooth composite surfaces for an unevenly spaced:3D data array, Comput. Aided Geom. Des. 10 (1993) 157– 171. [3] J. Peters, Local smooth surface interpolation: a classification, Comput. Aided Geom. Des. 7 (1990) 191–195. [4] Q. Duan, L. Wang, E.H. Twizell, A new bivariate rational interpolation based on function values, Inform. Sci. 166 (2004) 181–191. [5] Q. Duan, K. Djidjeli, W.G. Price, E.H. Twizell, A rational cubic spline based on funcion values, Comput. Graph. 22 (1998) 479–486. [6] Q. Duan, Y. Zhang, E.H. Twizell, A bivariate rational interpolation and the properties, Appl. Math. Comput. 179 (2006) 190–199. [7] Q. Duan, H. Zhang, Y. Zhang, Bounded Property and point control of a bavariate rational interpolating surface, Comput. Math. Appl. 52 (2006) 975–984. [8] Q. Duan, K. Djidjeli, W.G. Price, E.H. Twizell, Weighted rational cubic spline interpolation and its application, J. Comput. Appl. Math. 117 (2000) 121– 135. [9] Q. Duan, L. Wang, E.H. Twizell, A new weighted rational cubic interpolation and its approximation, Appl. Math. Comput. 168 (2005) 990–1003. [10] Q. Duan, L. Wang, E.H. Twizell, A new C2 rational interpolation based on function values and constrained control of the interpolant curves, Appl. Math. Comput. 161 (2005) 311–322. [11] M.P. George, Interpolation and Approximation by Polynomials, Springer-Verlag, 2003.