On Gauss-Lobatto integration on the triangle

0.27 MB·English

by  Yuan Xu

ON GAUSS-LOBATTO INTEGRATION ON THE TRIANGLE
YUANXU
1 Abstract. A recent result in [2] on the non-existence of Gauss-Lobatto cu-
1 baturerulesonthetriangleisstrengthenedbyestablishingalowerboundfor
0 the number of nodes of such rules. A method of constructing Lobatto type
2 cubaturerulesonthetriangleisgivenandusedtoconstructseveralexamples.

n
a
J
1. Introduction
5

Recentlyin[6],motivatedbyhp-finiteelementmethod,Gauss-Lobattocubature
]
A rule on the triangle
N (cid:52):={(x,y):0≤x,y,x+y ≤1}
. isstudied,whichrequires(n+1)(n+2)/2nodes,with(n−1)(n−2)/2nodesinthe
h
interior,n−1nodesineachsideand1pointateachvertex,andiscapableofexactly
t
a integrating polynomials of degree 2n−1. We call such a rule strict Gauss-Lobatto.
m
The main result of [6] shows that such rules do not exist. In this note we consider
[ cubature rules that have nodes in both interior and boundary of the triangle, and

establish a lower bound for the number of nodes, from which the non-existence
2
v of the strict Gauss-Lobatto rule follows immediately. We also study the structure
8 of rules that attain our lower bound and give a method for constructing cubature
4 rules of degree 2n−1 with n−1 nodes on each side and 1 node at each vertex.
8
The development is based on the observation that cubature rules with nodes on
4
the boundary can be constructed, by restricting to the class of bubble functions
.
7 (functions that vanish on the boundary), from cubature rules with nodes on the
0
interior. Thisleadstoabootstrappingschemefortransformingsomecubaturerules
0
with interior points into higher order rules with a specific number of nodes on the
1
: boundary. Several examples are constructed to illustrate the algorithm.
v
Xi 2. Results
r For n ∈ N , let Π2 denote the space of polynomials of (total) degree n in two
a 0 n
variables. It is known that
(cid:18) (cid:19)
n+2 (n+1)(n+2)
dimΠ2 = = .
n 2 2
Let W(x,y) be a non-negative weight function on the triangle (cid:52) with finite mo-
ments. AcubatureruleofprecisionswithrespecttoW isafinitesumthatsatisfies
(cid:90) N
(cid:88)
(2.1) f(x,y)W(x,y)dxdy = λ f(x ,y ), ∀f ∈Π2.
k k k s
(cid:52) k=1
Date:January6,2011.
2000 Mathematics Subject Classification. 65D32.
Key words and phrases. triangle,cubature,Gauss-Lobatto.
1
2 YUANXU
WechooseW tobetheJacobiweightW (x,y)=xαyβ(1−x−y)γ forα,β,γ >
α,β,γ
−1. The case α = β = γ = 0 corresponds to the constant weight. These weight
functions are often considered together with the orthogonal polynomials, called
Jacobi polynomials on the triangle, that are orthogonal with respect to them; see,
for example, [5, p. 86], and [1, 2] in connection with cubature rules.
For the Jacobi weight, it is known that the number of nodes N for (2.1) satisfies
(cid:40)
n(n+1) +(cid:98)n(cid:99) if s=2n−1,
(2.2) N ≥ 2 2
n(n+1) if s=2n−2.
2
This lower bound is classical for s = 2n−2 ([8]) and given in [1] for s = 2n−1,
whichagreeswithM¨oller’slowerboundforcentrallysymmetricweightfunctions[7]
as well as [5]. A cubature rule that attains the lower bound is naturally minimal,
meaning that it has the smallest number of nodes among all cubature rules of the
samedegree. Thelowerbound,however,ismostlikelynotsharp;thatis,aminimal
cubature could require more points than what the lower bound indicates. Minimal
cubature rules are sometimes called Gaussian cubature rules. Theirconstruction is
closely related to orthogonal polynomials of several variables. For discussion along
this line, see [3, 9] and references therein.
The cubature rules that we consider are of precision s and are of the form
(cid:90) (cid:88)N0 (cid:88)N1
(2.3) f(x,y)W (x,y)dxdy = λ f(x ,y )+ λ f(x ,0)
α,β,γ k,0 k,0 k,0 k,1 k,1
(cid:52) k=1 k=1
(cid:88)N2 (cid:88)N3
+ λ f(0,y )+ λ f(x ,1−x )
k,2 k,2 k,3 k,3 k,3
k=1 k=1
+µ f(0,0)+µ f(1,0)+µ f(0,1),
0 1 2
where (x ,y ) are distinct points in the interior of (cid:52), (x ,0), (0,y ), and
k,0 k,0 k,1 k,2
(x ,1−x ) are distinct points on the side y =0, x=0, and x+y =1 (but not
k,3 k,3
on the corners) of (cid:52), respectively, λ >0 and µ >0. Such a cubature has
k,j i
N :=N +N +N +N +3
0 1 2 3
nodes. The main result in [6] states that such a cubature rule does not exist if
s=2n−1 and
(n−2)(n−1)
N = , N =N =N =n−1,
0 2 1 2 3
which has a total number of nodes N = (n + 2)(n + 1)/2. This follows as an
immediate corollary of the following theorem.
Theorem 2.1. If a cubature rule of the form (2.3) exists with precision s=2n−1
or s=2n, then
(cid:40)
n(n−1) if s=2n−1,
(2.4) N ≥ 2
0 n(n−1) +(cid:98)n−1(cid:99) if s=2n,
2 2
(cid:40)
n(n−1) +(cid:98)n−1(cid:99) if s=2n−1,
(2.5) N +N ≥ 2 2 i=1,2,3.
0 i n(n+1) if s=2n,
2
GAUSS-LOBATTO INTEGRATION ON THE TRIANGLE 3
Proof. The cubature (2.3) exactly integrates degree s polynomials of the form
f(x,y)=xy(1−x−y)g(x,y) if
(cid:90) (cid:88)N0
(2.6) g(x,y)W (x,y)dxdy = λ∗ g(x ,y ), ∀g ∈Π2 ,
α+1,β+1,γ+1 k,0 k,0 k,0 s−3
(cid:52) k=1
where λ∗ = λ x y (1−x −y ), which is a cubature rule of precision
k,0 k,0 k,0 k,0 k,0 k,0
s−3 for the weight function W so that, by (2.2), N has to satisfy
α+1,β+1,γ+1 0
the lower bound in the inequality of (2.4). On the other hand, the cubature (2.3)
exactly integrates degree s polynomials of the form f(x,y)=x(1−x−y)g(x,y) if
∀g ∈Π2 ,
s−3
(cid:90) g(x,y)W (x,y)dxdy =(cid:88)N0 λ˜ g(x ,y )+(cid:88)N1 λ∗ g(x ,0)
α+1,β,γ+1 k,0 k,0 k,0 k,1 k,1
∆ k=1 k=1
where λ˜ = λ x (1−x −y ) and λ∗ = λ x (1−x ), which is a
k,0 k,0 k,0 k,0 k,0 k,1 k,1 k,1 k,1
cubature rule of precision s−2 for the weight function W (x,y), so that
α+1,β,γ+1
N +N satisfies the lower bound in (2.2), which gives the inequality of (2.5) for
0 1
i=1. Similarly, we can derive lower bound for N +N and N +N . (cid:3)
0 2 0 3
One naturally asks if there is any cubature rule that attains the lower bound in
thetheorem. Fors=2n−1, thisasksifthereisacubatureruleofprecision2n−1
with
(cid:22) (cid:23)
n(n−1) n−1
(2.7) N = and N = , i=1,2,3.
0 2 i 2
Weexpectthattheanswerisnegative. Aheuristicargumentcanbegivenasfollows:
Assume that N = n(n−1). Then the proof of the theorem shows that (2.6) is a
0 2
cubature of degree 2n−4 with n(n−1)/2 nodes, which is known to exist only for
small n. Assume that it does exist. We define a linear functional L , acting on
1
polynomials of one variable, by
(cid:90) (cid:88)N0
(2.8) L g := g(x)W (x,y)dxdy− λ x (1−x −y )g(x ),
1 α+1,β,γ+1 k,0 k,0 k,0 k,0 k,0
(cid:52) k=1
wherex , y andλ areasin(2.3). Applying (2.3)onpolynomialsoftheform
k,0 k,0 k,0
f(x,y)=g(x)x(1−x−y) shows that
(cid:88)N1
(2.9) L g = λ∗ g(x ), ∀g ∈Π ,
1 k,1 k,1 2n−3
k=1
whereλ∗ =λ x (1−x ). ThefunctionalL definedinequation(2.8)defines
k,1 k,1 k,1 k,1 1
a bilinear form [p,q] = L (pq) that could be indefinite ([q,q] is not necessarily
1 1 1
positive). If the bilinear form were positive definite on Π , then (2.9) could
2n−3
be regarded as a quadrature rule of N nodes and of degree 2n−3 for L and,
1 1
consequently, N ≥ n − 1 by the standard result in Gaussian quadrature rule,
1
which is stronger than the second equation of (2.7). Thus, in order for (2.7) to
hold, we would need L to be indefinite on Π , that is, L(q2) = 0 for some
1 2n−3
nonzero q ∈ Π , and we would need to require that L has a quadrature rule
n−1 1
of degree 2n−3 with N = (cid:98)n−1(cid:99) nodes. A simple count of variables (the nodes
1 2
and weights of the quadrature rule) and restraints (the polynomials that need to
be exactly integrated) shows that this is unlikely to happen, although it still might
4 YUANXU
as the equations are nonlinear. For a linear functional that defines an indefinite
bilinearform,thetheoryofGaussianquadraturerulebreaksdownsinceorthogonal
polynomials may not exist and, even they do, they may not have real or simple
zeros. In particular, we do not have a lower bound for the number of nodes of a
quadrature rule for such a linear functional.
The above argument indicates that if we want the cubature rule (2.3) to have
the smallest number of interior points, then the best that we can hope for will be
a cubature rule of degree 2n−1 that satisfies
n(n−1)
(2.10) N = , and N ≥n−1, 1≤i≤3.
0 2 i
Weshallcallacubaturerulethatattainsthelowerboundin(2.10)Gauss-Lobatto
cubature rule. Theproofofthelowerboundindicateshowsuchacubaturerulecan
be constructed. Let us also define linear functionals
(cid:90) (cid:88)N0
L g := g(y)W (x,y)dxdy− λ y (1−x −y )g(x ),
2 α,β+1,γ+1 k,0 k,0 k,0 k,0 k,0
(cid:52) k=1
(cid:90) (cid:88)N0
L g := g(x)W (x,y)dxdy− λ x y g(x ).
3 α+1,β+1,γ k,0 k,0 k,0 k,0
(cid:52) k=1
Then we can summarize the method of construction as follows.
Algorithm. We can follow the following procedure to construct a Gauss-Lobatto
cubature rule of form (2.3):
Step 1. Construct a cubature rule of degree s−3 for W in the form
α+1,β+1,γ+1
of (2.6) with all nodes in the interior of (cid:52), and define
λ∗
(2.11) λ = k,0 .
k,0 (x y (1−x −y ))
k,0 k,0 k,0 k,0
Step 2. Construct Gaussian quadrature rules (2.9) and
(cid:88)N2 (cid:88)N3
L g = λ∗ g(y ), L g = λ∗ g(x ), ∀g ∈Π ,
2 k,2 k,2 3 k,3 k,3 2n−3
k=1 k=1
with respect to the linear functional L ,L ,L , which gives the nodes of the cuba-
1 2 3
ture rule (2.3) on the boundary of the triangle, and define
λ∗ λ∗ λ∗
(2.12) λ = k,1 , λ = k,2 , λ = k,3 .
k,1 x (1−x ) k,2 y (1−y ) k,3 x (1−x )
k,1 k,1 k,2 k,2 k,3 k,3
Step 3. Finally, the weight µ , µ , and µ are determined by setting f(x,y) =
0 1 2
1,x,y in (2.3) and solve the resulted linear system of equations. (cid:3)
Proposition2.2. Assume,intheabovealgorithm,thatthelinearfunctionalL ,L
1 2
and L are positive definite on Π and all x ,y ,x are inside (0,1). Then
3 2n−3 k,1 k,2 k,3
the algorithm produces a cubature rule of degree 2n − 1 in the form (2.3). In
particular, if N = n(n−1), then the cubature rule is Gauss-Lobatto.
0 2
Proof. We need to show that the cubature rule constructed by the algorithm holds
for all f ∈Π2 , that is, (2.3) holds for all f ∈Π2 with s=2n−1. A moment of
2n−1 s
reflection shows that, with z =1−x−y, Π2 can be decomposed into a direct sum.
s
(2.13) Π2 =xyzΠ2 ⊕xzΠ [x]⊕yzΠ [y]⊕xyΠ [x]⊕Π2,
s s−3 s−2 s−2 s−2 1
GAUSS-LOBATTO INTEGRATION ON THE TRIANGLE 5
where Π [x] and Π [y] denote the space of polynomials of one variable in x-
s−2 s−2
variable and y-variable, respectively. If f ∈ xyzΠ2 , then (2.3) reduces to (2.6),
s−3
which holds by our construction. If f ∈ xzΠ [x], then (2.3) reduces to (2.9)
s−2
for L , which holds by our construction. The same holds for f ∈ yzΠ [y] and
1 s−2
f ∈ xyΠ [x], whereas for f ∈ Π2, the cubature is verified by Step 3. Thus, by
s−2 1
(2.13), the cubature rule holds for all f ∈Π . (cid:3)
2n−1
Itshouldbepointedoutthat,aslongaswehaveacubatureruleofdegree2n−4
with all nodes in the interior of (cid:52) in the Step 1, regardless if it is a minimal one,
the Step 2 and Step 3 could be carried out and Proposition 2.2 applies. Thus, the
algorithmcanbeusedtoconstructcubaturerulesofdegree2n−1intheform(2.3)
with N =N =N =n−1.
1 2 3
Let us comment on how the steps in the algorithm can be realized.
For Step 1, a cubature with the specification can be constructed by solving
moment equations, that is, solving the system of equations formed by setting
g(x,y) = xiyj for i+j ≤ s−3 in (2.6) for x , y and λ∗ . There have been
k,0 k,0 k,0
a number of papers based on this method, see e.g. [11] for the latest result and
further references. Another approach is to use a characterization of the cubature
rules that attain lower bound in (2.4), which is given in terms of common zeros
of certain orthogonal polynomials and can be used to find cubature rules of lower
order, see [7, 9]. Not all cubature rules obtained via either methods work for our
purpose, since we require that all nodes are inside the domain.
For Step 2, one needs to check that L ,L ,L are positive linear functionals on
1 2 3
Π . Oncetheyare,thestandardprocedureofconstructingGaussianquadrature
2n−3
rules applies. In particular, we can apply the standard algorithm to generate a
sequenceoforthogonalpolynomialsuptodegreen−1withrespecttoL inductively;
i
the nodes of the Guassian quadrature rule of degree 2n−3 for L are the zeros of
i
the orthogonal polynomial of degree n−1 with respect to L .
i
In the following we give several examples of Gauss-Lobatto cubature rules of
degree 5 and 7 for the unit weight function. The minimal cubature rules for the
degree 5 and 7 have nodes 7 and 12, respectively ([4]). For Gauss-Lobatto rules,
the number of nodes are necessarily larger, as seen in (2.7).
Example 1. There is a cubature formula of degree 5 with 12 nodes, 3 interior, 2
on each side and 1 at each corner of the triangle. To illustrate our procedure, we
shall present the nodes and weights in steps.
Thethreeinteriorpointsandweightsaregiveninthefirsttable. Whilethenodes
arethoseofacubatureruleofdegree2fortheweightfunctionW (x,y)=xy(1−
1,1,1
x−y), the weights are relates to those of the latter cubature rule by (2.11). These
nodes are common zeros of quasi-orthogonal polynomials (see [10] for definition) of
degree2whichwerefoundbysolvingthenonlinearsystemofequationsinTheorem
4.1 of [10].
x y λ
k,0 k,0 k,0
0.15881702219143 0.19201873632215 0.101342396527698
0.56219234596964 0.19201873632215 0.117181247909596
0.22100936816107 0.55798126367785 0.118066904793533
6 YUANXU
where the nodes are the common zeros of following three polynomials of degree 2:
√ √ √
19− 105+x(−91+ 105+112x)+2(−7+ 105)y,
√ √ √
49−3 105−175x+5 105x+112x2−84y+4 105y+224xy,
√ √ √
154−6 105−301x+11 105x+112x2−609y+7 105y+560xy+560y2.
The nodes on the edges of the triangle, but not on the corners, are the nodes of
GaussianquadraturerulesforL ,L ,L ,respectively,wheretheweightsarerelated
1 2 3
to those of Gaussian quadrature rules by (2.12). They are given in the following
three tables:
x λ
k,1 k,1
0.3931870086016 0.02991955921794
0.8595419130359 0.01756588222187
where the nodes are zeros of the orthogonal polynomials p of degree 2 for L ,
1 1
1 (cid:16) √ (cid:17) 1 (cid:16) √ (cid:17)
p (x)=x2+ −469−9 105 x+ 889+61 105 ;
1 448 4480
y λ
k,2 k,2
0.4305843026985 0.02290932968619
0.7924406473476 0.02022650113138
where the nodes are zeros of the orthogonal polynomials p of degree 2 for L ,
2 2
3 (cid:16) √ (cid:17) 3 (cid:16) √ (cid:17)
p (x)=x2+ −29+ 105 x+ 63+ 105 ;
2 46 644
x λ
k,3 k,3
0.2629899118578 0.02514330117112
0.7030163143652 0.03109870484395
where the nodes are zeros of the orthogonal polynomials p of degree 2 for L ,
3 3
1 (cid:16) √ (cid:17) 1 (cid:16) √ (cid:17)
p (x)=x2+ −10997+51 105 x+ 665−9 105 .
3 10843 3098
Finally, the weights µ ,µ ,µ are given by
1 2 3
µ µ µ
1 2 3
0.0081170837035 0.00326155091683 0.00516753787639 (cid:3)
The nodes of this Lobatto cubature rule are depicted in Figure 1.
Our next two examples are Gauss-Lobatto cubature rules that are symmetric
in the sense that the nodes are invariant under the symmetric group of (cid:52), or the
permutation of (x,y,1−x−y). A symmetric Gauss-Lobatto rule for the constant
weight takes the form
(cid:90) (cid:88)M0
(2.14) f(x,y)dxdy = A [f(u ,v )+f(v ,w )+f(w ,u )]
k k,0 k,0 k,0 k,0 k,0 k,0
(cid:52) k=1
(cid:88)M1
+ B [f(u ,0)+f(0,1−u )+f(1−u ,u )]
k k,1 k,1 k,1 k
k=1
+C[f(0,0)+f(1,0)+f(0,1)],
where w =1−u −v .
k,0 k,0 k,0
GAUSS-LOBATTO INTEGRATION ON THE TRIANGLE 7
Figure 1. Nodes of Gauss-Lobatto cubature rules of degree 5
Both examples are constructed by following the steps in the algorithm. The
symmetrymakestheconstructionmucheasier,sinceweonlyneedtoconsiderpoly-
nomials that are symmetric under the symmetric group of (cid:52). In particular, the
cubatureruleforW (x,y)inStep1canbefoundbysolvingthereducedmoment
1,1,1
equations of symmetric polynomials. We shall skip details and only list the nodes
and weights of these two cubature rules as formulated in (2.14). Their nodes are
depicted in Figure 2.
Example 2. Symmetric Gauss-Lobatto cubature rules of degree 5 with 12 nodes.
This formula is in the form of (2.14) with M = 1 and M = 2. The nodes and
0 1
weights are given below:
√ √
u =v = 1 (7− 7), A = 7 (cid:0)14− 7(cid:1),
1,0 1,0 21 1 720
u = 1 (cid:16)21−(cid:113)21(cid:0)4√7−7(cid:1)(cid:17), u = 1 (cid:16)21+(cid:113)21(cid:0)4√7−7(cid:1)(cid:17),
1,1 42 1,2 42
√ √
B =B = 1 (cid:0)7+4 7(cid:1), C = 1 (cid:0)8− 7(cid:1).
1 2 720 720
(cid:3)
Example 3. Symmetric Gauss-Lobatto cubature rules of degree 7 with 18 nodes.
This formula is in the form of (2.14) with M = 2 and M = 3. The nodes and
0 1
weights are given below:
√ √
u =v = 1 (5− 7), u =v = 1 (5+ 7),
1,0 1,0 18 2,0 2,0 18
√ √
A = 1 (cid:0)1141−94 7(cid:1), A = 1 (cid:0)1141+94 7(cid:1),
1 17640 2 17640
√ √
u = 1(cid:0)3− 3(cid:1), u = 1, u = 1(cid:0)3+ 3(cid:1),
1,1 6 2,1 2 3,1 6
B = 3 , B = 4 , B = 3 , C = 1 .
1 280 2 315 3 280 315
(cid:3)
These cubature rules appear to be new (see the list in [4]). Their numbers of
nodes are more than the minimal given in (2.2). The existence of higher order
Gauss-Lobatto rules depends on the existence of minimal cubature rules of degree
2n−4 for the weight W . The latter cubature rules most likely do not exist for
1,1,1
n ≥ 6, but the algorithm can still be applied to produce cubature formulas with
n−1 points in each side of the triangle.
8 YUANXU
Figure 2. Nodes of symmetric Gauss-Lobatto cubature rules of
degree 5 and 7
It should be mentioned that there are cubature rules of degree 2n−1 in the
form of (2.3) that have fewer than n−1 points on each side, naturally with more
interior points according to (2.5), such rules cannot be constructed directly by our
algorithm. It is possible to modify the algorithm, however, since L for such rules
i
cannot be positive definite on Π but it must be positive definite on a subspace
n−1
of Π .
n−1
Finally, let us mention that our algorithm can also be modified for constructing
cubature rules of degree 2n. This means a cubature rule of degree 2n − 3 for
W in Step 1, and a quadrature for L of degree 2n−2 in Step 2. The
α+1,β+1,γ+1 i
quadrature of degree 2n−2 is generated by a quasi-orthogonal polynomial of the
formq :=p +αp ,whereαisafreeparameterwhichcanbefixedbyrequiring,
n n n−1
say, q (1/2) = 0, which means fixing the middle point on the corresponding side
n
of the triangle as a node of the cubature. We have tried this construction for
cubature rules of degree 6 with 4 interior points, 3 on each sides and 1 at each
vertex, starting with a cubature rule for W of degree 3. The Lobatto type
1,1,1
cubature rule of degree 6 that we obtained, however, has three negative weights.
Acknowledgment. The author thanks both referees for their careful reading of
the manuscript, especially David Day for his extensive and thoughtful suggestions.
References
[1] H. Berens and H. J. Schmid, On the number of nodes of odd degree cubature formulae for
integrals with Jacobi weight on a simplex, Numerical Integration (Bergen, 1991), 37–44,
NATOAdv.Sci.Inst.Ser.CMath.Phys.Sci.,357,KluwerAcad.Publ.,Dordrecht,1992.
[2] H. Berens, H. Schmid and Y. Xu, On two-dimensional definite orthogonal systems and on
lowerboundforthenumberofassociatedcubatureformulas,SIAMJ.Math.Anal.26(1995),
468–487.
[3] R.Cools,I.P.Mysovskikh,andH.J.Schmid,Cubatureformulaeandorthogonalpolynomials,
J. Comp. Appl, Math127(2001)121-152.
[4] R. Cools and P. Rabinowitz, Monomial cubature rules since “Stroud”: a compilation. J.
Comput. Appl. Math.48(1993),no.3,309–326.
[5] C.F.DunklandY.Xu,OrthogonalPolynomialsofSeveralVariables EncyclopediaofMath-
ematicsanditsApplications81,CambridgeUniversityPress,Cambridge,2001.
[6] B.T.Helenbrook,Ontheexistenceofexplicithp-finiteelementmethodusingGauss-Lobatto
integrationonthetriangle,SIAM J. Numer. Anal.47(2009),1304-1318.
GAUSS-LOBATTO INTEGRATION ON THE TRIANGLE 9
[7] H.Mo¨ller,KubaturformelnmitminimalerKnotenzahl,Numer. Math. 25(1976),185–200.
[8] A. H. Stroud, Approximate calculation of multiple integrals, Prentice-Hall, Inc., Englewood
Cliffs,N.J.,1971.
[9] Y.Xu,Commonzerosofpolynomialsinseveralvariablesandhigherdimensionalquadrature,
PitmanResearchNotesinMathematicsSeries,Longman,Essex,1994.
[10] Y.Xu,Onzerosofmultivariatequasi-orthogonalpolynomialsandGaussiancubatureformu-
lae,SIAM J. Math. Anal.25(1994),991-1001.
[11] L. Zhang, T. Cui, and H. Liu, A set of symmetric quadratures on triangles and tetrahedra,
J. Computational Math.,27(2009),89-96.
Department of Mathematics, University of Oregon, Eugene, Oregon 97403-1222.
E-mail address: [email protected]