Rev. Anal. Num´er. Th´eor. Approx., vol. 34 (2005) no. 2, pp. 195–205 ictp.acad.ro/jnaat
BIDIMENSIONAL INTERPOLATION OPERATORS
OF FINITE ELEMENT TYPE AND DEGREE OF EXACTNESS TWO
DANIELA ROS¸CA∗
Abstract. For a given arbitrary triangulation ofR2, we construct an interpo- lating operator which is exact for the polynomials in two variables of total degree
≤2. This operator is local, in the sense that the information around an inter- polation node are taken from a small region around this point. We study the remainder of the interpolation formula.
MSC 2000. 41A63, 41A05, 41A25, 41A80, 47A57.
Keywords. Two-dimensional interpolation operator, degree of exactness.
1. PRELIMINARIES
Given a set of V distinct points in R2,we construct a triangulation T and denote by ∆, ∆⊆R2, the set covered by the triangles ofT. Then, for each of the triangles inT, some functions will be associated in the following way. Con- sider the triangleT ∈ T with the verticesM1(x1, y1), M2(x2, y2), M3(x3, y3) and the number
DT =
x1 y1 1 x2 y2 1 x3 y3 1
.
We define the functionsA, B, C :R2→R,depending on the triangle T, A(x, y) = (x3−x2)(y−y3)−(yD 3−y2)(x−x3)
T
= (x3−x2)(y−y2)−(yD 3−y2)(x−x2)
T ,
B(x, y) = (x1−x3)(y−y1)−(yD 1−y3)(x−x1)
T
= (x1−x3)(y−y3)−(yD 1−y3)(x−x3)
T ,
C(x, y) = (x2−x1)(y−y2)−(yD 2−y1)(x−x2)
T
= (x2−x1)(y−y1)−(yD 2−y1)(x−x1)
T
∗Technical University of Cluj-Napoca, Dept. of Mathematics, Str. Daicoviciu 15, 400020 Cluj-Napoca, Romania, e-mail: [email protected].
and X, Y :R2→R,
X(x, y) =x, Y(x, y) =y, for all (x, y)∈R2. Definition 1. The functions fT
M1, gT
M1, hT
M1 : R2 → R, associated to the vertexM1 of the triangle T, defined by
fMT
1 = 2C3+ 2B3−3C2−3B2+ 1 +α1ABC, gMT
1 = (x3−x1)(C3−CB2−2C2) + (x2−x1)(B3−BC2−2B2) +X−x1+β1ABC,
hTM
1 = (y3−y1)(C3−CB2−2C2) + (y2−y1)(B3−BC2−2B2) +Y −y1+γ1ABC,
are called shape functions. Analogously, for the vertices M2 and M3, the shape functions are defined by circular permutations of the functions A, B, C, the parameters α1, β1, γ1 being replaced with α2, β2, γ2 for the vertex M2 and respectively with α3, β3, γ3 for M3.
Some immediate properties of these functions are given in the following proposition.
Proposition 1. The following statements are true for i, j ∈ {1,2,3}.
1) fT
Mi(Mj) =δij, ∂f
T Mi
∂x (Mj) = 0, ∂f
T Mi
∂y (Mj) = 0, 2) gT
Mi(Mj) = 0, ∂g
T Mi
∂x (Mj) =δij, ∂g
T Mi
∂y (Mj) = 0, 3) hT
Mi(Mj) = 0, ∂h
T Mi
∂x (Mj) = 0, ∂h
T Mi
∂y (Mj) =δij,
4) Along the edges of the triangle M1M2M3, the functions fT
Mi, gT
Mi, hT
Mi
depend only on the corresponding vertices.
We also give some preliminary results, which will be used in the sequel.
These results are summarized in the following proposition.
Proposition2. LetT =M1M2M3 be a triangle and letfMTi, gTMi, hTMi be the shape functions defined in Definition1fori= 1,2,3.The following statements are true.
1)
3
P
i=1
fMTi = 1−ABC(12 +
3
P
i=1
αi), 2) P3
i=1
xifMT
i +gTM
i =X−ABC P3
i=1
βi+αixi+ 4xi, 3)
3
P
i=1
yifMTi+hTMi =Y −ABC
3
P
i=1
γi+αiyi+ 4yi, 4) P3
i=1
x2ifMT
i + 2xigTM
i = X2 −ABCh6(x1x2 +x2x3 +x3x1) + P3
i=1
2βixi +αix2i −2x2ii,
5)
3
P
i=1
y2ifMTi + 2yihTMi = Y2 −ABCh6(y1y2 +y2y3 +y3y1) +
3
P
i=1
2γiyi
+αiy2i −2yi2i, 6) P3
i=1
xiyifMT
i+yigTM
i+xihTM
i =XY −ABCh3(x2y1+x1y2+x3y1+x1y3 +x2y3+x3y1)−2 (x1y1+x2y2+x3y3) +
3
P
i=1
γixi+βiyi+αixiyi
i. Proof. The proof of the above formulas needs some elementary but quite long calculations which are not given here because of lack of space.
2. THE INTERPOLATION OPERATOR
As in [3], to a given function ϕ : R2 → R, ϕ ∈ C1(∆), we associate the piecewise polynomial function Pϕ: ∆→R,
(1) (Pϕ)(x, y) =
V
X
i=1
ϕ(Mi)Fi(x, y) +∂ϕ∂x(Mi)Gi(x, y) +∂ϕ∂y(Mi)Hi(x, y), for (x, y)∈∆,where
Fi(x, y) =
fT
Mi(x, y), if (x, y) is inside or on the edges of the triangleT ∈ T which hasMi as vertex, 0, on the triangles which do not containMi
and analogous definitions for the functionsGi and Hi.
The functions Fi, Gi, Hi are continuous on ∆ due to the property 4) of Proposition 1, so the piecewise polynomial function Pϕ is continuous on ∆.
In [3] we also gave the graphs of some particular functions Fi, Gi.Hi. It is also immediate thatPϕsatisfies the interpolating conditions (2) (Pϕ)(Mi) =ϕ(Mi), ∂(Pϕ)∂x (Mi) = ∂ϕ∂x(Mi), ∂(Pϕ)∂y (Mi) = ∂ϕ∂y(Mi), fori= 1, . . . , V.
Other properties of the interpolating operator P : C1(∆) → C1(∆), de- fined in (1), are given in the following proposition.
Proposition3. The operatorP reproduces the polynomials in two variables of global degree two if and only if the following conditions are satisfied.
3
X
i=1
αi =−12,
3
X
i=1
αixi+βi =−4
3
X
i=1
xi,
3
X
i=1
αiyi+γi =−4
3
X
i=1
yi,
3
X
i=1
αix2i + 2βixi =−6(x1x2+x2x3+x3x1) + 2
3
X
i=1
x2i,
3
X
i=1
αiy2i + 2γiyi =−6(y1y2+y2y3+y3y1) + 2
3
X
i=1
yi2,
3
X
i=1
αixiyi+βiyi+γixi =−3(x2y1+x1y2+x3y1+x1y3+x2y3+x3y1)
+ 2
3
X
i=1
xiyi.
Proof. The proof is immediate if we use Proposition 2 and the fact that the expression A(x, y)B(x, y)C(x, y) cannot be zero for all (x, y) ∈∆. Also, we have to remark that, for a fixed (x, y) ∈ ∆, the sum in (1) contains at most three nonzero terms, namely those corresponding to the triangle which
contains the point (x, y).
As an immediate consequence of this proposition we deduce the following result.
Proposition 4. Let T be a triangulation. If the parameters αi, βi, γi, i= 1,2,3,are solutions of the system of equations given in Proposition 3, then the following identities are true for(x, y)∈∆:
V
X
i=1
Fi(x, y) = 1,
V
X
i=1
(xi−x)Fi(x, y) +Gi(x, y) = 0,
V
X
i=1
(yi−y)Fi(x, y) +Hi(x, y) = 0,
V
X
i=1
(xi−x)2Fi(x, y) + 2(xi−x)Gi(x, y) = 0,
V
X
i=1
(yi−y)2Fi(x, y) + 2(yi−y)Hi(x, y) = 0,
V
X
i=1
(xi−x)(yi−y)Fi(x, y) + (yi−y)Gi+ (xi−x)Hi(x, y) = 0.
The question now is if there exist such operators which reproduces the polynomials of global degree two. More precisely, we have to decide if the system of six equation, given in Proposition 3, with the unknownsαi, βi, γi, i= 1,2,3,is solvable for all xi, yi, i= 1,2,3.
If we make the change of variables αi =λi−4, i= 1,2,3,
β1 =µ1+ 3(x1−x2), β2=µ2+ 3(x2−x3), β3 =µ3+ 3(x3−x1), γ1 =δ1+ 3(y1−y2), γ2=δ2+ 3(y2−y3), γ3 =δ3+ 3(y3−y1), the system becomes
λ1+λ2+λ3= 0,
λ1x1+λ2x2+λ3x3+µ1+µ2+µ3 = 0, λ1y1+λ2y2+λ3y3+δ1+δ2+δ3 = 0,
λ1x21+λ2x22+λ3x23+ 2(µ1x1+µ2x2+µ3x3) = 0, λ1y21+λ2y22+λ3y32+ 2(δ1y1+δ2y2+δ3y3) = 0,
λ1x1y1+λ2x2y2+λ3x3y3+µ1y1+µ2y2+µ3y3+δ1x1+δ2x2+δ3x3 = 0, and is always solvable, being homogeneous. Moreover, this system is undeter- mined.
Remark 1. The same homogeneous system can be obtained making the change of variables
αi =λi−4, i= 1,2,3,
β1 =µ1+ 3(x1−x3), β2=µ2+ 3(x2−x1), β3 =µ3+ 3(x3−x2), γ1 =δ1+ 3(y1−y3), γ2=δ2+ 3(y2−y1), γ3 =δ3+ 3(y3−y2).
In the following we will restrict ourselves to the zero solution of the homo- geneous system, more precisely to the operator P defined in (1), with
α1=α2 =α3 =−4,
β1 = 3(x1−x2), β2 = 3(x2−x3), β3= 3(x3−x1), γ1= 3(y1−y2), γ2= 3(y2−y3), γ3 = 3(y3−y1).
In order to study the rest Rϕ of the interpolation formula ϕ=Pϕ+Rϕ,
we need to establish the following result.
Lemma 1. Let T =M1M2M3, Mi(xi, yi), i = 1,2,3, be a triangle. Then, for all (x, y)∈T, the following inequalities are true.
1) 0≤fMT
1(x, y)≤1, 2)
gMT
1(x, y)≤maxn1627|x2−x1|,1627|x3−x1|o, 3) hTM1(x, y)≤maxn1627|y2−y1|,1627|y3−y1|o. Proof. 1) The inequalities were already proved in [2].
2) Using the equality x−x1−(x3 −x1)C−(x2−x1)B = 0,the function gTM1 can be written as
gMT1 = (1−B−C) (rC(1−C+B) +sB(1−B−2C)), wherer =x3−x1, s=x2−x1.
Case 1. x2=x3.Denoting
φ(B, C) = (1−B−C)(C(1−C) +B(1−B)−BC),
we have to find the extremes ofφwhen 0≤B ≤1, 0≤C≤1, B+C ≤1.The stationary points of the function φ are (B, C) ∈ {(0,1),(1,0),(ρ, ρ)}, where ρ is a root of the equation 9ρ2 −7ρ+ 1 = 0, namely ρ1,2 = (7±√
13)/18.
The inequality B+C≤1 is not satisfied by the stationary point (ρ1, ρ1). At the other stationary points, the functionφtake the values φ(1,0) =φ(0,1) = 0, φ(ρ2, ρ2) =r·(35 + 13√
13)/486.
Then, as in [2],1 we can prove that, on the edges of the triangle T, the functiongMT1 takes values between 0 and 4r27.
In conclusion, in this case we have gMT1(x, y)≤ 1627|r|,for all (x, y)∈T.
Case 2. x2 6= x3. Making the transform ω : R2 → R2 described by the functions
(3) u=u(x, y) =C−B, v=v(x, y) =B+C, (u, v)∈R2, the triangleT maps into a domain denoted byU and the functiongMT
1(x, y) = gTM1 ω−1(u, v)=ψ(u, v) becomes
(4) ψ(u, v) = 12(1−v)r(v+u)(1−u) +s(v−u)(1−3v+u2 ). The stationary points of the function ψare
(5) (u, v)∈ {(1,1),(−1,1),(ρ(s−r), ρ(2r−s) + 1)},
whereρis a root of the equation 3(4s2−8rs+r2)ρ2+ 4(r−2s)ρ+ 1 = 0,that is
ρ1,2 = 4s−2r∓p4s2+ 8rs+r2−1, if 4s2−8rs+r2 6= 0, ρ = (4(2s−r))−1, if 4s2−8rs+r2 = 0.
We have to decide which of the stationary points (u, v) are situated inside the domainU. The stationary points (1,1) and (−1,1) are situated on the border of U and ψ(1,1) =ψ(−1,1) = 0.
Subcase 2a). 4s2 −8rs+r2 6= 0. In this case, at the stationary points Pρ1,2((s−r)ρ1,2,(2r−s)ρ1,2+ 1), B and C have the values
B = v−u2 = 12(1 + (3r−2s)ρ1,2), C = u+v2 = 12(1 +rρ1,2).
1The restrictions to the edges are the same for our function and for the function considered in [2].
The condition that a stationary point Pρ is situated inside the domain U is equivalent to B ∈ [0,1], C ∈[0,1], B+C ≤1. In Appendix A we proved that, ifPρis situated inside U,then
|ψ(Pρ)| ≤maxn1627|r|,1627|s|o, forρ=ρ1,2. On the edges, we have
gT
M1|M
1M2 = s B(1−B)2, and take values between 0 and 4s27, gT
M1|M
1M3 = r C(1−C)2, and take values between 0 and 4r27, gT
M1|M
2M3 = 0.
In conclusion, in this subcase we have
(6) gMT1(x, y)≤maxn16|r|27 ,16|s|27 o.
Subcase 2b). 4s2 −8rs+r2 = 0. In this case ρ = (4(2s−r))−1, s =
1±
√ 3 2
r and at the stationary point Pρ((s−r)ρ,(2r−s)ρ+ 1), B and C take the values
B = 2±
√ 3
8 , C = 7±
√ 3 16 . The only stationary point situated inside U is P3+
√ 3 16 ,11−3
√ 3 16
, where the functionψ take the value ψ(P) = 19+11
√ 3
256 r 'r·0.1486. . . On the edges, |gMT
1|take values between 0 and 1 +
√ 3 2
4
27|r|< 278 |r|.
In conclusion, in this subcase we also have
gMT1(x, y)≤maxn1627|r|,1627|s|o.
3) The proof is analogously with 2).
3. THE INTERPOLATION FORMULA
In this section we study the interpolation formula
(7) ϕ=Pϕ+Rϕ,
withP defined in (1), by proving the following theorem.
Theorem 1. Let ϕ∈C3(int ∆). Then we have (8) kRϕk∞= sup
(x,y)∈∆
|ϕ(x, y)−(Pϕ)(x, y)| ≤329 +√
2K3L3max, where
(9) K3 = sup
(x,y)∈int ∆
n ∂
3ϕ
∂x3(x, y),∂x∂32ϕ∂y(x, y),∂x∂y∂3ϕ2(x, y),∂∂y3ϕ3(x, y)o and Lmax is the length of the greatest edge of the triangles of T.
Proof. We follow the ideas in [3] and write some Taylor formulas around the point (x, y)∈int ∆.For all 1≤i≤V we have
ϕ(Mi) = ϕ(x, y) + (xi−x)∂ϕ∂x(x, y) + (yi−y)∂ϕ∂y(x, y) +2!1 h(xi−x)∂x∂ + (yi−y)∂y∂ i(2)ϕ(x, y) +Ri(x, y),
∂ϕ
∂x(Mi) = ∂ϕ∂x(x, y) + (xi−x)∂∂x2ϕ2(x, y) + (yi−y)∂x∂y∂2ϕ(x, y) +Si(x, y),
∂ϕ
∂y(Mi) = ∂ϕ∂y(x, y) + (xi−x)∂x∂y∂2ϕ(x, y) + (yi−y)∂∂y2ϕ2(x, y) +Ti(x, y), where
Ri(x, y) = 3!1 h(xi−x)∂x∂ + (yi−y)∂y∂i(3)ϕ(c1i, c2i), Si(x, y) = 2!1 h(xi−x)∂x∂ + (yi−y)∂y∂i(2) ∂ϕ∂x(d1i, d2i), Ti(x, y) = 2!1 h(xi−x)∂x∂ + (yi−y)∂y∂i(2) ∂ϕ∂y(e1i, e2i)
with (c1i, c2i),(d1i, d2i),(e1i, e2i) situated ‘between’ the points (x, y) andMi(xi, yi).
Replacing these expressions into (1) we obtain, using the identities given in Proposition 4,
Pϕ(x, y) =ϕ(x, y) +
V
X
i=1
Ri(x, y)Fi(x, y) +Si(x, y)Gi(x, y) +Ti(x, y)Hi(x, y).
On the other hand, for all i= 1, . . . , V we have
|Ri(x, y)| ≤ 16K3(|x−xi|+|y−yi|)3.
If (x, y) is situated inside a triangleT, then the sumPVi=1Ri(x, y)Fi(x, y) has at most three terms which are not zero, namely those corresponding to the vertices of the triangleT. In this case we can write2
V
X
i=1
|Ri(x, y)| ≤
√ 2 3 K3
3
X
i=1
h(x−xτ(i))2+ (y−yτ(i))2i
3
2 ≤√
2K3L3max, where (xτ(i), yτ(i)) are the vertices of the triangle T.
If (x, y) is situated on an edgeE, then at least two terms of this sum are nonzero, namely those which correspond to the end-points of the edge E. In this case we have
V
X
i=1
|Ri(x, y)| ≤ 2
√ 2
3 K3L3max. Further, for all i= 1, . . . , V we have
|Si(x, y)| ≤ 12K3(|x−xi|+|y−yi|)2
2We use the inequality (a+b)2≤2(a2+b2) fora=|x−xτ(i)|, b=|y−yτ(i)|.
and using Lemma 1 we obtain
V
X
i=1
|Si(x, y)| · |Gi(x, y)| ≤ 1627
3
X
i=1
|Sτ(i)(x, y)| · max
j=1,2,3|xτ(i)−xτ(j)|
≤ 1627LmaxK3 3
X
i=1
(x−xτ(i))2+ (y−yτ(i))2
≤ 169K3L3max, and analogously
V
X
i=1
|Ti(x, y)| · |Hi(x, y)| ≤ 169K3L3max,
for (x, y) situated inside a triangleT. When (x, y) is situated on an edge with the end-pointsMk, Mj, we have
Gk(x, y) =Hj(x, y) = 0.
Combining all the above inequalities, we finally obtain that (10) |(Pϕ)(x, y)−ϕ(x, y)| ≤329 +√
2K3L3max,
whence the conclusion.
4. APPENDIX
With the notations of Lemma 1, we have to decide which of the stationary points Pρ1 and Pρ2, ρ1,2 = 4s−2r∓√
4s2+ 8rs+r2−1 are situated inside the domainU.
We consider two cases.
Case 1. s6= 0.Denoting t= rs,condition B ∈[0,1] reduces to
−1≤ 3t−2
4−2t−√
t2+8t+4 ≤1,
for ρ1 whens >0 and forρ2 when s <0 (case A)
−1≤ 3t−2
4−2t+√
t2+8t+4 ≤1,
for ρ1 whens <0 and forρ2 when s >0 (case B) and is fulfilled when
t∈h−4 + 2√
3,0i∪h17−
√ 97
12 ,∞i, in case A t∈h−4 + 2√
3,17+
√97 12
i
, in case B.
ConditionC ∈[0,1] reduces to
−1≤ t
4−2t−√
t2+8t+4 ≤1,in case A
−1≤ t
4−2t+√
t2+8t+4 ≤1,in case B
and is fulfilled when
t∈−∞,−4−2√
3i∪h−4 + 2√ 3,4−
√ 10 2
i∪h34,∞ in case A, t∈−∞,−4−2
√
3i∪h−4 + 2√ 3,4+
√10 2
i
in case B.
Finally, condition B+C≤1 means
2t−1 4−2t−√
t2+8t+4 ≤0, in case A,
2t−1 4−2t+√
t2+8t+4 ≤0, in case B.
and is satisfied for
t∈−∞,−4−2√
3i∪h−4 + 2√
3,12i∪h4−2√
3,∞ in case A, t∈−∞,−4−2√
3i∪h−4 + 2√
3,12i∪h4 + 2√
3,∞ in case B.
Intersecting all the intervals and denoting I1=h−4 + 2√
3,0i, I2 =h34,∞, I3 = h−4 + 2√
3,12i,we find:
Ifs >0, then Pρ1 ∈U fort∈I1∪I2, Pρ2 ∈U fort∈I3. Ifs <0, then Pρ1 ∈U fort∈I3,
Pρ2 ∈U fort∈I1∪I2. Then, the values of ψ at the stationary points are
ψ(ρ(s−r), ρ(2r−s) + 1) =−14(s−2r)2ρρ2(r2−8rs+ 4s2) +ρ(2r−4s) + 1, ρ=ρ1,2.Using the equality ρ2(4s2−8rs+r2) =−13(1 + 4ρ(r−2s)),we find
ψ(Pρ) =−16(s−2r)2ρ[1 +ρ(r−2s)], meaning
ψ(Pρ1,2) = 16(s−2r)(1−2t)(t−2±√
t2+8t+4) (2t−4±√
t2+8t+4)2 . We also need the inequalities
|s−2r| ≤3·max{|r|,|s|}, if rs <0, (11)
|s−2r| ≤2·max{|r|,|s|}, if rs >0.
(12)
Definingξ1 :I1∪I2 →R,
ξ1(t) = (1−2t)(t−2+√
t2+8t+4) (2t−4+√
t2+8t+4)2 , we have
ξ01(t) = 3(8+t−3
√
t2+8t+4) (4−2t−√
t2+8t+4)3 >0 soξ1 will be increasing onI1 and on I2.Since
ξ1(−4 + 2√ 3) =
√3−5
16 , ξ1(0) = 0, ξ134=−169, lim
t→∞ξ1(t) =−49,
using (11) and (12) we finally obtain
(13) |ψ(Pρ1)| ≤maxn1627|r|,1627|s|o. Further, defining ξ2 :I3 →R,
ξ2(t) = (1−2t)(t−2−√
t2+8t+4) (2t−4−√
t2+8t+4)2 , we have
ξ20(x) = 3(8+t+3
√
t2+8t+4) (4−2t+√
t2+8t+4)3 >0 and therefore the function is increasing on I3.Since
ξ2(−4 + 2√ 3) =
√3−5
16 , ξ2(0) =−19, ξ2(12) = 0, using again (11) and (12) we have
(14) |ψ(Pρ2)| ≤ 5−
√3
32 max{|r|,|s|}<maxn1627|r|,1627|s|o.
Case 2. s= 0.In this case ρ is a root of the equation 3r2ρ2+ 4rρ+ 1 = 0, meaning ρ1 = −1/r, ρ2 = −1/(3r). Conditions B ∈ [0,1], C ∈ [0,1] and B+C ∈[0,1] are satisfied only by the stationary point Pρ2 and the value of ψ at this point isψ(Pρ2) = 4r/27,whence|ψ(Pρ2)|< 1627|r|.
REFERENCES
[1] Petrila, T. and Gheorghiu, C.I., Metode element finit ¸si aplicat¸ii, Ed. Academiei, 1987.
[2] Ros¸ca, D.,Bounds of some shape functions, Aut. Comput. Appl. Math.,13, no. 1, pp.
183–189, 2004.
[3] Rosca, D.,Finite element operators for scattered data, in M. Ivan (ed.), Mathematical Analysis and Approximation Theory, Mediamira Science Publisher, pp. 229–238, 2005.
[4] Stancu, D.D., Coman, Gh. and Blaga, P., Analiz˘a numeric˘a ¸si teoria aproxim˘arii, Presa Universitar˘a Clujean˘a, 2002.
[5] Zienkiewicz, O.C. andTaylor, R.L.,The finite element method, vol.2:Solid Mechan- ics, Butterworth-Heinemann, 2000.
Received by the editors: February 9, 2005.